#------------------------------------------------------------- # DIAGNÓSTICO GEE PARA REGRESSÃO DE POISSON #------------------------------------------------------------- # Com esse ajuste, calculamos algumas medidas para análise de # diagnóstico e construímos os gráficos conforme # Tan, Qu e Kutner (1997) e Venezuela (2003). #------------------------------------------------------------- # # Descrição dos argumentos da função: # # modelo: objeto do tipo gee(Y~X, family=poisson(link=log), # ...). # # dados: banco de dados utilizado no ajuste do modelo. # # umaJanela: argumento lógico que dá a opção # de imprimir todos os gráficos da análise de diagnóstico # numa mesma janela ou um gráfico em cada janela. Por # default, essa opção é TRUE. # # selGraficos: vetor de uns e zeros que dá a opção de não # imprimir todos os gráficos da análise de diagnóstico. # Para não imprimir um gráfico dentre os quatro indicados, # basta colocar 0 na posição correspondente do gráfico. # Alteração na configuração default c(1,1,1,1) deve ser # utilizada junto com umaJanela=FALSE. # # identifica: vetor de quantidades de observações a serem # identificadas em cada gráfico da análise de diagnóstico. # Por exemplo, para identificar um ponto no 2º gráfico e 1 # no 4º gráfico, basta fazer identifica=c(0,1,0,1). # #------------------------------------------------------------- diag_gee_poisson <- function(modelo, dados, umaJanela=T, selGraficos = c(1,1,1,1), identifica=c(0,0,0,0)) { X <- model.matrix(as.formula(paste("~ ", modelo$call$formula[3])), dados) y <- modelo$y beta <- coef(modelo) R <- modelo$work mi <- fitted(modelo) individuo <- modelo$id repet <- dim(modelo$work)[1] ue <- modelo$nobs/repet N <- nrow(X) p <- ncol(X) #Matriz C <- A * Delta #Ligação canônica -> Delta=Identidade A <- diag(mi,N) C <- A #Matriz Omega - variância e covariância de y Omega <- matrix(0,N,N) invOmega <- matrix(0,N,N) l <- 1 while (l0) identify(individuo,h,labels=labelsGraf,n=identifica[1]) #if (!umaJanela) savePlot("graf1.jpeg", type="jpeg") } if (selGraficos[2]) { plot(hue,xlab="Unidade Experimental", ylab="H por Unidade Experimental", pch=16) abline(cut,0,lty=2) if (identifica[2]>0) identify(hue, n=identifica[2]) #if (!umaJanela) savePlot("graf2.jpeg", type="jpeg") } if (selGraficos[3]) { plot(individuo,cd,xlab="Unidade Experimental", ylab="Distância de Cook", pch=16) if (identifica[3]>0) identify(individuo,labels=labelsGraf,cd,n=identifica[3]) #if (!umaJanela) savePlot("graf3.jpeg", type="jpeg") } if (selGraficos[4]) { plot(individuo,rsd,xlab="Unidade Experimental", ylab="Residuo Padronizado", pch=16) if (identifica[4]>0) identify(individuo,labels=labelsGraf,rsd,n=identifica[4]) #if (!umaJanela) savePlot("graf4.jpeg", type="jpeg") } saida = list(ResPadronizado=rsd, DistCook=cd, h=h, hue=hue) return(invisible(saida)) }