#=========================================
#     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),col=1)
 lines(x,dchisq(x,3),col=2)
 lines(x,dchisq(x,4),col=3)
 lines(x,dchisq(x,7),col=4)
 lines(x,dchisq(x,10),col=5)
 lines(x,dchisq(x,15),col=6)
 legend(15,0.2, legend = c("gl=1","gl=3","gl=4","gl=7","gl=10","gl=15"), 
 lty=c(1,1,1,1,1,1),col = c(1,2,3,4,5,6),  bty="n", cex = 1.0)
#
 pchisq(3.84,1)
 1-pchisq(3.84,1)
 qchisq(0.95,1)
 gera<-sort(rchisq(50000,3))
 hist(gera,breaks=80,freq=F,xlab="x")
 lines(x<-seq(.1,15,.1),dchisq(x,3),col=2)
#====================================================
#
# 1b. Explorando a Distribuição Binomial
#     Testes para uma única amostra: exato
#     qui-quadrado, MC e Bootstrap
#=====================================================
#
 help(Binomial)
 x<-0:20
 prob<-round(dbinom(x,20,0.5),4)
 plot(x,prob,type="h")
#
#     Exemplo Ecologia: n=20, x=12 # número de fêmeas
#
 barplot(c(12,8),names.arg =c("fêmea","macho"),ylab="frequência")
 barplot(c(12,8)*5,names.arg =c("fêmea","macho"),ylab="porcentagem",
 col=c("pink","blue"))
#
#   Teste Exato - Uma Amostra - H_0: \pi = 0.5
#
 binom.test(12,20)  # testa H_0:p=0,5 para n=20 e x=12
# valor-p = 1 - P(8 < X < 12) sob H_0: 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.15,0.90,0.01)
 x2<-(((12-(20*x))^2)/(20*x))+((((20*x)-12)^2)/(20*(1-x)))
 dados1<-cbind(x,x2)
# dados1 
 plot(x,x*0,pch="",ylim=range(c(0,10)),xlim=range(c(0.2,0.9)),
         xlab="pi",ylab="X2",bty="n",
       main="Int. 95% de Confiança baseado em X^2")
 lines(x,x2)
 abline(h =3.84, untf = FALSE,col=2)
 arrows(0.382,1.5,0.382,-0.3,length=0.12,col="blue")
 arrows(0.782,1.5,0.782,-0.3,length=0.12,col="blue")
#
#  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)
#
#  Intervalo de Confiança baseado no MC
#    gerar sob os dados \pi = 0.6
#
 n<-20
 xsim<-rbinom(20000,n,12/20)
 xic<-sort(xsim)
 inf<-xic[1000]
 sup<-xic[19000]
 inf/n
 sup/n
#  IC --> (0,4;0,8) ----extremamente discreto
#
# 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/40]
 sup_p<-p.boot[b-(b/40)]
 inf_p
 sup_p
#
#
# 1.c Inferência Bayesiana - Priori U(0,1)
#
  x<-rbeta(10000,13,9) # beta (13,9) a posteriori
  hist(x,xlab="pi",ylab="frequência")
  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.
#
# Usando o pacote TeachingDemos
 library("TeachingDemos")
# Encontrando o IC HPD
 hpd(qbeta,shape1=13,shape2=9)
#
#
#  1d. Explorando a Função de Verossimilhança
#        Amostra Binomial n=20 (x=12) 
#
 n<-20
 phi<-seq(0.01,1,0.01)
 vero<-(phi^12)*((1-phi)^8)
 plot(phi,vero,type="l")
#
# Log da Razão de Verossimilhanças
 RV<-vero/((0.6^12)*(0.4^8))
 logRV<-log(RV)
 plot(phi,RV,type="l")
 plot(phi[logRV>-3],logRV[logRV>-3],type="l")
 plot(phi[logRV>-3],-2*logRV[logRV>-3],type="l")
 abline(h =3.84, untf = FALSE,col=2)
#
# Verossimilhança Aproximada (Normal)
 a<- -(0.5*log(2*phi))
 lognor<- -(0.5*log(20*phi*(1-phi)))-((12-(20*phi))^2/
 (2*20*phi*(1-phi)))+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))
  plot(phi[logRV>-3],logRV[logRV>-3],type="l",ylab="logRV",xlab="phi")
 lines(phi[logRV>-3],logRVnor[logRV>-3],lty=2,col=2)
#
# Intervalos de Confiança (90% do TRV)
  h=-(1.96^2)/2
  h
  abline(h=h)
# Verossimilhanca Verdadeira
  LI<-(12*log((20*0.383)/12))+ (8*log((20*(1-0.383))/8))
  LI
# LI=0.383
  LS<-(12*log((20*0.793)/12))+ (8*log((20*(1-0.793))/8))
  LS
# LS=0.793
# Verossimilhança Aproximada (normal)
  theta=0.386
  LI<- 0.5*(-log(20*theta*(1-theta))-(((12-(20*theta))^2)/
  (20*theta*(1-theta)))+log(12*0.4))
  LI
 #
# Teste da RV
  phi<-0.5
  vero<-(phi^12)*((1-phi)^8)
  RV<-vero/((0.6^12)*(0.4^8))
  logRV<--2*log(RV)
  1-pchisq(logRV,1)
#
#================================================= 
# 2. Testes para Tabela 2 por 2
#=====================================================
#
#     2a. Descrevendo a Tabela  
#==================================================
#     Exemplo AZT
#  
  azt<-matrix(c(1,16,144,121),nc=2)
  dimnames(azt) <- list(Grupo = c("AZT", "Placebo"), 
                  Óbito = c("Sim", "Não")) 
  azt
  margin.table(azt)
  margin.table(azt,2)
  margin.table(azt,1)
  prop.table(azt,1)
  barplot(prop.table(azt,1),beside=T, ylim= c(0,1.3),legend=T)
  require(vcd) 
  tile(prop.table(azt,1))
  spineplot(prop.table(azt,1))
#    
#
#     2b. Estatística Qui-Quadrado
#==================================================
#
#   Exemplo AZT
#
   Qp<-chisq.test(azt,correct=F) # teste usual
   Qp
#
#   Dados Fernando Milton (ortopedia)
 dados<-matrix(c(11,14,9,166),nc=2)
 dimnames(dados) <- list(Injeção = c("Inj. Sim", "Inj. Não"), 
                  Flexão = c("Dif. Flex: Sim", "Dif. Flex: Não")) 
 dados
 prop.table(dados,1)
 barplot(prop.table(dados,1),beside=T, ylim= c(0,1.2),legend=T)
#  
 Qp<-chisq.test(dados,correct=F) # teste usual
 Qp
# Qpc<-chisq.test(dados,correct=T) # teste corrigido
# Qpc
#
#    2c. Estatística da RV
#=====================================================
#
#  Exemplo AZT
  grupo <- factor(rep(c("AZT", "Placebo"), c(2,2)))
  hist  <- factor(rep(c("Vivo", "Morto"),2))
  freq <- c(144,1,121,16)
  azt.data <- data.frame(grupo, hist, freq)
  azt.data
  glm.out1 <- glm(freq~grupo+hist, family=poisson)
  glm.out1
  1-pchisq(glm.out1$deviance,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(glm.out2$deviance,1)
#
#  2d. Teste Computacionalmente Intensivo
#       MC sob H_0
#========================================================================
#
  Qpcmc<-chisq.test(dados,correct=F,simulate.p.value = T, 
  B = 2000) # MC teste corrigido
  Qpcmc
#
#
#  2e. Teste Exato de Fisher
#=================================================================
# Exemplo Fernando Milton 
 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)
#
#   3. Medidas de Associação
#=========================================
#
# 3a. Razão de Chances = RC e IC95%(RC)
#==========================================
#
# 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
#
#  Obs.: IC para RC baseado no modelo logístico
#    Duas binomias---> modelo saturado com 2 gl
#
  placebo<-c(1,0)
  attack<-c(189,104)
  n<-c(189+10845,104+10933)
  fit<-glm(attack/n~placebo, weight=n, family=binomial)
  summary(fit)
# obter IC baseado na verossimilhança perfilada
#     este IC deve ser preferido ao Wald
  exp(confint(fit)) 
#
#   Usando o Pacote Epitools
#
  require(epitools)
  oddsratio(dados)
  azt<-matrix(c(144,121,1,16),nc=2)
  azt
  oddsratio (azt, method = "wald") # comum
  oddsratio (azt)
  oddsratio (azt, method = "midp")
  oddsratio (azt, method = "fisher")
  oddsratio (azt, method = "wald")
 oddsratio (azt, method = "small")
#
# 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
#=============================================
#    Efeito de Propaganda Eleitoral
 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
 K<-Kappa(x)
 confint(K)
 summary(K)
 print(K, CI = TRUE)
#  
# 5. Tabelas r por s
#============================================================
#
 dados<-matrix(c(14,29,117,84,37,20),nc=3)
 dados
 dimnames(dados) <- list(Sexo = c("Menino", "Menina"), 
                  Dif.Idade = c("-9 a -1", "0 a 5", "5 a 15"))
#Gráfico de barras colorido para a proporção de menino
 p1=dados[1,1]/sum(dados[,1])
 p1
 p2=dados[1,2]/sum(dados[,2])
 p2
 p3=dados[1,3]/sum(dados[,3])
 p3
 barplot(c(p1,p2,p3), names.arg=c("-9 a -1", "0 a 5", "5 a 15"), xlab="Difernça de idade", 
  ylab="proporção menino", col=c("pink", "green", "blue"))
#
 prop.table(dados,1)
 barplot(prop.table(dados,1),beside=T, ylim= c(0,1.2),legend=T)
#
 Qp<-chisq.test(dados,correct=F) # teste usual
 Qp
 Qp$residuals
 Qs<-chisq.test(dados,simulate.p.value = T, 
 B = 50000) 
 Qs
#
#   Teste da Razão de Verossimilhanças
  genero <- factor(rep(c("menino", "menina"), c(3,3)))
  diferenca  <- factor(rep(c("negativa", "+", "++"),2))
  freq <- c(14,117,37,29,84,20)
  nature.data <- data.frame(genero, diferenca, freq)
  nature.data
  glm.out1 <- glm(freq~genero+diferenca, family=poisson, data=nature.data)
  glm.out1
  1-pchisq(glm.out1$deviance,2)
#
# Residuos
  names(Qp)
  Qp$residuals
  Qp$stdres
#
 require(epitools)
 dados1<-matrix(c(117,84,14,29),nc=2)
 dados1
 oddsratio(dados1)
 dados2<-matrix(c(37,20,14,29),nc=2)
 dados2
 oddsratio(dados2)
#
# 5. Estatística Escore Médio = QS e p-valor (tabela 2 x 3)
#    Exemplos: Nature e  Suely (Cap. 4)
#=========================================================
#   Exemplo Artificial dos slides
#
 dados<-matrix(c(10,90,20,80,40,60),nc=3)
 dados
 prop.table(dados,1)
 barplot(prop.table(dados,1),beside=T, ylim= c(0,1.2),legend=T)
 Qp<-chisq.test(dados,correct=F) # teste usual
 Qp
#
 dados1<-matrix(c(10,90,40,60,20,80),nc=3)
 dados1
 prop.table(dados1,1)
 barplot(prop.table(dados1,1),beside=T, ylim= c(0,1.2),legend=T)
 Qp<-chisq.test(dados1,correct=F) # teste usual
 Qp
#
 dados
# dados<-dados1
 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
#
#  Exemplo: Nature
# Já que a diferença de idades está ordenada, podemos
# realizar um teste de tendência. Ou seja, será que 
# probabilidade de ocorrer menina aumenta com a diferença
# de idade entre Pai e Mãe (idade do pai - idade
# da mãe)? Os escores são os valores médio das classes (-5; 2,5;10) 
#
 dados<-matrix(c(14,29,117,84,37,20),nc=3)
 dados
 escore<-c(-5,2.5,10)
 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 Vendas Produto X
#
#  Cidade A
 dados1<-matrix(c(60,320,140,480),nc=2)
 dados1
 Qp<-chisq.test(dados1,correct=F) # teste usual
 Qp
 OR<-(dados1[1,1]*dados1[2,2])/(dados1[1,2]*
 dados1[2,1])
 OR
#
#   Cidade B
 dados2<-matrix(c(640,180,160,20),nc=2)
 dados2
 Qp<-chisq.test(dados2,correct=F) # teste usual
 Qp
 OR<-(dados2[1,1]*dados2[2,2])/(dados2[1,2]*
 dados2[2,1])
 OR
#
#  Combinado (Cidade A + Cidade B) 
 dados3<-matrix(c(700,500,300,500),nc=2)
 dados3
 Qp<-chisq.test(dados3,correct=F) # teste usual
 Qp
 OR<-(dados3[1,1]*dados3[2,2])/(dados3[1,2]*
 dados3[2,1])
 OR
#
#    Teste de Mantel-Haenszel
 tab<-array(c(60,320,140,480,640,180,160,20),dim=c(2,2,2))
 tab
 mantelhaen.test(tab, correct=T)
#
#  Exemplo Clínico - Comparação de Medicamentos - Giolo, p. 71
#
 tab<-array(c(29,14,16,31,37,24,8,21),dim=c(2,2,2))
 tab
 mantelhaen.test(tab, correct=F)
#
# 7. Estatística QSMH e valor p em tabelas 2x3
#     Exemplo Giolo (Cap. 5, p.92, Tabela 5.4)
#     Teste de Tendência e MH
#================================================
#
 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)
#
#=========================================================================================