dmax.glm <- function(modelo=fit.model,iden=double(ncol(model.matrix(modelo))),nome=seq(along = model.matrix(modelo)[,1]),coefs=NA) { # # Descrição e detalhes: # A saída terá gráficos do dmax para cada coeficiente de um ajuste de um modelo linear generalizado ou binomial negativa. # # A influência local consiste em procurar pontos que sob pequenas perturbações causam variações muito grandes nos # resultados. O dmax é o autovetor que corresponde ao maior autovalor da matriz do processo de perturbações. Para # maiores detalhes veja Paula (2003, págs.50-54 e 65-66). O critério foi o de destacar observações maiores do que duas # vezes a média de todos os dmax's absolutos. Na influência local utiliza-se a matriz de informação de Fisher observada. # Como estamos utilizando a matriz esperada, os resultados obtidos são apenas aproximados para MLG's que não sejam ajustados # com a ligação canônica e para a binomial negativa. Na binomial negativa também não estamos considerando o parâmetro # theta, apenas os parâmetros do preditor linear. No caso da binomial negativa o correto seria programar as expressões # obtidas por Svletiza (2002, pág.23) que são as da matriz observada e também possui o termo theta. # # Os dados devem estar disponíveis pelo comando attach( ). # # Argumentos obrigatórios: # modelo: deve-se informar o objeto onde está o ajuste de um modelo linear generalizado, caso não seja informado, a # função procurará o ajuste no objeto fit.model; # # Argumentos opcionais: # iden: caso deseje, informe o número de observações que irá querer destacar em cada gráfico. O vetor deve # conter números inteiros. A ordem que deve ser informada é a mesma das variáveis da opção var, caso seja # utilizada, ou deve ter a mesma ordem da matriz modelo. Os componentes do vetor iguais a 0 indicam que não se # quer que identifique pontos, se for um inteiro positivo irá automaticamente nos gráficos pertinentes permitir # que identifiquemos o número de pontos solicitados e qualquer outro valor (negativo ou decimal) parará nos # gráficos e irá solicitar que especifiquemos o número de pontos a ser destacado. O padrão é c(0,...,0) caso não # se entre com nada e c(-1,...,-1) caso se entre com qualquer coisa que não satisfaça os requisitos citados # de ser número inteiro, não negativo e de ter o mesmo comprimento da opção var ou da matriz modelo; # nome: esse argumento só é utilizado caso algum dos componentes do vetor da opção iden não seja 0. Caso não # seja informado nada, os pontos identificados serão os números da ordem em que estão no banco de dados. # Caso se queira, pode-se informar um vetor de nomes ou de identificações alternativas. Obrigatoriamente # esse vetor deve ter o mesmo comprimento do banco de dados; # coefs: se nada for informado para essa opção, será feito um gráfico para cada coeficiente do modelo. Se o seu modelo # tiver muitos coeficientes o que faria com que a função fizesse muitos gráficos pequenos ou se não quiser # o gráfico para todos os coeficientes utilize essa opção para especificar todos os coeficientes que você quer que # sejam avaliados. Por ex., utilizando 2:9 fará com que a função avalie os primeiros 8 coeficientes com exceção # do intercepto. A função irá utilizar a mesma ordem impressa no summary do ajuste do modelo. # # Autor: Frederico Zanqueta Poleto , arquivo disponível em http://www.poleto.com # # Referências: # PAULA, G. A. (2003). Modelos de Regressão com apoio computacional. IME-USP, São Paulo. [Não publicado, # disponível em http://www.ime.usp.br/~giapaula/Book.pdf] # SVETLIZA, C. F. (2002). Modelos Não-Lineares com Resposta Binomial Negativa. Tese de Doutorado. IME-USP, São Paulo. # # Exemplo: # dmax.glm(ajuste,iden=-1,nome=estados) # if(class(modelo)[1] != "glm" & class(modelo)[1] != "negbin") { stop(paste("\nA classe do objeto deveria ser glm ou negbin e nao ",class(modelo),"!!!\n")) } X <- model.matrix(modelo) n <- nrow(X) p <- ncol(X) r<-resid(modelo,type="pearson") w <- modelo$weights W <- diag(w) if(is.na(coefs)[1]) { coefs<-1:p } p <- length(coefs) if(p>length(iden)) { iden<-rep(-1,p) } if (p>2) { if (p>8) { par(mfrow=c(3,ceiling(p/3))) } else { par(mfrow=c(2,ceiling(p/2))) } } else { par(mfrow=c(1,ceiling(p))) } j<-0 for(i in coefs) { j<-j+1 expl<-"" X1<-X[,i] X2<-X[,-i] Fis <- t(X2)%*%W%*%X2 V <- solve(Fis) H <- sqrt(W)%*%X2%*%V%*%t(X2)%*%sqrt(W) v<-sqrt(W)%*%X1-H%*%sqrt(W)%*%X1 dmax<-v*r dmax<-dmax/sqrt(as.double(t(dmax) %*% dmax)) #cut<-2/sqrt(n) cut<-2*mean(abs(dmax)) plot(dmax,xlab="Índice", ylab="dmax",main=dimnames(X)[[2]][i],ylim=c(min(-cut,dmax),max(cut,dmax)), pch=16) abline(-cut,0,lty=2) abline(cut,0,lty=2) while ( (!is.numeric(iden[j])) || (round(iden[j],0) != iden[j]) || (iden[j] < 0) ) { cat("Digite o num.de pontos a ser identificado (0=nenhum) e para continuar\n") out <- readline() iden[j]<-as.numeric(out) } if(iden[j]>0) {identify(dmax,n=iden[j],labels=nome)} } par(mfrow=c(1,1)) cat("\n") }