diag.mult <- function(modelo=fit.model,iden=NULL,nome=seq(along = model.matrix(modelo)[,1]),sub=F) { # # Descrição e detalhes: # A maioria das medidas utilizadas foram discutidas por Silva (2003). # 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. Quando sub=F, a medida # é composta pela média dos h's de cada submodelo e pelo resíduo componente do desvio. Essa proposta não é comentada por # Silva (2003) e desconhecemos se já foi utilizada por alguém anteriormente. Quando sub=T, será utilizado a distância de # Cook das observações em cada sub-modelo separadamente. Serão destacadas 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. Quando sub=F, # a medida é composta pela média dos h's de cada submodelo e pelo resíduo componente do desvio. Essa proposta não é # comentada por Silva (2003) e desconhecemos se já foi utilizada por alguém anteriormente. Quando sub=T, a medida # para cada sub-modelo é construída com o h e o resíduo de Pearson do próprio sub-modelo. 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 avaliando a maior curvatura. O dmax é o autovetor que corresponde ao maior autovalor # da matriz do processo de perturbações (Silva, 2003, utiliza o mesmo conceito do dmax descrito, mas chamando-o por # lmax). Para maiores detalhes veja Paula (2003, págs.50-54 e 65-66). Quando sub=F, o dmax é construído utilizando o # modelo inteiro e quando sub=T, cada sub-modelo separadamente. O critério foi o de destacar observações maiores do # que duas vezes a média de todos os dmax's; # 4º) Influência Local Total. A influência local total avalia a curvatura na direção de cada observação, ao invés de avaliar # a maior curvatura. Quando sub=F, a medida é em relação ao modelo inteiro e quando sub=T, as medidas serão obtidas de # cada sub-modelo. O critério foi o de destacar observações maiores do que duas vezes a média de todos os valores # obtidos. Nos referiremos a essa medida por tmax; # 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. Quando sub=F, será utilizado o h médio dos # sub-modelos e quando sub=T, será utilizado o h de cada sub-modelo separadamente. # 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. # 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 h2=diag(H2), em que # H2=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 H2 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=h2. 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 h2, 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. # Quando sub=F, serão utilizadas as medidas h2 médias de todos os sub-modelos e quando sub=T, serão utilizadas as # medidas h de cada sub-modelo; # 7º) Pontos Alavanca 3. Aqui será utilizada a medida de alavanca proposta por Silva (2003, pág.41). Quando sub=F, será # tomada a média das medidas de todos os sub-modelos. Essa medida será denotada por medida m; # 8º) 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. Quando sub=F, será utilizado o resíduo componente do desvio e quando sub=T, # será utilizado o resíduo de Pearson. # # 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 multinomial feito pela biblioteca # nnet, função multinom, 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. Se sub=T, deve-se especificar um vetor # com 8x(J-1) posições, sendo as 8 primeiras últimas referente aos gráficos do 1º sub-modelo e assim por diante; # 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; # sub: se for utilizada a opção padrão (sub=F, False), serão utilizadas medidas globais do ajuste do modelo. Em sub=T, # será feito um gráfico de diagnósticos para cada sub-modelo, parando nos intervalos para caso se queira copiar # o gráfico para algum programa. Para saber as diferenças nas medidas em sub=T e sub=F, veja a seção descrição e # detalhes. # # A função retorna os seguintes valores: Res, Di, Ci, Dmax, Tmax, h, h2 e m. # # Autor: Frederico Zanqueta Poleto , arquivo disponível em http://www.poleto.com # # Referências: # 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] # SILVA, J. A. P. (2003). Métodos de Diagnóstico em Modelos Trinomiais. Dissertação de Mestrado. IME-USP, São Paulo. # # Exemplos: # diag.mult(ajuste,iden=c(1,5,2,0,4,3,0,0),nome=estados) # diag.mult(ajuste,iden=-1,sub=T) # if(class(modelo)[1] != "multinom") { stop(paste("\nA classe do objeto deveria ser multinom (biblioteca nnet) e nao ",class(modelo),"!!!\n")) } X <- model.matrix(modelo) n <- nrow(X) p <- ncol(X) ntot<-modelo$weights mu<-modelo$fitted.values y<-do.call("as.array",list(as.name(dimnames(attr(modelo$terms,"factors"))[[1]][1]))) niveis<-modelo$lab J<-length(niveis) if(sub==F) { if(is.null(iden)) { iden<-rep(0,8) } else { if(length(iden)<8) { iden<-rep(-1,8) } } } else { if(is.null(iden)) { iden<-rep(0,(J-1)*8) } else { if(length(iden)<((J-1)*8)) { iden<-rep(-1,(J-1)*8) } } } if(is.na(dim(y)[2])) { y2<-matrix(0,n,J) for(j in 1:J) { y2[y==niveis[j],j]<-1 } y<-y2 } #IO<-modelo$Hess #if(is.null(IO)) { # modelo<-multinom(y~X-1,Hess=T,trace=F) # IO<-modelo$Hess #} #VarO<-solve(IO) S<-kronecker(diag(J-1),X) V<-matrix(0,(J-1)*n,(J-1)*n) for(j in 1:(J-1)) { for(k in 1:(J-1)) { if(j == k) { V[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)]<-diag(mu[,j]*(1-mu[,j])) } else { V[(n*(j-1)+1):(n*j),(n*(k-1)+1):(n*k)]<-diag(-mu[,j]*mu[,k]) } } } IF<-t(S)%*%V%*%S VarF<-solve(IF) M<-S%*%VarF%*%t(S)%*%V H<-matrix(0,(J-1)*n,(J-1)*n) h<-matrix(0,n,J-1) h2<-matrix(0,n,J-1) m<-matrix(0,n,J-1) for(j in 1:(J-1)) { VT<-V[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)] H[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)]<-sqrt(VT)%*%X%*%solve(t(X)%*%VT%*%X)%*%t(X)%*%sqrt(VT) h[,j]<-diag(H[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)]) h2[,j]<-diag( X%*%solve(t(X)%*%VT%*%X)%*%t(X) ) m[,j]<-diag( M[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)] ) } hmed<-apply(h,1,mean) h2med<-apply(h,1,mean) mmed<-apply(m,1,mean) if(sub==F) { r<-numeric(n) for(j in 1:J) { r<-r+ifelse( y[,j]==0,0,y[,j]*log(y[,j]/(ntot*mu[,j])) ) } r<-sqrt(2*r/(1-hmed)) di<-(hmed/(1-hmed))*(r^2) ci <- sqrt( ((n-p)*hmed) / (p*(1-hmed)) )*abs(r) u<-matrix(0,(J-1)*p,n) for(j in 2:J) { u[((j-2)*p+1):((j-1)*p),]<-t(X)%*%diag(as.numeric( y[,j]-ntot*mu[,j] )) } A<-t(u)%*%VarF%*%u dmax<-abs(eigen(A)$vec[,1]/sqrt(eigen(A)$val[1])) tmax<-2*abs(diag(A)) } else { r<-matrix(0,n,J-1) dmax<-matrix(0,n,J-1) tmax<-matrix(0,n,J-1) di<-matrix(0,n,J-1) ci<-matrix(0,n,J-1) for(j in 1:(J-1)) { r[,j]<-(y[,j]-ntot*mu[,j])/sqrt(ntot*mu[,j]*(1-mu[,j])) A<-diag(r[,j])%*%H[(n*(j-1)+1):(n*j),(n*(j-1)+1):(n*j)]%*%diag(r[,j]) dmax[,j]<-abs(eigen(A)$vec[,1]/sqrt(eigen(A)$val[1])) tmax[,j]<-2*abs(diag(A)) r[,j]<-r[,j]/sqrt(1-h[,j]) di[,j]<-(h[,j]/((1-h[,j])*p))*(r[,j]^2) ci[,j]<-sqrt( ((n-p)*h[,j]) / (p*(1-h[,j])) )*abs(r[,j]) } } par(mfrow=c(2,4)) if(sub==F) { plot(di,xlab="Índice", 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(di,n=iden[1],labels=nome)} plot(ci,xlab="Índice", 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(ci,n=iden[2],labels=nome)} plot(dmax,xlab="Índice", 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(dmax,n=iden[3],labels=nome)} plot(tmax,xlab="Índice", ylab="tmax",main="Influência Local Total", ylim=c(0,max(tmax,2*mean(tmax))), pch=16) abline(2*mean(tmax),0,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(tmax,n=iden[4],labels=nome)} plot(hmed,xlab="Índice",ylab="Medida h",main="Pontos Alavanca 1",ylim=c(0,max(hmed,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(hmed,n=iden[5],labels=nome)} plot(h2med,xlab="Índice",ylab="Medida h Modificada",main="Pontos Alavanca 2",ylim=c(0,max(h2med,2*mean(h2med))),pch=16) abline(2*mean(h2med),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(h2med,n=iden[6],labels=nome)} plot(mmed,xlab="Índice",ylab="Medida m",main="Pontos Alavanca 3",ylim=c(0,max(mmed,2*p/n)),pch=16) abline(2*p/n,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(mmed,n=iden[7],labels=nome)} plot(r,xlab="Índice",ylab="Resíduo Componente do Desvio",main="Pontos Aberrantes", ylim=c(0,max(r)+1), pch=16) abline(2,0,lty=2) 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(r,n=iden[8],labels=nome)} } else { iden2<-iden for(j in 1:(J-1)) { iden<-iden2[((j-1)*8+1):(j*8)] plot(mu[,(j+1)],di[,j],xlab="Valor Ajustado", ylab="Distância de Cook",main="Influência na Locação", ylim=c(0,max(di[,j],2*mean(di[,j]))), pch=16) abline(2*mean(di[,j]),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(mu[,(j+1)],di[,j],n=iden[1],labels=nome)} plot(mu[,(j+1)],ci[,j],xlab="Valor Ajustado", ylab="Distância de Cook Modificada",main="Influência Locação/Escala", ylim=c(0,max(ci[,j],2*mean(ci[,j]))), pch=16) abline(2*mean(ci[,j]),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(mu[,(j+1)],ci[,j],n=iden[2],labels=nome)} plot(mu[,(j+1)],dmax[,j],xlab="Valor Ajustado", ylab="dmax",main="Influência Local", ylim=c(0,max(dmax[,j],2*mean(dmax[,j]))), pch=16) abline(2*mean(dmax[,j]),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(mu[,(j+1)],dmax[,j],n=iden[3],labels=nome)} plot(mu[,(j+1)],tmax[,j],xlab="Valor Ajustado", ylab="tmax",main="Influência Local Total", ylim=c(0,max(tmax[,j],2*mean(tmax[,j]))), pch=16) abline(2*mean(tmax[,j]),0,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(mu[,(j+1)],tmax[,j],n=iden[4],labels=nome)} plot(mu[,(j+1)],h[,j],xlab="Valor Ajustado",ylab="Medida h",main="Pontos Alavanca 1",ylim=c(0,max(h[,j],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(mu[,(j+1)],h[,j],n=iden[5],labels=nome)} plot(mu[,(j+1)],h2[,j],xlab="Valor Ajustado",ylab="Medida h Modificada",main="Pontos Alavanca 2",ylim=c(0,max(h2[,j],2*mean(h2[,j]))),pch=16) abline(2*mean(h2[,j]),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(mu[,(j+1)],h2[,j],n=iden[6],labels=nome)} plot(mu[,(j+1)],m[,j],xlab="Valor Ajustado",ylab="Medida m",main="Pontos Alavanca 3",ylim=c(0,max(m[,j],2*p/n)),pch=16) abline(2*p/n,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(mu[,(j+1)],m[,j],n=iden[7],labels=nome)} plot(mu[,(j+1)],r[,j],xlab="Valor Ajustado",ylab="Resíduo de Pearson Padronizado",main="Pontos Aberrantes", ylim=c(min(r[,j])-1,max(r[,j])+1), pch=16) abline(2,0,lty=2) abline(-2,0,lty=2) 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(mu[,(j+1)],r[,j],n=iden[8],labels=nome)} if (j < (J-1)) { cat(paste("Sub-modelo ",niveis[j+1]," OK. Tecle para o proximo sub-modelo.\n",sep="")) out <- readline() } else { cat(paste("Sub-modelo ",niveis[j+1]," OK.\n",sep="")) } } } par(mfrow=c(1,1)) list(Res=r,Di=di,Ci=ci,Dmax=dmax,Tmax=tmax,h=h,h2=h2,m=m) }