diag.nb <- function(modelo=fit.model,iden=c(0,0,0,0,0,0,0,0),nome=seq(along = model.matrix(modelo)[,1]),res="D") { # # Descrição e detalhes: # A saída terá oito gráficos: # 1º) Influência na Locação. O gráfico feito é das distâncias de Cook contra os valores ajustados. Utilizou-se o critério # de destacar observações maiores do que duas vezes a média de todas as distâncias obtidas; # 2º) Influência Locação/Escala. A medida C, que é um aperfeiçoamento do DFFIT e também é conhecida como distância de Cook # modificada, foi utilizada para medir a influência das observações nos parâmetros de locação e escala. O critério # foi o de destacar observações maiores do que duas vezes a média de todas as distâncias obtidas; # 3º) Influência Local. 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. 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. Também não estamos # considerando o parâmetro theta, apenas os parâmetros do preditor linear. 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; # 4º) Função de Ligação. O gráfico feito é do preditor linear ajustado (eta) contra a variável dependente ajustada (z). # Segundo McCullagh e Nelder (1989, pág.401), o padrão esperado é de uma linha reta. Para funções de ligação da # família potência uma curvatura dos pontos acima da reta sugere uma ligação de uma potência maior que a utilizada, # enquanto que uma curvatura abaixo da reta uma potência menor. Conforme sugerido, é adicionado uma reta suavizada # pelo método lowess robusto e também uma linha tracejada com inclinação de 45º. O método deve ser utilizado com # cautela, uma vez que por ex.para a binomial ele não é informativo; # 5º) Pontos Alavanca 1. Para os MLG's a matriz H=sqrt(W)%*%X%*%solve(t(X)%*%W%*%X)%*%t(X)%*%sqrt(W) é interpretada como # sendo a matriz de projeção da solução de mínimos quadrados da regressão linear de z contra X com pesos W. Assim, # sugere-se utilizar h=diag(H) para detectar a presença de pontos alavanca no modelo de regressão normal ponderado. Um # ponto é considerado alavanca (leverage) quando este exerce uma forte influência no seu valor ajustado. O critério foi # o de destacar observações maiores do que duas vezes a média de todos os h’s, que nesse caso resume-se a duas vezes # o número de parâmetros do modelo, dividido pelo número de observações. # A medida de alavanca h depende da ligação através dos pesos. Como o gráfico de alavanca é feito dos h's contra as # médias, temos uma idéia do que esperar dependendo da ligação escolhida. # Na ligação canônica (log), o peso é mu/(mu/theta+1). # Na ligação raiz quadrada (sqrt) o peso é 4/(mu/theta+1). # Na ligação identidade (mu) o peso é o inverso da média 1/((mu^2)/theta+mu). # Os pesos da ligação log são baixos para valores pequenos de mu e crescem conforme mu aumenta. O crescimento # é suave para valores pequenos de theta e são mais "bruscos" conforme theta aumenta. Os pesos da ligação raiz e # identidade são altos para pequenos valores de mu e diminuem conforme mu aumenta. A diferença é que o caimento # para a ligação identidade é mais suave que o da raiz quadrada. # Caso se queira ter uma idéia visual dos pesos, use os comandos abaixo para desenhar as curvas para 4 valores de # theta e mu de 1 a 10: # mu<-seq(1,10,0.1) # par(mfrow=c(2,2)) # for(t in c(.5,1,5,10)) { # w1<-mu/(mu/t+1) # w2<-4/(mu/t+1) # w3<-1/((mu^2)/t+mu) # plot(cbind(mu,mu,mu),cbind(w1,w2,w3),type="n",xlab="mu",ylab="peso") # lines(mu,w1,col=1,lty=1,lwd=3) # lines(mu,w2,col=3,lty=1,lwd=3) # lines(mu,w3,col=4,lty=1,lwd=3) # legend(x=min(mu),y=max(c(w1,w2,w3)),legend=c("log","raiz","ident."),col=c(1,3,4),lty=c(1,1,1),lwd=c(3,3,3)) # title(paste("theta=",t)) # } # 6º) Pontos Alavanca 2. Comentamos no quinto gráfico que a medida h tem uma forte dependência dos pesos conforme a # ligação escolhida. Assim, sugerimos uma medida h modificada que se dá por hm=diag(HM), em que # HM=X%*%solve(t(X)%*%W%*%X)%*%t(X). A idéia por trás dela é de tentar eliminar a forte dependência com os pesos e # facilitar a detecção dos pontos alavanca. Obviamente a medida não elimina toda a dependência, uma vez que HM ainda # depende de W, mas podemos utilizá-la adicionalmente à medida h já tradicionalmente sugerida. Vale notar que se for # escolhida uma ligação em que W seja uma matriz identidade, então h=hm. Hosmer e Lemeshow (2000, págs.169 a 173) # discutem as duas medidas para a regressão logística. Parece não existir estudos muito aprofundados de hm, mas # Hosmer e Lemeshow discutem críticas feitas por Pregibon à medida h com relação a dependência de W e críticas feitas # por Lesaffre em tentar ignorar a informação contida em W. Por fim, sugerem que ambas sejam utilizadas com cautela; # 7º) Pontos Aberrantes. Um ponto é aberrante (discrepante, outlier) se o seu valor estiver mal ajustado pelo modelo. # Adicionamos linhas tracejadas em -2 e 2. Assim, se os resíduos escolhidos forem aproximadamente normais, # esperamos que cerca de 5% dos pontos possam estar um pouco fora desses limites. Mesmo sem normalidade, o gráfico # serve para inspecionar valores mal ajustados. Deve-se no entando tomar cuidado pois sabemos que por ex.o resíduo # de Pearson é assimétrico. Esse gráfico serve como indicação para detectar valores aberrantes marginalmente. # Devido o desconhecimento da distribuição dos resíduos e se o objetivo for detectar valores conjuntamente aberrantes # deve-se construir o gráfico de envelopes; # 8º) Função de Variância. McCullagh e Nelder (1989, pág.400) sugere o gráfico dos resíduos absolutos contra os valores # ajustados (ou contra os valores ajustados transformados em escala constante) para checar se a função de variância # adotada é adequada. O padrão esperado é de não encontrarmos nenhuma tendência. Funções de variância erradas irão # resultar em tendências dos resíduos com a média. Tendências positivas indicam que a função de variância está # crescendo muito devagar com a média, então deve-se aumentar a potência (no caso de uma função de variância da # família potência). Uma linha suavizada pelo método lowess robusto é adicionada para ajudar na procura de tendências. # # Os dados devem estar disponíveis pelo comando attach( ). # # Argumentos obrigatórios: # modelo: deve-se informar o objeto onde está o ajuste do modelo com distribuição binomial negativa, 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 8 # posições de números inteiros. A ordem que deve ser informada é a mesma em que os gráficos são feitos. 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 respectivos permitir que identifiquemos o número de pontos solicitados e qualquer # outro valor (negativo ou decimal) parar nos gráficos e solicitar que especifiquemos o número de pontos a ser # destacado. O padrão é c(0,0,0,0,0,0,0,0) caso não se entre com nada e c(-1,-1,-1,-1,-1,-1,-1,-1) caso se entre # com qualquer coisa que não seja um vetor de 8 posições, como por ex.-1; # 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 (índices). # 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; # res: permite-se a escolha dos resíduos que serão utilizados nos gráficos de pontos aberrantes, da função # de ligação e na medida de influência na locação e na escala. As opções dos resíduos são: "Q" quantil (ver Dunn e # Smyth, 1996), "D" componente do desvio e "P" Pearson padronizado. A opção padrão é a "D". # # A função retorna os seguintes valores: ResQuantil, ResCompDesv, ResPearsonStd, Di, Ci, Dmax e h. # # Autor: Frederico Zanqueta Poleto , arquivo disponível em http://www.poleto.com # # Referências: # DUNN, K. P., and SMYTH, G. K. (1996). Randomized quantile residuals. J. Comput. Graph. Statist. 5, 1-10 # [http://www.statsci.org/smyth/pubs/residual.html e http://www.statsci.org/smyth/pubs/residual.ps] # HOSMER, D. W. e LEMESHOW, S. (2000). Applied Logistic Regression. John Wiley & Sons, New York. # MCCULLAGH, P. e NELDER, J. A. (1989). Generalized Linear Models. 2ª ed. Chapman and Hall, London. # 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. # # Exemplos: # diag.nb(ajuste,iden=c(1,5,2,0,4,3,0,0),nome=estados) # diag.nb(ajuste,iden=-1) # if(class(modelo)[1] != "negbin") { stop(paste("\nA classe do objeto deveria ser binomial negativa e nao ",class(modelo)[1],"!!!\n")) } if(length(iden)<8) { iden<-c(-1,-1,-1,-1,-1,-1,-1,-1) } X <- model.matrix(modelo) n <- nrow(X) p <- ncol(X) fi <- modelo$theta w <- modelo$weights W <- diag(w) Fis <- t(X)%*%W%*%X V <- solve(Fis) H <- sqrt(W)%*%X%*%V%*%t(X)%*%sqrt(W) h <- diag(H) HM <- X%*%V%*%t(X) hm <- diag(HM) #para evitar divisão por 0 ao studentizar os residuos, mas tentando manter o valor exagerado da alavanca h[round(h,15)==1]<-0.999999999999999 y <- modelo$y m <- predict(modelo,type="response") pl <- predict(modelo) adj <- pl+residuals(modelo,type="working") #variável dependente ajustada rp <- resid(modelo,type="pearson") ts <- rp/sqrt(1-h) td <- resid(modelo,type="deviance")/sqrt(1-h) pi<-fi/(m+fi) a<-ifelse( y>0 , pbeta(pi,fi,pmax(y,1)) , 0) b<-pbeta(pi,fi,y+1) rq<-qnorm( runif(n=n,min=a,max=b) ) di <- (h/((1-h)*p))*(ts^2) A <- diag(rp)%*%H%*%diag(rp) dmax <- abs(eigen(A)$vec[,1]/sqrt(eigen(A)$val[1])) if(res=="Q") { cat("Ao utilizar o residuo Quantil para distribuicoes discretas, sugere-se plotar pelo menos 4 graficos para evitar conclusoes viesadas pela aleatoriedade que esta sendo incluida.\n") tipor<-"Resíduo Quantil" r<-rq } else { if(res=="D") { tipor<-"Resíduo Componente do Desvio" r<-td } else { if(res=="P") { tipor<-"Resíduo de Pearson Padronizado" r<-ts } else { stop(paste("\nVoce nao escolheu corretamente um dos residuos disponiveis!!!\n")) } } } ci <- sqrt( ((n-p)*h) / (p*(1-h)) )*abs(r) par(mfrow=c(2,4)) plot(m,di,xlab="Valor Ajustado", ylab="Distância de Cook",main="Influência na Locação", ylim=c(0,max(di,2*mean(di))), pch=16) abline(2*mean(di),0,lty=2) while ( (!is.numeric(iden[1])) || (round(iden[1],0) != iden[1]) || (iden[1] < 0) ) { cat("Digite o num.de pontos a ser identificado (0=nenhum) e para continuar\n") out <- readline() iden[1]<-as.numeric(out) } if(iden[1]>0) {identify(m,di,n=iden[1],labels=nome)} plot(m,ci,xlab="Valor Ajustado", ylab="Distância de Cook Modificada",main="Influência Locação/Escala", ylim=c(0,max(ci,2*mean(ci))), pch=16) abline(2*mean(ci),0,lty=2) while ( (!is.numeric(iden[2])) || (round(iden[2],0) != iden[2]) || (iden[2] < 0) ) { cat("Digite o num.de pontos a ser identificado (0=nenhum) e para continuar\n") out <- readline() iden[2]<-as.numeric(out) } if(iden[2]>0) {identify(m,ci,n=iden[2],labels=nome)} plot(m,dmax,xlab="Valor Ajustado", ylab="dmax",main="Influência Local", ylim=c(0,max(dmax,2*mean(dmax))), pch=16) abline(2*mean(dmax),0,lty=2) while ( (!is.numeric(iden[3])) || (round(iden[3],0) != iden[3]) || (iden[3] < 0) ) { cat("Digite o num.de pontos a ser identificado (0=nenhum) e para continuar\n") out <- readline() iden[3]<-as.numeric(out) } if(iden[3]>0) {identify(m,dmax,n=iden[3],labels=nome)} plot(adj,pl,xlab="Variável Dependente Ajustada",ylab="Preditor Linear Ajustado",main="Função de Ligação", pch=16) lines(lowess(adj,pl)) abline(a=0,b=1,lty=2) while ( (!is.numeric(iden[4])) || (round(iden[4],0) != iden[4]) || (iden[4] < 0) ) { cat("Digite o num.de pontos a ser identificado (0=nenhum) e para continuar\n") out <- readline() iden[4]<-as.numeric(out) } if(iden[4]>0) {identify(adj,pl,n=iden[4],labels=nome)} plot(m,h,xlab="Valor Ajustado",ylab="Medida h",main="Pontos Alavanca 1",ylim=c(0,max(h,2*p/n)),pch=16) abline(2*p/n,0,lty=2) while ( (!is.numeric(iden[5])) || (round(iden[5],0) != iden[5]) || (iden[5] < 0) ) { cat("Digite o num.de pontos a ser identificado (0=nenhum) e para continuar\n") out <- readline() iden[5]<-as.numeric(out) } if(iden[5]>0) {identify(m,h,n=iden[5],labels=nome)} plot(m,hm,xlab="Valor Ajustado",ylab="Medida h Modificada",main="Pontos Alavanca 2",ylim=c(0,max(hm,2*mean(hm))),pch=16) abline(2*mean(hm),0,lty=2) while ( (!is.numeric(iden[6])) || (round(iden[6],0) != iden[6]) || (iden[6] < 0) ) { cat("Digite o num.de pontos a ser identificado (0=nenhum) e para continuar\n") out <- readline() iden[6]<-as.numeric(out) } if(iden[6]>0) {identify(m,hm,n=iden[6],labels=nome)} plot(m,r,xlab="Valor Ajustado",ylab=tipor,main="Pontos Aberrantes", ylim=c(min(r)-1,max(r)+1), pch=16) abline(2,0,lty=2) abline(-2,0,lty=2) while ( (!is.numeric(iden[7])) || (round(iden[7],0) != iden[7]) || (iden[7] < 0) ) { cat("Digite o num.de pontos a ser identificado (0=nenhum) e para continuar\n") out <- readline() iden[7]<-as.numeric(out) } if(iden[7]>0) {identify(m,r,n=iden[7],labels=nome)} plot(m,abs(r),xlab="Valor Ajustado",ylab=paste(tipor," Absoluto",sep=""),main="Função de Variância", pch=16) lines(lowess(m,abs(r))) while ( (!is.numeric(iden[8])) || (round(iden[8],0) != iden[8]) || (iden[8] < 0) ) { cat("Digite o num.de pontos a ser identificado (0=nenhum) e para continuar\n") out <- readline() iden[8]<-as.numeric(out) } if(iden[8]>0) {identify(m,abs(r),n=iden[8],labels=nome)} par(mfrow=c(1,1)) list(ResQuantil=rq,ResCompDesv=td,ResPearsonStd=ts,Di=di,Ci=ci,Dmax=dmax,h=h) }