#
#   Dados de Constipação - Caroline - META
#
#=================================================================================
#
# No banco de dados foram incluídos 151 recém-nascidos 
# acompanhados nos primeiros 12 meses de vida. Foram previstos um 
# total de 1812 (151*12) medições mensais. Entretanto, sete crianças 
# perderam o acompanhamento em algum momento, resultando em um total 
# de 1751 medidas. Ou seja, tivemos 61 valores perdidos ao longo do 
# acompanhamento de 12 meses.
#
# Este estudo apresenta cinco variáveis respostas: (1) Binárias: Dificuldade 
# para evacuar (difevac), Esforço evacuatório (esfevac), Dor ao evacuar (dorevac), 
#  Tipo de fezes (Tipofezcod) e 
# (2) Contagem: Freqüência evacuatória/semana (freqevac). 
#
# Variável temporal: idade (em dias (idadedias) ou em meses (idademeses)).
#
# Foram registradas uma covariáveis fixa (sexo) e as seguintes 
# dependentes do tempo:
#  1- Aleitamento maternal (2): uma com dois níveis (materno e outros) e 
#     outra com 3 (materno, misto e artificial (mamadeira)).
#  2- Dieta (binária:s/n): cereais; frutas; vegetais, legumin (leguminosas), leitev (leite de vaca),
#      leitec (leite de cabra), formpart (fórmula de partida), formsegui (fórmula de seguimento), 
#      carnes, açudoces (açúcar e doces), sangue nas fezes, dirreia, constag (constipação).
#
#==================================================================================
#
#   Nesta análise preliminar vamos considerar duas respostas: (1) Binária: Dificuldade 
#           para evacuar (difevac) e (2) Contagem: Freqüência evacuatória/semana (freqevac)
#           e cinco covariáveis: (1) idade (idademeses); (2) sexo (0-menina e 1- menino); (3) consumo de cereais (cereais)
#           (4) consumo de carne (carnes) e (5) aleitamento materno (aleitamcod3(0-peito, 1-misto e 2-artificial)).   
#
#
 require(nlme)
 require(lme4)
 require (geepack)
 dados<-read.table('banco_dados.txt',dec=",",header=T)
 names(dados)
# [1] "ID"          "Paciente"    "sexo"        "idadedias"   "idademeses" 
# [6] "aleitamcod2" "aleitamcod3" "cereais"     "frutas"      "vegetais"   
# [11] "legumin"     "leitev"      "leitec"      "formpart"    "formsegui"  
# [16] "carnes"      "açudoces"    "sangue"      "dirreia"     "constag"    
# [21] "difpevac"    "esfevac"     "dorevac"     "Tipofezcod"  "freqevac"   
#[26] "tipofeze" 
 head(dados)
 is.numeric(dados$idademeses)
 is.numeric(dados$freqevac)
 summary(dados)
 dim(dados)
#
#  Observe: (1) que o Banco já está no formato longo e (2) Paciente é
#              a variável indicadora de paciente
#==============================================================================
# Análise Descritiva
#
 length(unique(dados$Paciente)) # qtos individuos no estudo 
 table(dados$Paciente) # qtos individuos por tempo
 as.data.frame(table(dados$Paciente)) #mostra qtas vezes cada individuo aparece 
 dados[dados$Paciente==108,] # Paciente 108 somente uma medida com 53 dias
#
 max(table(dados$Paciente)) #numero maximo de medidas por individuo
 min(table(dados$Paciente)) #numero minimo de medidas por individuo
# Qtos individuos por sexo? 
 tam<-function(x) {length(unique(x))}
 tam(dados$sexo)
 tapply(dados$Paciente,dados$sexo, FUN=tam) #mostra qtos indivíduos por sexo
#  
#  Vamos criar uma idade discreto somente para termos uma ideia
#    preliminar do comportamento das variáveis ao longo do acompanhamento
#    
  dados$idade<-round(dados$idademeses)
  table(dados$idade)
  tam(dados$idade)
#
# Resposta binária: difpevac
#
 y=dados$difpevac
 x<-table(y,dados$idade)
 x
 prop<-x[2,]/(x[1,]+x[2,]) 
 prop
 eixox<-c(1:13)
 plot(eixox,prop,type="l",xlab="idade (meses)", ylab="proporção")
 title(main="Dif. Evacuar")
# Gráfico na escala do logit
 plot(eixox,log(prop/(1-prop)),type="l",xlab="idade (meses)", ylab="logit")
 title(main="Dif. Evacuar (escala logit)")
#
#  Estratificando por sexo
#
 dados1<-subset(dados,sexo==0)
 dados2<-subset(dados,sexo==1)
 par(mfrow=c(1,2))
# 
 y1=dados1$difpevac
 x1<-table(y1,dados1$idade)
 x1
 prop1<-x1[2,]/(x1[1,]+x1[2,]) 
 prop1
#
 y2=dados2$difpevac
 x2<-table(y2,dados2$idade)
 x2
 prop2<-x2[2,]/(x2[1,]+x2[2,]) 
 prop2
#
 eixox<-c(1:13)
 plot(lowess(eixox,prop1),type="l",xlab="idade (meses)", ylab="proporção",ylim=c(0,.40))
 title(main="Dif. Evacuar")
 lines(lowess(eixox,prop2),type="l",lty=2)
 legend(2, .35, legend = c("meninas","meninos"), 
 lty = c(1,2), bty = "n", cex = 1)
# Gráfico na escala do logit
 plot(eixox,log(prop1/(1-prop1)),type="l",xlab="idade (meses)", ylab="logit")
 title(main="Dif. Evacuar (escala logit)")
 lines(eixox,log(prop2/(1-prop2)),type="l",xlab="idade (meses)", ylab="logit",lty=2)
 title(main="Dif. Evacuar (escala logit)")
 legend(2, -2, legend = c("meninas","meninos"), 
 lty = c(1,2), bty = "n", cex = 1)
#
#
#    GEE
#
#  Simetria Composta
# 
 out<-geeglm(difpevac~sexo+legumin+idademeses+sangue+constag+cereais+factor(aleitamcod3)+sexo:idademeses+
      carnes+factor(aleitamcod3)*idademeses, 
      family=binomial(link="logit"), data=dados, id=Paciente,
      corstr ="exchangeable", 
      std.err="san.se")
  summary(out)
 out1<-geeglm(difpevac~sexo+legumin+idademeses+sangue+constag+cereais+factor(aleitamcod3)+
      carnes, 
      family=binomial(link="logit"), data=dados, id=Paciente,
      corstr ="exchangeable", 
      std.err="san.se")
  summary(out1)
 out2<-geeglm(difpevac~sexo+legumin+idademeses+sangue+constag+factor(aleitamcod3), 
      family=binomial(link="logit"), data=dados, id=Paciente,
      corstr ="exchangeable", 
      std.err="san.se")
  summary(out2)
#
#
#   Modelo Misto
#
#
 m1 <- glmer(difpevac~(1|Paciente)+sexo+aleitamcod2+legumin+idademeses+sangue+
 constag,family=binomial,data=dados)
 summary(m1)
#
#
# Resposta contagem: freqevac
#
 y=dados$freqevac
 min(y)
 max(y)
#
# Gráfico na escala log
 plot(dados$idade,log(y),type="l",xlab="idade (meses)", ylab="log contagem",lwd=0.2)
 title(main="Dif. Evacuar (escala logit)")
 lines(lowess(dados$idade,log(y)),col=2, lwd=3)
#
#
#    GEE
#
#  Simetria Composta
# 
 out<-geeglm(difpevac~sexo+aleitamcod2+legumin+idademeses+sangue+constag, 
      family=binomial(link="logit"), data=dados, id=Paciente,
      corstr ="exchangeable", std.err="san.se")
  summary(out)
#
#
#   Modelo Misto
#
#
 m1 <- glmer(difpevac~(1|Paciente)+sexo+aleitamcod2+legumin+idademeses+sangue+
 constag,family=binomial,data=dados)
 summary(m1)

#