#========================================= # ANÁLISE DE DADOS CATEGÓRICOS # Profa.: SUELY RUIZ GIOLO # Prof.: Enrico A. Colosimo # PRIMEIRA PARTE: Modelo Binomial e # Tabelas de Contingência #========================================= # #========================================= # 1- Modelo Binomial #========================================= # # # 1a. Explorando a Distribuição Qui-quadrado #========================================= # help(dchisq) x<-seq(1,30,.1) plot(x,x*0,pch="",ylim=range(c(0,0.2)),xlim=range(c(0,30)), xlab="x",ylab="f(x)",bty="n") lines(x,dchisq(x,1),lty=1) lines(x,dchisq(x,3),lty=2) lines(x,dchisq(x,5),lty=3) lines(x,dchisq(x,8),lty=4) lines(x,dchisq(x,10),lty=5) lines(x,dchisq(x,15),lty=6) pchisq(3.84,1) 1-pchisq(3.84,1) qchisq(0.95,1) gera<-sort(rchisq(50000,3)) hist(gera,breaks=80,freq=F) lines(x,dchisq(x,3),lty=2) # # # 1b. Explorando a Distribuição Binomial # Testes para uma única amostra: exato # qui-quadrado, MC e Bootstrap #========================================= # help(Binomial) x<-0:20 prob<-dbinom(x,20,0.5) round(prob,4) plot(x,prob,type="h") # # Teste Exato - Uma Amostra # Exemplo: n=20, x=12 # binom.test(12,20) # testa H_0:p=0.5 para n=20 e x=12 # Achar P(8 < X < 12) for X Binomial(20,0.5) pvalor<-1-sum(dbinom(9:11, 20, 0.5)) pvalor # # Teste qui-quadrado # ?chisq.test dados<-c(12,8) chisq.test(dados,correct=F) # # Intervalo de Confiança baseado no teste qui-quadrado # x<-seq(0.05,0.995,0.005) x2<-(((12-(20*x))^2)/(20*x))+((((20*x)-12)^2)/(20*(1-x))) dados1<-cbind(x,x2) plot(x,x2) dados1 # Inf=0.382 e sup=0.782 # # Teste via Simulação MC # chisq.test(dados,correct=F,simulate.p.value=T,B=2000) # n<-20 xsim<-rbinom(20000,n,.5) xsq<-numeric() for(i in 1:20000) {xsq[i]<-2*((xsim[i]-(n/2))^2/(n/2))} sim.p.value<-length(xsq[xsq>=.8])/20000 # 0.8=2*(12-10)/10 sim.p.value # Diferenca do MC para o Qui-quadrado gera<-rchisq(200000,1) hist(xsq,breaks=200,freq=F,xlim=range(c(0,10))) x<-seq(0.01,10,.001) lines(x,dchisq(x,1),lty=1,col=2) # # Teste via bootstrap: amostrando de Bernoulli com reposição # # amostra observada: 12 sucessos em 20 trials b <- 20000 # número de amostras bootstrap n <- 20 chi.boot<-NULL x<- c(rep(1,12),rep(0,8)) # for (j in 1:b) { # começo do for do bootstrap x.boot<-NULL i. <- sample(n,replace=T) for (i in 1:n){ x.boot[i] <- x[i.[i]] } chi.boot[j] <- 2*((sum(x.boot)-2-10)^2)/10 # -2 corrige pra gerar sob H_0 } # fim do for do bootstrap # hist(chi.boot) boot.p.value<-length(chi.boot[chi.boot>=.8])/b # boot.p.value # # fim do for do bootstrap # # Intervalo de Confiança bootstrap # p.boot<-NULL for (j in 1:b) { # começo do for do bootstrap x.boot<-NULL i. <- sample(n,replace=T) for (i in 1:n){ x.boot[i] <- x[i.[i]] } p.boot[j] <- sum(x.boot)/n # } p.boot<-sort(p.boot) inf_p<-p.boot[b/20] sup_p<-p.boot[b-(b/20)] inf_p sup_p # # 1c. Explorando a Função de Verossimilhança # Amostra Binomial n=20 (x=12) # n<-20 theta<-seq(0.01,1,0.01) vero<-(theta^12)*((1-theta)^8) plot(theta,vero,type="l") # # Log da Razão de Verossimilhanças RV<-vero/((0.6^12)*(0.4^8)) logRV<-log(RV) plot(theta,RV,type="l") plot(theta[logRV>-3],logRV[logRV>-3],type="l") # # Verossimilhança Aproximada (Normal) a<- -(0.5*log(2*pi)) lognor<- -(0.5*log(20*theta*(1-theta)))-((12-(20*theta))^2/ (2*20*theta*(1-theta)))+a logRVnor<- (lognor - (-(0.5*log(20*0.6*(1-0.6)))-((12-(20*0.6))^2/ (2*20*0.6*(1-0.6)))+a)) lines(theta[logRV>-3],logRVnor[logRV>-3],lty=2) # # Intervalos de Confiança (90% do TRV) abline(h=-1.353) # Verossimilhanca Verdadeira LI<-(12*log((20*0.4169)/12))+ (8*log((20*(1-0.4169))/8)) LI LS<-(12*log((20*0.7658)/12))+ (8*log((20*(1-0.7658))/8)) LS # Verossimilhança Aproximada (normal) LI<- 0.5*(-log(20*0.4189*(1-0.4189))-(((12-(20*0.4189))^2)/(20*0.4189*(1-0.4189)))+log(12*0.4)) LI LS<- 0.5*(-log(20*0.7643*(1-0.7643))-(((12-(20*0.7643))^2)/(20*0.7643*(1-0.7643)))+log(12*0.4)) LS # # Teste da RV theta<-0.5 vero<-(theta^12)*((1-theta)^8) RV<-vero/((0.6^12)*(0.4^8)) logRV<--2*log(RV) 1-pchisq(logRV,1) # # # 1.d Inferência Bayesiana - Priori U(0,1) # x<-rbeta(10000,13,9) # beta (13,9) a posteriori hist(x) y<-sort(x) y[250] # 0.3822948 y[9750] # 0.7829416 # I.C. percentilico z<-(13*9)/(22*22*23) # variância a posteriori dif<-1.96*sqrt(z) (13/22)+dif # 0.7918472 (13/22)-dif # 0.3899709 # I.C. assintotico. # #================================================= # 2. Testes para Tabela 2 por 2 # 2a. Estatística Qui-Quadrado #================================================== # Dados Fernando Milton (ortopedia) dados<-matrix(c(11,14,9,166),nc=2) dados Qp<-chisq.test(dados,correct=F) # teste usual Qp Qpc<-chisq.test(dados,correct=T) # teste corrigido Qpc Qpcmc<-chisq.test(dados,correct=F,simulate.p.value = T, B = 50000) # MC teste corrigido Qpcmc # # 2b. Estatística da RV #================================== # grupo <- factor(rep(c("caso", "controle"), c(2,2))) hist <- factor(rep(c("sim", "nao"),2)) freq <- c(14,4,233,281) cancer.data <- data.frame(grupo, hist, freq) cancer.data glm.out1 <- glm(freq~grupo+hist, family=poisson, data=cancer.data) glm.out1 1-pchisq(7.656,1) # # Fernando Milton # injecao <- factor(rep(c("sim", "nao"), c(2,2))) flexao <- factor(rep(c("sim", "nao"),2)) freq <- c(11,14,9,166) fernando.data <- data.frame(injecao, flexao, freq) fernando.data # Outra forma de apresentar os dados mas nao para os modelos (ov<-xtabs(freq~injecao+flexao)) # Usando Modelo de Poisson glm.out2 <- glm(freq~injecao+flexao, family=poisson, data=fernando.data) glm.out2 1-pchisq(24.79,1) # # 2c. Teste Exato de Fisher #============================== dados1<-matrix(c(11,14,9,166),nc=2) fisher.test(dados1) # # Experimento do Chá dados2<-matrix(c(3,1,1,3),nc=2) fisher.test(dados2) # Outro Exemplo dados3<-table(x1,x2) # cruzando as variaveis dados3<-matrix(c(6,3,2,5), nc=2) fisher.test(dados3) # # Fernando Milton fisher.test(dados) # #3. Medidas de Associação #========================================= # # 3a. Odds Ratio = OR e IC95%(OR) #============================== # # dados<-matrix(c(189,104,10845,10933),nc=2) # aspirina dados<-matrix(c(11,14,9,166),nc=2) # Fernando Milton dados OR<-(dados[1,1]*dados[2,2])/(dados[1,2]*dados[2,1]) OR vf<-(1/dados[1,1])+(1/dados[1,2])+(1/dados[2,1]+(1/dados[2,2])) vf dpf<-sqrt(vf) dpf z<-qnorm(0.975) li<-exp(log(OR)-z*dpf) li ls<-exp(log(OR)+z*dpf) ls # # 3b. Risco Relativo = RR e IC95%(RR) #================================== # # dados<-matrix(c(189,104,10845,10933),nc=2)# aspirina dados<-matrix(c(11,14,9,166),nc=2) # Fernando Milton dados p11<-(dados[1,1]/(sum(dados[1,]))) p21<-(dados[2,1]/(sum(dados[2,]))) RR<-p11/p21 RR vf1<-((1-p11)/(sum(dados[1,])*p11)) + ((1-p21)/(sum(dados[2,])*p21)) dpf1<-sqrt(vf1) z<-qnorm(0.975) li<-exp(log(RR)-z*dpf1) li ls<-exp(log(RR)+z*dpf1) ls # # 4. Dados Pareados: # 4a. Teste de Mcnemar #============================================= # dados<-matrix(c(380,150,70,400),nc=2) # dados<-matrix(c(20,10,5,10),nc=2) dados mcnemar.test(dados,correct=F) # # 4b. Concordância: Estatística Kappa #================================================= # # Obs: baixar e instalar: vcd_0.1-3.2.zip (http://www.r-project.org) # require(vcd) x<-c(71,13,6,59) x<-matrix(x,2,2) x Kappa(x) # Exemplo Livro Suely x<-c(38,5,0,1,33,11,3,0,10,14,5,6,3,7,3,10) x<-matrix(x,4,4) Kappa(x) # # 5. Tabelas r por s #============================================================ # dados<-matrix(c(14,29,117,84,37,20),nc=3) dados Qp<-chisq.test(dados,correct=F) # teste usual Qp Qs<-chisq.test(dados,simulate.p.value = T, B = 50000) Qs # genero <- factor(rep(c("menino", "menina"), c(3,3))) diferenca <- factor(rep(c("negativa", "+", "++"),2)) freq <- c(14,117,37,29,84,20) fernando.data <- data.frame(genero, diferenca, freq) fernando.data glm.out1 <- glm(freq~genero+diferenca, family=poisson, data=fernando.data) glm.out1 1-pchisq(7.656,1) # # Residuos Qp$residuals # # 5. Estatística Escore Médio = QS e p-valor (tabela 2 x 3) # Exemplo Suely (Cap. 4) #========================================================= # dados<-matrix(c(13,29,7,7,21,7),nc=3) dados escore<-c(1,2,3) fb1<-(sum(dados[1,]*escore))/sum(dados[1,]) fb2<-(sum(dados[2,]*escore))/sum(dados[2,]) esp<-(c(sum(dados[,1]),sum(dados[,2]),sum(dados[,3])))/sum(dados) mua<-sum(escore*esp) va<-sum((escore-mua)^2*esp) vbf1<-((sum(dados) - sum(dados[1,]))/(sum(dados[1,])*(sum(dados)-1)))*va QS = ((fb1-mua)^2)/vbf1 QS gl<-nrow(dados)-1 p<-1-pchisq(QS,gl) p # # 6. Mantel Haenszel = QMH em tabelas 2 x 2, ORMH e IC(ORMH) #=========================================================== # # Exemplo Câncer de Ovário # tab<-array(c(54,6,62,10,74,3,49,14),dim=c(2,2,2)) tab mantelhaen.test(tab, correct=F) # # 7. Estatística QSMH e valor p em tabelas 2x3 # Exemplo Suely (Cap. 5, p.76, Tabela 5.3) #================================================ # dados<-matrix(c(6,19,7,10,5,7,2,0,16,6,5,1),nc=3) dados escore<-c(0,1,2) fb11<-(sum(dados[1,]*escore))/sum(dados[1,]) fb21<-(sum(dados[3,]*escore))/sum(dados[3,]) c(fb11,fb21) fm1<-sum(c(sum(dados[1,]),sum(dados[3,]))*c(fb11,fb21)) esp1<-(c(sum(dados[1:2,1]),sum(dados[1:2,2]),sum(dados[1:2,3])))/sum(dados[1:2,]) mu1<-sum(escore*esp1) esp2<-(c(sum(dados[3:4,1]),sum(dados[3:4,2]),sum(dados[3:4,3])))/sum(dados[3:4,]) mu2<-sum(escore*esp2) mu<-sum(c(sum(dados[1,]),sum(dados[3,]))*c(mu1,mu2)) v1<- sum(((escore-mu1)^2)*esp1) v2<- sum(((escore-mu2)^2)*esp2) vfma<-(sum(dados[1,])*sum(dados[2,])*v1)/(sum(dados[1:2,])-1) vfmb<-(sum(dados[3,])*sum(dados[4,])*v2)/(sum(dados[3:4,])-1) vfm<- sum(c(vfma,vfmb)) QSMH<-((fm1-mu)^2)/vfm p<-1-pchisq(QSMH,1) round(c(QSMH,p),digits=5) # # 8. Modelos Log-lineares - três Dimensões #========================================== # 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 # glm.out1<-glm(formula=freq ~alcool + cigarro + maconha, family=poisson, data=alunos) glm.out1 glm.out1<-glm(formula=freq ~alcool*cigarro*maconha, family=poisson, data=alunos) anova(glm.out1, test="chisq") # # 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. #============================================================================== # # Exemplo da Lista 3: Exerc. 6 # 1- meio de comunicação (tv, jornal, rádio) # 2- idade (< 40 e > 40) # 3- Escolaridade (1o, 2o e sup.) # tipo<- factor(rep(rep(c("tv","tv","jornal","jornal","radio","radio")),3)) idade<-factor(rep(c("mais","menos"),9)) escola<-factor(rep(c("1","2","sup"),c(6,6,6))) freq<-c(52,14,27,8,22,18,34,13,23,11,6,15,15,8,35,12,7,3) folha<-data.frame(tipo,idade,escola,freq) folha glm.out0<-glm(formula=freq ~1, family=poisson, data=folha) glm.out1<-glm(formula=freq ~tipo + idade + escola, family=poisson, data=folha) anova(glm.out0,glm.out1,test="Chisq") glm.out11<-glm(formula=freq ~tipo * escola+ idade, family=poisson, data=folha) glm.out2<-glm(formula=freq ~tipo + escola + idade + tipo:idade + tipo:escola, family=poisson, data=folha) anova(glm.out11,glm.out2,test="Chisq") glm.out21<-glm(formula=freq ~tipo + escola + idade + tipo:escola + idade:escola, family=poisson, data=folha) glm.out21<-glm(formula=freq ~tipo + escola + idade + tipo:idade + tipo:escola + idade:escola, family=poisson, data=folha) glm.out3<-glm(formula=freq ~tipo * escola * idade, family=poisson, data=folha) #=========================================================================================