#=========================================
#     ANÁLISE DE DADOS CATEGÓRICOS
#     Profa.: SUELY RUIZ GIOLO
#     Prof.: Enrico A. Colosimo
#     Terceira Parte: Modelos: Logístico
#                     
#================================================================
#  1. Regressão Logística Dicotômica
#================================================================
#
# 1.1 Exemplo - Capítulo 7 - Tabela 7.1
# ===============================================================
#  
# Vamos utilizar os dados brutos
#    n=100 Coorte de 100 indivíduos
#     Y (resposta): 0/1 (incidência de DC)
#     X: idade
#==================================================================
   idade<-round(c(runif(10,20,29),runif(15,30,34),runif(12,35,39),
   runif(15,40,44),runif(13,45,49),runif(8,50,54),runif(17,55,59),
   runif(10,60,69)))
   idade
   y<-c(rep(1,1),rep(0,9),rep(1,2),rep(0,13),rep(1,3),rep(0,9),
   rep(1,5),rep(0,10),rep(1,6),rep(0,7),rep(1,5),rep(0,3),
   rep(1,13),rep(0,4),rep(1,8),rep(0,2))
   y
   sum(y)
#  Gráfico de Dispersão
   plot(idade,y,xlab="idade",ylab="ocorrência de DC")
#
#  Ajuste de um Modelo Linear ---> Inadequado
#
   outlm<-lm(y~idade)
   summary(outlm)
   par(mfrow=c(2,2))
   plot(outlm)
   summary(lm(abs(residuals(outlm))~fitted(outlm))) # verificando homocedasticidade
   shapiro.test(outlm$resid) # verificando normalidade
#
   par(mfrow=c(1,1))
   plot(idade,y,abline(outlm$coefficients[1], outlm$coefficients[2]),
   xlab="idade", ylab="ocorrência de DC")
#
#  Vamos agrupar as idades por faixa etária
#
  resim<-c(1,2,3,5,6,5,13,8)
  resnao<-c(9,13,9,10,7,3,4,2)
  idade<-c(25,32,38,43,47,53,57,65)
  dados<-cbind(resim, resnao,idade)
  dados
  dados<-as.data.frame(dados)
# *Outra Forma de Entrar com os Dados
#
  x<-c(rep(25,10),rep(32,15),rep(38,12),
  rep(43,15),rep(47,13),rep(53,8),rep(57,17),
  rep(65,10))
  y<-c(rep(1,1),rep(0,9),rep(1,2),rep(0,13),
  rep(1,3),rep(0,9),rep(1,5),rep(0,10),
  rep(1,6),rep(0,7),rep(1,5),rep(0,3),
  rep(1,13),rep(0,4),rep(1,8),rep(0,2))
  dados1<-cbind(y,x)
  dados1<-as.data.frame(dados1)
  dados1
#
# Gráfico dos Dados Agrupado
#
  theta<-resim/(resim+resnao)
  dados2<-cbind(resim, resnao,idade,theta)
  dados2
  plot(idade,theta,ylim=range(0,0.9),xlab="idade",ylab="E(Y|x)",pch=16)
  abline(outlm$coefficients[1], outlm$coefficients[2])
#
# Forma similar mostrando que o modelo logístico é o default
  ajust<-glm(as.matrix(dados[,c(1,2)])~idade,family=binomial(link="logit"),data=dados)
  summary(ajust)
  anova(ajust,test="Chisq")
#
 ajust1<-glm(y~x,family="binomial",data=dados1)
 summary(ajust1)
 anova(ajust1,test="Chisq")
#
# Gráfico com os dados, ajuste logístico (vermelho) e ajuste linear
#
 idade<-c(25,32,38,43,47,53,57,65)
 plot(idade,theta,ylim=range(0,0.9),xlab="idade",ylab="E(Y|x)",pch=16)
 idade<-20:70
 modajust<-(exp(-5.123+0.1058*idade))/(1+ exp(-5.123+0.1058*idade))
 modajust
 lines(idade,modajust,col=2)
 outlm<-lm(y~x)
 abline(outlm$coefficients[1], outlm$coefficients[2])
#
 ajust$fitted.values
 ajust$y
 ajust$residuals
 dev<-residuals(ajust,type='deviance')
 dev
 QL<-sum(dev^2)
 QL
#   Valor-p para o teste de adequação do Desvio 
 p1<-1-pchisq(QL,ajust$df.residual)
 p1
# 
 rpears<-residuals(ajust,type='pearson')
 rpears
 QP<-sum(rpears^2)
 QP
#   Valor-p para o teste de adequação do Qui-quadrado 
 p2<-1-pchisq(QP,ajust$df.residual)
 p2
 theta<-resim/(resim+resnao)
# 
 require(hnp)
 hnp(ajust1, sim=19, conf=1)
# Intervalo de 95% de confiança para RC
 exp( ajust$coefficients[2])
 LS<-exp(ajust$coefficients[2] + 1.96* 0.02337) 
 LI<-exp(ajust$coefficients[2] - 1.96* 0.02337)
 LI
 LS 
#
# Apresentação dos Resultados e Diagnóstico
#
 require(LogisticDx)
 dx(ajust)
 OR(ajust)
 plot(ajust)
#
#================================================================================
# 1.2 Exemplo 2
#
#      Somente Sexo  
#       Y: doença coronariana
#       n= 79
#=====================================================================================
 resim<-c(12,30)
 resnao<-c(21,15)
 sexo<-c(1,0)# 1 - feminino e 0 - masculino
 dados<-cbind(resim, resnao,sexo)
 dados
 theta<-resim/(resim+resnao)
 theta
#
#   Razão de Chances
#
  RC<- (12/21)/(30/15)
  RC
# [1] 0.2857143
#
  RC1<-1/RC
  RC1
#[1] 3.5
#
 dados<-data.frame(dados)
 ajust<-glm(as.matrix(dados[,c(1,2)])~sexo,family=binomial(link="logit"),data=dados)
 ajust
 summary(ajust)
 rc<- exp(ajust$coefficients[2])
 rc
#
# 1.3 Exemplo 1 - Tabela 7.7 - Pg 135 
# ==================================================
#  Resposta: doença coronariana
#   X1: sexo (0/1)
#   X2: ecg (0/1) alterado: sim/não
#
#===========================================================================
 resim<-c(4,8,9,21)
 resnao<-c(11,10,9,6)
 sexo<-c(1,1,0,0)
 ecg<-c(0,1,0,1)
 dados<-cbind(resim, resnao,sexo,ecg)
 dados
 theta<-resim/(resim+resnao)
 theta
 dados<-as.data.frame(dados)
 ajust<-glm(as.matrix(dados[,c(1,2)])~sexo+ecg,family=binomial,data=dados)
 ajust
 summary(ajust)
 anova(ajust,test="Chisq")
 names(ajust)
 ajust$residuals
 ajust$fitted.values
 theta
 ajust$y
 dev<-residuals(ajust,type='deviance')
 dev
 QL<-sum(dev^2)
 QL
 p1<-1-pchisq(QL,1)
 p1
 rpears<-residuals(ajust,type='pearson')
 rpears
 QP<-sum(rpears^2)
 QP
 p2<-1-pchisq(QP,2)
 p2
 logito<-log(ajust$fitted.values/(1-ajust$fitted.values))
 logito
 odds<-exp(ajust$coefficients)
 odds
 1/odds
 lower<-exp(ajust$coefficient - (1.96* 0.498))
 upper<-exp(ajust$coefficient + (1.96* 0.498))
 lower
 upper
 1/lower
 1/upper
#  
# Encontrando o Deviance (null e residual)
#
 resim<-c(1,0)
 resnao<-c(49,50)
 x<-c(1,-1)
 dados<-cbind(resim, resnao,x)
 dados
 dados<-as.data.frame(dados)
 attach(dados)
 ajust<-glm(as.matrix(dados[,c(1,2)])~x,family=binomial, data=dados)
 ajust
#
# Null Deviance
  2*((49*log(49))-(50*log(50))+(100*log(100))-(99*log(99)))
# [1] 1.396396
  temp<-exp(ajust$coefficients[1]+ajust$coefficients[2])
  pi1<-temp/(temp+1)
  temp1<-exp(ajust$coefficients[1]-ajust$coefficients[2])
  pi2<-temp1/(temp1+1)
  lcur<-log(pi1)+(49*log(1-pi1))+(50*log(1-pi2))
#
# Residual Deviance
   -2*(lcur-(49*log(49))+(50*log(50)))
#
#
# 1.4 Exemplo 2 - Capítulo 7
# ===========================
#
 resim<-c(78,101,68,40,54,34)
 resnao<-c(28,11,46,5,5,6)
 diag<-c(1,1,1,0,0,0)
 tratA<-c(1,0,0,1,0,0)
 tratB<-c(0,1,0,0,1,0)
 int1<-diag*tratA
 int2<-diag*tratB
 dados<-cbind(resim, resnao,diag,tratA,tratB,int1,int2)
 dados
 dados<-as.data.frame(dados)
 attach(dados)
 ajust1<-glm(as.matrix(dados[,c(1,2)])~diag+tratA+tratB+int1+int2,family=binomial
                      (link="logit"),data=dados)
 ajust1
 summary(ajust1)
 anova(ajust1)
 ajust<-glm(as.matrix(dados[,c(1,2)])~diag+tratA+tratB,family=binomial(link="logit"),data=dados)
 ajust
 ajust$fitted.values
 ajust$y
 dev<-residuals(ajust,type='deviance')
 dev
 QL<-sum(dev^2)
 QL
 p1<-1-pchisq(QL,2)
 p1
 rpears<-residuals(ajust,type='pearson')
 rpears
 QP<-sum(rpears^2)
 QP
 p2<-1-pchisq(QP,2)
 p2
 logito<-log(ajust$fitted.values/(1-ajust$fitted.values))
 logito
 odds<-ajust$fitted.values/(1-ajust$fitted.values)
 odds
#
# 1.5 Exemplo 3 - Capítulo7
# ===========================
#
 dc<-c(0,0,0,1,0,1,0,0,0,0,0,0,1,0,1,1,0,0,0,0,1,1,0,0,0,0,1,1,0,0,1,1,0,0,1,1,1,0,1,
       1,0,1,0,0,0,1,1,0,1,1,0,1,1,0,0,1,1,0,0,0,1,1,1,1,1,1,1,1,1,0,1,1,1,1,0,1,1,1)
 sexo<-c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
         1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
 ecg<-c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,
        1,1,1,1,1,1,1,1,1,1,2,2,2,2,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2)
 idade<-c(28,34,38,41,44,45,46,47,50,51,51,53,55,59,60,32,33,35,39,40,42,44,45,46,48,50,
         52,52,54,55,59,59,32,37,38,38,42,43,43,44,46,48,49,49,52,53,54,55,57,46,48,57,
         60,30,34,36,38,39,42,45,45,45,46,48,57,57,59,60,63,35,37,43,47,48,49,58,59,60)
#
 ajust1<-glm(dc~sexo+ecg+idade+sexo*ecg+sexo*idade+ecg*idade+sexo*ecg*idade, family=binomial(link="logit"))
 ajust1
 summary(ajust1)
 anova(ajust1,test="Chisq")
 ajust2<-glm(dc~sexo+ecg+idade,family=binomial(link="logit"))
 ajust2
 summary(ajust2)
 anova(ajust2, test="Chisq")
 cbind(dc,sexo,ecg,idade,ajust2$fitted.values)
 dev<-residuals(ajust2,type='deviance')
 dev
 plot(dev)
 rpears<-residuals(ajust2,type='pearson')
 rpears
 plot(rpears)
#
#=====================================================================
#   2 - Tópicos Especiais - Regressão Logística Binária
#============================================================
#   
#  2.1 Teste de Hosmer e Lemeshow
#
#
 hosmerlem = function(y, yhat, g=10) {
  cutyhat = cut(yhat,
  breaks = quantile(yhat, probs=seq(0,
  1, 1/g)), include.lowest=TRUE)
  obs = xtabs(cbind(1 - y, y) ~ cutyhat)
  expect = xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
  chisq = sum((obs - expect)^2/expect)
  P = 1 - pchisq(chisq, g - 2)
  return(list(chisq=chisq,p.value=P))
  }
#
#==========================================================================
#
#  2.2  Construção do Nomograma
#       Útil na interpretação dos resultados
#============================================================
# Referência: Harrell (2000)
# Necessário: pacote Design
# Exemplo de interpretação do Nomograma em
# http://www.prostate-cancer.org/education/riskases/Tisman_UsingNomograms1.html
#============================================================
#
 require(Design)
 ajust<-lrm(dc~factor(sexo)+factor(ecg)+idade)
 ddist <- datadist(idade, sexo, ecg)
 options(datadist='ddist')
 nomogram(ajust,fun=plogis,funlabel="probabilidade", 
 fun.at=c(.01,.05,.1,.25,.5,.75,.90,.95,.99), xfrac=.45)
#
#===============================================================
# 2.3  Q-QPLOT com envelope simulado (código: envel.bino)
#
#==============================================================
#
 fit.model<-ajust2   
 par(mfrow=c(1,1))
 X <- model.matrix(fit.model)
 n <- nrow(X)
 p <- ncol(X)
 w <- fit.model$weights
 W <- diag(w)
 H <- solve(t(X)%*%W%*%X)
 H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
 h <- diag(H)
 td <- resid(fit.model,type="deviance")/sqrt(1-h)
 e <- matrix(0,n,100)
#
 for(i in 1:100){
  dif <- runif(n) - fitted(fit.model)
  dif[dif >= 0 ] <- 0
  dif[dif<0] <- 1
  nresp <- dif
  fit <- glm(nresp ~ X, family=binomial)
  w <- fit$weights
  W <- diag(w)
  H <- solve(t(X)%*%W%*%X)
  H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
  h <- diag(H)
  e[,i] <- sort(resid(fit,type="deviance")/sqrt(1-h))}
  e1 <- numeric(n)
  e2 <- numeric(n)
#
 for(i in 1:n){
  eo <- sort(e[i,])
  e1[i] <- eo[5]
  e2[i] <- eo[95]}
  med <- apply(e,1,mean)
  faixa <- range(td,e1,e2)
  par(pty="s")
  qqnorm(td,xlab="Percentis da N(0,1)",
  ylab="Componente da deviance", ylim=faixa, pch=16)
  par(new=T)
#
 qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1)
 par(new=T)
 qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1)
 par(new=T)
 qqnorm(med,axes=F,xlab="", ylab="", type="l",ylim=faixa,lty=2)
#
#===================================================================
# 2.4 - Usando o brglm (verossimilhanca penalizada - Schemper 02)
#        Verossimilhança Monótona
#        Mudando a entrada de dados (Binomial e Bernoulli)
#==================================================================
#  Ajustes sem penalização de duas formas.
  n<-c(21,17,17)
  x<-c(4,2,0)
  groups<-factor(1:3)
#
  dados1<-cbind(n,groups,x,n-x)
  aj<-glm(dados1[,c(3,4)]~groups,binomial)
  summary(aj)
  D1<-aj$null.deviance-aj$deviance
  D1
#
# Mudando a entrada de dados: Binomial --> Bernoulli
# 
  dados2<-data.frame(n,groups,x,n-x)
  dados3<-rbind(cbind(case=1,dados2[rep(1:3,dados2$x),1:3]),cbind(case=0,dados2[rep(1:3,dados2$n-dados2$x),1:3]))
  ajj<-glm(dados3$case~dados3$groups,binomial)
  summary(ajj)
#
  D2<-ajj$null.deviance-ajj$deviance
  D2
#
#D1 e D2 são identicos, logo tanto faz eu ajustar aj ou ajj para meu objetivo, estou trabalhando com aj.
#
# Ajustes penalisados.
#
  require(brglm)
  n<-c(21,17,17)
  x<-c(4,2,0)
  groups<-factor(1:3)
#
  dados1<-cbind(n,groups,x,n-x)
#
  ai<-brglm(dados1[,c(3,4)]~groups,binomial,method="brglm.fit")
  summary(ai)
  D11<-ai$null.deviance-ai$deviance
  D11
#   
#===================================================================
# 2.5 - Usando a Var robusta (Sanduíche)
#        
#==================================================================
#
#I'm not sure what you mean here but I include a worked example. Caveat:
#The data I use are cross-section data with an overly simplified set of
#regressors. So none of this makes sense for the application - but it shows
#how to use the commands.
#
## load AER package which provides the example data
## and automatically loads "lmtest" and "sandwich"
 library("AER")
 data("PSID1976", package = "AER")
#
## fit a simple logit model and obtain marginal Wald tests
## for the coefficients and an overall chi-squared statistic
 m <- glm(participation ~ education, data = PSID1976, family = binomial)
 summary(m)
 anova(m, test = "Chisq")
#
## replicate the same statistics with coeftest() and lrtest()
 coeftest(m)
 lrtest(m)
#
## the likelihood ratio test is asymptotically equivalent
## to the Wald test leading to a similar chi-squared test here
 waldtest(m)
#
## obtain HAC-corrected (Newey-West) versions of the Wald tests
 coeftest(m, vcov = NeweyWest)
 waldtest(m, vcov = NeweyWest)
#
# Instead of NeweyWest other covariance estimators (e.g., vcovHAC, kernHAC,
# etc.) can also be plugged in.             
#
#================================================================================
# 3. Regressão Logística Politômica
#====================================================================
#
# 3.1 Regressão Logística Ordinal
#========================================================================
 set.seed(1)
 dados <- data.frame(grupo=sample(3, 100, rep=T), sexo=sample(2, 100, rep=T))
 require(VGAM)
 require(Design)
 modelo<-rrvglm(grupo ~ sexo, data=dados, multinomial)
 summary(modelo)
#
#  3.2 Exemplo: modelo de odds proporcionais
# ==========================================
#
 require(MASS)
 melhora<-rep(c("ac","alg","nenh"),4)
 sexo<-c(1,1,1,1,1,1,0,0,0,0,0,0)
 trat<-c(1,1,1,0,0,0,1,1,1,0,0,0)
 Freq<-c(16,5,6,6,7,19,5,2,7,1,0,10)
 artrite<-cbind(melhora,sexo,trat)
 artrite<-as.data.frame(artrite)
 attach(artrite)
 options(contrasts = c("contr.treatment", "contr.poly"))
 ajust1 <- polr(melhora ~ sexo + trat + sexo*trat, weights = Freq, data=artrite)
 ajust1
 summary(ajust1)
 ajust2 <- polr(melhora ~ sexo + trat, weights= Freq, data = artrite)
 ajust2
 summary(ajust2)
 ajust2$fitted.values
#
# *****************************************************************************
# * Obs: inverter sinais dos parâmetros dos efeitos e manter os dos interceptos 
# *****************************************************************************
#
#===========================================
# 3.3. Regressão Logística Condicional
#===================================
#
# 15.1 Exemplo: estudo retrospectivo
# ==================================
#
 skin<-read.table("http://www.est.ufpr.br/~suely/CE215/Dados/skin.txt",h=T)
 attach(skin)
 require(survival)
 model1<-clogit(melhora~trat+sexo+idade+grauini+strata(clinica))
 model1
 summary(model1)
 plot(model1$residuals, pch=16)
 model2<-clogit(melhora~trat+grauini+strata(clinica))
 model2
 summary(model2)
 plot(model2$residuals, pch=16,ylab="residuos",xlab="i")
#
#
# 15.2 Exemplo: estudo crossover
# ==============================
#
 cross<-read.table("http://www.est.ufpr.br/~suely/CE215/Dados/cross.txt",h=T)
 attach(cross)
#
# preparando os dados para análise
# ----------------------------------
#
 n<-sum(freq)
 m<-dim(cross)[2]
 k<-dim(cross)[1]
 cross1<-matrix(0,n,m)
 cross2<-as.data.frame(cross1)
 count<-c(0,freq)
 for(j in 1:k){
  for(i in (sum(count[1:j])+1):(sum(count[1:(j+1)]))){
   cross2[i,] <- cross[j,]}
  }
 names(cross2)<-names(cross)
 obs<-1:300
 cross2$obs<-obs
 cross3<-as.data.frame(rbind(cross2,cross2))
 i<-order(cross3$obs)
 cross4<-cross3[i,]
#
#-----------------------------------------------------------
# idade: 1 se adulto e 0 se jovem, F = 1 e U = 2
# sequencias: AB = 1, AP = 2, BA= 3, BP = 4, PA = 5 e PB = 6
# criando variaveis dummies
#-----------------------------------------------------------
#
 periodo<-rep(c(1,0),300)    #1 se periodo1 e 0 se periodo2
 cross4$periodo<-periodo
 drogaA<-c(rep(c(1,0),50),rep(c(0,0),50),rep(c(0,1),50),
           rep(c(0,1),50),rep(c(1,0),50),rep(c(0,0),50))
 cross4$drogaA<-drogaA
 drogaB<-c(rep(c(0,1),50),rep(c(1,0),50),rep(c(0,0),50),
           rep(c(1,0),50),rep(c(0,0),50),rep(c(0,1),50))
 cross4$drogaB<-drogaB
 resA<-c(rep(c(0,1),50),rep(c(0,0),50),rep(c(0,0),50),
         rep(c(0,0),50),rep(c(0,1),50),rep(c(0,0),50))
 cross4$resA<-resA
 resB<-c(rep(c(0,0),50),rep(c(0,1),50),rep(c(0,0),50),
         rep(c(0,1),50),rep(c(0,0),50),rep(c(0,0),50))
 cross4$resB<-resB
 attach(cross4)
 resp<-rep(0,600)
 for(i in 1:600){
  ifelse(p1[i]==1 & p2[i]==1, resp[i]<-1,resp[i])
  ifelse(p1[i]==1 & p2[i]==2 & periodo[i]==1, resp[i]<-1,resp[i])
  ifelse(p1[i]==1 & p2[i]==2 & periodo[i]==0, resp[i]<-0,resp[i])
  ifelse(p1[i]==2 & p2[i]==1 & periodo[i]==1, resp[i]<-0,resp[i])
  ifelse(p1[i]==2 & p2[i]==1 & periodo[i]==0, resp[i]<-1,resp[i])
  ifelse(p1[i]==2 & p2[i]==2, resp[i]<-0, resp[i])
 }
 cross4$resp<-resp
 gpidade<-periodo*idade;
 cross4$gpidade<-gpidade
#
# Usando arquivo cross4 para ajustar modelos
# --------------------------------------------
#
 attach(cross4)
 require(survival)
 model1<-clogit(resp~periodo+drogaA+drogaB+gpidade+resA+resB+strata(obs),data=cross4)
 model1
 summary(model1)
 plot(model1$residuals, pch=16)
 model2<-clogit(resp~periodo+drogaA+drogaB+gpidade+strata(obs),data=cross4)
 model2
 summary(model2)
 plot(model2$residuals, pch=16)
#
 model3<-clogit(resp~periodo+drogaA+drogaB+strata(obs),data=cross4)
 model3
 summary(model3)
#
# Testando Ho: tau1 = tau2
# --------------------------
#
 model3$var
 vardif<-model3$var[2,2]+model3$var[3,3]-2*(model3$var[2,3])
 teste<-((1.408-0.296)/sqrt(vardif))^2
 teste
 1-pchisq(teste,1)
#
# 15.3 Exemplo: estudo caso-controle
# ==================================
#
 match<-read.table("http://www.est.ufpr.br/~suely/CE215/Dados/match.txt",h=T)
 attach(match)
 require(survival)
 model1<-clogit(cc~hvb+est+hip+id+nest+strata(par),data=match)
 model1
 model2<-clogit(cc~hvb+est+strata(par),data=match)
 model2
 summary(model2)
 plot(model2$residuals, pch=16)
#========================================================================================