#========================================= # ANÁLISE DE DADOS CATEGÓRICOS # Profa.: SUELY RUIZ GIOLO # Prof.: Enrico A. Colosimo # Segunda Parte: Modelos: log-lineares, # e Poisson #========================================= # # 1. Modelos Log-lineares - três Dimensões #========================================== # # 1.0 Exemplo Ilustrativo # require(vcd) # Independência Mútua teste10<-array(c(2,2,2,2,2,2,2,2),c(2,2,2)) mosaic(teste10) teste11<-array(c(3,1,3,1,3,1,3,1),c(2,2,2)) teste11 mosaic(teste11) # Independência Marginal (AB,C) teste20<-array(c(1,3,3,1,1,3,3,1),c(2,2,2)) teste20 mosaic(teste20) teste21<-array(c(1,3,1,3,3,1,3,1),c(2,2,2)) teste21 mosaic(teste21) # Independência Condicional teste30<-array(c(3,1,3,1,3,3,1,1),c(2,2,2)) teste30 mosaic(teste30) # # # 1.1 Exemplo: Alcool, Cigarro e Maconha # # dados <-array(c(911,44,3,2,538,456,43,279),c(2,2,2)) dimnames(dados) <- list(maconha = c("Sim", "Não"), cigarro = c("Sim", "Não"), alcool= c("Sim", "Não")) dados mosaic(dados) # dados1<-matrix(dados[c(1:2,5:6)],2) prop.table(dados1,1) # alcool<-factor(rep(rep(c("Sim","Nao"),c(2,2)),2)) cigarro<-factor(rep(c("Sim","Nao"),2)) maconha<-factor(rep(c("Sim","Nao"),c(4,4))) freq <-c(911,44,3,2,538,456,43,279) alunos<-data.frame(alcool,cigarro,maconha,freq) alunos # # alcool cigarro maconha freq # 1 Sim Sim Sim 911 # 2 Sim Nao Sim 44 # 3 Nao Sim Sim 3 # 4 Nao Nao Sim 2 # 5 Sim Sim Nao 538 # 6 Sim Nao Nao 456 # 7 Nao Sim Nao 43 # 8 Nao Nao Nao 279 # # Modelo (A,C,M) glm.out1<-glm(formula=freq ~alcool + cigarro + maconha, family=poisson, data=alunos) summary(glm.out1) glm.out1$fitted p1<-1-pchisq(glm.out1$deviance,glm.out1$df.residual) pearson1<-summary.lm(glm.out1)$residuals sum(pearson1^2) # anova(glm.out1, test="Chisq") anova(glm.out1, test="LRT") # glm.out2<-glm(formula=freq ~alcool*cigarro*maconha, family=poisson, data=alunos) anova(glm.out2, test="LRT") # # Modelo (AC,M) glm.out3<-update(glm.out1, .~.+alcool:cigarro) summary(glm.out3) # # Modelo (AC,AM,CM) glm.out3<-update(glm.out1, .~.+alcool:cigarro + alcool:maconha + cigarro:maconha) summary(glm.out3) # Resíduos --> os padronizados devem ser preferidos pearson<-summary.lm(glm.out3)$residuals # Pearson residuos sum(pearson^2) std.resid<-rstandard(glm.out3,type="pearson") # residuo padronizado expected<-fitted(glm.out3) cbind(freq,expected,pearson, std.resid) # # freq expected pearson std.resid # 1 911 910.38317 0.02044342 0.6333249 # 2 44 44.61683 -0.09234564 -0.6333249 # 3 3 3.61683 -0.32434090 -0.6333250 # 4 2 1.38317 0.52447895 0.6333250 # 5 538 538.61683 -0.02657821 -0.6333249 # 6 456 455.38317 0.02890528 0.6333249 # 7 43 42.38317 0.09474777 0.6333249 # 8 279 279.61683 -0.03688791 -0.6333249 # # O resíduo std.resid deve ser preferido por ser aproximado N(0,1) # Conclusão: Os resultados sugerem que o modelo com todas as interações de # segunda ordem mas sem a de terceira ordem é o melhor modelo. Isto significa # que nenhum par de variáveis é condicionalmente independente dada a terceira # mas que a associação entre duas variáveis é idêntica no dois níveis da # terceira. # #============================================================== # # 2. Modelo de Poisson #============================================================== # 2.1 Exemplo 1: Estudo: Câncer não-melanoma de pele # # Duas cidades (Mineapolis e Dallas) em 1994 # no. de câncer de pele não melanoma # por faixa etaria (8) # 15-24, 25-34, .....,>84 # controlado pela população (offset) #========================================================= # Entrada e Preliminares dos Dados # contagem<-c(1,16,30,71,102,130,133,40,4,38,119,221,259,310,226,65) cidade<-c(rep(0,8),rep(1,8)) faixa<-c(rep(c(20,30,40,50,60,70,80,90),2)) pop<-c(172675,123065,96216,92051,72159,54722,32185,8328,181343, 146207,121374,111353,83004,55932,29007,7538) # #=========================================================== # Descrição dos Dados # # table(cidade) table(cidade,faixa) y<-log(contagem/pop) plot(faixa[cidade==0],y[cidade==0],ylim=c(-13,-4), main="Minneapolis (preto) vs Dallas (vermelho)", xlab="idade", ylab="log(contagem/pop)") # Minneapolis points(faixa[cidade==1],y[cidade==1],col=2) # Dallas # plot(faixa[cidade==1],y[cidade==1],main="Dallas") # Dallas # #================================================================= # Modelo de Poisson # require(hnp) m1<-glm(contagem~cidade + offset(log(pop)),family=poisson) m1 summary(m1) 1-pchisq(m1$deviance,m1$df.residual) # pearson<-summary.lm(m1)$residuals # Pearson residuos x2<-sum(pearson^2) 1-pchisq(x2,m1$df.residual) plot(m1) # Incluindo Faixa Etária Contínua m2<-glm(contagem~cidade+faixa + offset(log(pop)),family=poisson) m2 summary(m2) 1-pchisq(m2$deviance,m2$df.residual) # pearson<-summary.lm(m2)$residuals # Pearson residuos x2<-sum(pearson^2) 1-pchisq(x2,m2$df.residual) # hnp(m2, sim=19, conf=1) plot(m2) # Incluindo Faixa Etária Contínua com termo quadrático m3<-glm(contagem~cidade*faixa + I(faixa^2)+offset(log(pop)),family=poisson) m3 summary(m3) 1-pchisq(m3$deviance,m3$df.residual) # pearson<-summary.lm(m3)$residuals # Pearson residuos x2<-sum(pearson^2) 1-pchisq(x2,m2$df.residual) # hnp(m3, sim=19, conf=1) plot(m3) # # Incluindo Faixa Etária Categórica m4<-glm(contagem~cidade+factor(faixa) + offset(log(pop)),family=poisson) m4 summary(m4) 1-pchisq(m4$deviance,m4$df.residual) # pearson<-summary.lm(m4)$residuals # Pearson residuos x2<-sum(pearson^2) 1-pchisq(x2,m3$df.residual) # plot(m3) # # Q-QPLOT com envelope simulado (código: envel.bino) # hnp(m4, sim=19, conf=1) # # #========================================================================== # Modelo de Poisson/Binomial Negativo # Exemplo 2: Partos Cesarianos # Comparação de Hospitais privados e publicos # Tamanho do hospital (offset) # #========================================================= # Entrada e Preliminares dos Dados # partos<-c(236,739,970,2371,309,679,26,1272,3246,1904,357,1080, 1027,28,2507,138,502,1501,2750,192) hosp<-c(0,rep(1,5),0,rep(1,6),0,1,0,rep(1,4)) ces<-c(8,16,15,23,5,13,4,19,33,19,10,16,22,2,22,2,18,21,24,9) # pos<-seq(1:20) tx<-ces/partos plot(pos[hosp==0],log(tx)[hosp==0],ylim=c(-5,-1),xlim=c(0,20), main="Taxa de Cesarianas",xlab="Hosp: Pub. (vermelho) e Priv. (preto)",ylab="log(Taxa)") points(pos[hosp==1],log(tx)[hosp==1],col=2,pch=16) # #========================================================== # Explorando Regressão Linear SImples # mod1<-lm(ces~partos+hosp) # dois regressores plot(mod1) summary(mod1) # Absurdo o coef. para hospital é + # mod2<-lm(sqrt(ces)~partos+hosp) # dois regressores plot(mod2) summary(mod2) # o coef. para hosp é menor mas continua + # mod3<-lm(ces/partos~hosp) # um regressor plot(mod3) summary(mod3) # agora é negativo # #================================================================= # Modelo de Poisson # m1<-glm(ces~hosp+offset(log(partos)),family=poisson) summary(m1) # o ajuste náo é adequado mas o coef. é consistente e - 1-pchisq(m1$deviance,m1$df.residual) # pearson<-summary.lm(m1)$residuals # Pearson residuos x2<-sum(pearson^2) 1-pchisq(x2,m1$df.residual) # plot(m1) plot(residuals(m1)~fitted(m1)) abline(h=0) # problema com a homocedasticidade hnp(m1, sim=19, conf=1) # plot(log(fitted(m1)),log((ces-fitted(m1))^2), xlab=expression(hat(mu)),ylab=expression((y-hat(mu))^2)) abline(0,1) ## gráfico para verificar se a média = variância # Incluindo um parâmetro de dispersão (dp<-sum(residuals(m1,type="pearson")^2)/m1$df.res) summary(m1,dispersion=dp) # Teste F é mais indicado do que o Qui-quadrado drop1(m1,test="F") # Outra alternativa: variância robusta (estimador sanduiche) library(sandwich) library(lmtest) coeftest(m1, vcov=sandwich) # #=============================================================== # Q-QPLOT com envelope simulado (código: envel.bino) # #============================================================== # require(hnp) hnp(m1, sim=19, conf=1) # #============================================================================ # Modelo Binomial Negativo require(MASS) modn<-glm.nb(ces~offset(log(partos))+hosp) summary(modn) 1-pchisq(modn$deviance,modn$df.residual) # pearson<-summary.lm(modn)$residuals # Pearson residuos x2<-sum(pearson^2) 1-pchisq(x2,modn$df.residual) # hnp(modn, sim=19, conf=1) plot(modn) #=================================================================== m2<-glm(ces~hosp+partos,family=poisson) summary(m2) # o ajuste náo é adequado mas o coef. é consistente e - plot(m2) plot(residuals(m2)~fitted(m2)) abline(h=0) # problema com a homocedasticidade # ======================================================================