#=================================================
# Influência da Menarca no Ganho de Peso Corporal
#  Fitzmaurice et al. (2004) p. 217
#=================================================
# 
# 162 garotas do MIT Grown and Development Study
#
#  Os dados representam um subconjunto do banco de dados
#
#  No início do estudo todas as garotas estavam na pré-menarca
#  e não eram obesas. É conhecido que ganho de peso ocorrem
#  antes ou próximo da menarca. Acredita-se que o ganho de peso
# estacione quatro anos após a menarca. 
#
#  Objetivo: estudar o comportamento do peso de garotas a partir
#            da menarca. 
#
#  Variáveis: time: tempo após a menarca (pode ser negativo)
#             idade na menarca
#             pbf: porcentagem de gordura corporal 
#    
#==============================================
#       Carregando pacotes exigidos
#=============================================
 require(nlme)
 source("longit.r") 
#
#================================================
#        Leitura dos dados 
#==================================================
 dados<-read.table("fat.txt",header=T)
 attach(dados)
 head(dados)
#
 length(unique(ID)) # qtos individuos no estudo 
 as.data.frame(table(ID)) #mostra qtas vezes cada individuo aparece
 max(table(ID)) #numero maximo de medidas por individuo
 min(table(ID)) #numero minimo de medidas por individuo
#
#=============================================== 
#     Gráfico de perfis
#================================================
 datalong = as.ldframe(dados,ID,time)
 plot(datalong,pbf)
 lines(lowess(time,pbf),col="red",lwd=2)
 abline(v=0,lty=2)
#
#================================================== 
#       Ajuste de modelos
#===================================================
#
#=====================================================
#               Modelo 1
#   Modelo com knot no tempo=0 (duas inclinações: 
#    antes e pós menarca
#=====================================================
#
 time0 = ifelse(time>0,time,0)
# time0
 out<-lme(pbf~time+time0, random= ~time+time0|ID,data=dados)
 summary(out)
 names(out)
 plot(out)
#
#=========================================================
#         Análise de resíduos
#==========================================================
#
#   Resíduos não-transformados
 res = out$residuals[,1]
 pred = out$fitted[,1]
 length(res)
#
#   Resíduos transformados
 g11 = 6.777960^2
 g22 = 1.277076^2
 g33 = 1.658200^2
 g12 = 0.292*6.777960*1.277076
 g13 = -0.544*6.777960*1.658200
 g23 = -0.827*1.277076*1.658200
#
 G = matrix(c(g11,g12,g13,g12,g22,g23,g13,g23,g33),3,3,byrow=T)
 sigma2 = (out$sigma)^2
#
 Z1 = rep(1,1049)
 Z2 = time
 Z3 = time0
 Z = cbind(ID,Z1,Z2,Z3)
 resID = cbind(res,ID)
 predID = cbind(pred,ID)
 timeID = cbind(time,ID)
 rest = NULL
 predt = NULL
 timet = NULL
 for (i in 1:162){
 covYi = Z[ID==i,2:4]%*%G%*%t(Z[ID==i,2:4])+sigma2*diag(sum(ID==i))
 rest = c(rest,as.vector(solve(chol(covYi))%*%resID[ID==i,1]))
 predt = c(predt,as.vector(solve(chol(covYi))%*%predID[ID==i,1]))
 timet = c(timet,as.vector(solve(chol(covYi))%*%timeID[ID==i,1]))
 }
#
#   Histogramas com curva normal
 par(mfrow=c(1,2))
 x = seq(min(res)-1,max(res)+1,by=1)
 hist(res,prob=T,xlab="Resíduos",ylab="Densidade",main="")
 lines(x,dnorm(x,mean=mean(res),sd=sd(res)))
 x = seq(min(rest)-1,max(rest)+1,by=0.05)
 hist(rest,prob=T,xlab="Resíduos transformados",ylab="Densidade",main="")
 lines(x,dnorm(x,mean=mean(rest),sd=sd(rest)))
# 
#   Gráfico de Probabilidade Normal dos Resíduos
 par(mfrow=c(1,2))
 qqnorm(res,xlab="Quantis da Normal Padrão",ylab="Quantis dos Resíduos",main="")
 qqline(res)
 qqnorm(rest,xlab="Quantis da Normal Padrão",ylab="Quantis dos Resíduos Transformados",main="")
 qqline(rest)
# 
#    Gráfico de dispersão dos Resíduos versus Valores ajustados
 par(mfrow=c(1,2))
 plot(pred,res,xlab="Valores preditos",ylab="Resíduos")
 abline(h=0)
 lines(lowess(pred,res),type="l",lwd=2,col="red")
 plot(predt,rest,xlab="Valores preditos transformados",ylab="Resíduos transformados")
 abline(h=0)
 lines(lowess(predt,rest),type="l",lwd=2,col="red")
# 
#    Gráfico de dispersão dos Resíduos versus Variável Tempo
 par(mfrow=c(1,2))
 plot(time,res,xlab="Tempo em relação à menarca",ylab="Resíduos")
 abline(h=0)
 lines(lowess(time,res),type="l",lwd=2,col="red")
 plot(timet,rest,xlab="Tempo transformado",ylab="Resíduos transformados")
 abline(h=0)
 lines(lowess(timet,rest),type="l",lwd=2,col="red")
# 
 plot(Variogram(out,form=~time,resType="n", robust=T))
#
#=============================================================
#      Modelo 2
#==============================================================
 time2 = ifelse(time>0,time^2,0)
 out2<-lme(pbf~time+time0+time2, random= ~time+time0+time2|ID,data=dados)
 summary(out2)
 plot(out2)
 intervals(out2)
 out01<-update(out,method="ML")
 out21<-update(out2,method="ML")
 anova(out01,out21)
#
#===================================================================
#        Análise de resíduos
#====================================================================
#
# Resíduos não-transformados
 res = out2$residuals[,1]
 pred = out2$fitted[,1]
 length(res)
#   Resíduos transformados
 g11 = 6.9324298^2
 g22 = 1.3162990^2
 g33 = 2.3171645^2
 g44 = 0.3423537^2
 g12 = 0.332*6.9324298*1.3162990
 g13 = -0.597*6.9324298*2.3171645
 g14 = 0.273*6.9324298*0.342353
 g23 = -0.503*1.3162990*2.3171645
 g24 = -0.385*1.3162990*0.3423537
 g34 = -0.554*2.3171645*0.3423537
#
 G = matrix(c(g11,g12,g13,g14,g12,g22,g23,g24,g13,g23,g33,g34,g14,g24,g34,g44),4,4,byrow=T)
 sigma2 = (out2$sigma)^2
#
 Z1 = rep(1,1049)
 Z2 = time
 Z3 = time0
 Z4 = time2
 Z = cbind(ID,Z1,Z2,Z3,Z4)
 resID = cbind(res,ID)
 predID = cbind(pred,ID)
 timeID = cbind(time,ID)
 rest = NULL
 predt = NULL
 timet = NULL
 for (i in 1:162){
 covYi = Z[ID==i,2:5]%*%G%*%t(Z[ID==i,2:5])+sigma2*diag(sum(ID==i))
 rest = c(rest,as.vector(solve(chol(covYi))%*%resID[ID==i,1]))
 predt = c(predt,as.vector(solve(chol(covYi))%*%predID[ID==i,1]))
 timet = c(timet,as.vector(solve(chol(covYi))%*%timeID[ID==i,1]))
 }
#
#   Histogramas com curva normal
 par(mfrow=c(1,2))
 x = seq(min(res)-1,max(res)+1,by=1)
 hist(res,prob=T,xlab="Resíduos",ylab="Densidade",main="")
 lines(x,dnorm(x,mean=mean(res),sd=sd(res)))
 x = seq(min(rest)-1,max(rest)+1,by=0.05)
 hist(rest,prob=T,xlab="Resíduos transformados",ylab="Densidade",main="")
 lines(x,dnorm(x,mean=mean(rest),sd=sd(rest)))
#
#    Gráfico de Probabilidade Normal dos Resíduos
 par(mfrow=c(1,2))
 qqnorm(res,xlab="Quantis da Normal Padrão",ylab="Quantis dos Resíduos",main="")
 qqline(res)
 qqnorm(rest,xlab="Quantis da Normal Padrão",ylab="Quantis dos Resíduos Transformados",main="")
 qqline(rest)
#
#    Gráfico de dispersão dos Resíduos versus Valores ajustados
 par(mfrow=c(1,2))
 pred = out$fitted[,1]
 plot(pred,res,xlab="Valores preditos",ylab="Resíduos")
 abline(h=0)
 lines(lowess(pred,res),type="l",lwd=2,col="red")
 plot(predt,rest,xlab="Valores preditos transformados",ylab="Resíduos transformados")
 abline(h=0)
 lines(lowess(predt,rest),type="l",lwd=2,col="red")
#
#     Gráfico de dispersão dos Resíduos versus Variável Tempo
 par(mfrow=c(1,2))
 plot(time,res,xlab="Tempo em relação à menarca",ylab="Resíduos")
 abline(h=0)
 lines(lowess(time,res),type="l",lwd=2,col="red")
 plot(timet,rest,xlab="Tempo transformado",ylab="Resíduos transformados")
 abline(h=0)
 lines(lowess(timet,rest),type="l",lwd=2,col="red")
#
 plot(Variogram(out2,form=~time+time0+time2,resType="n", robust=T))
#
# Verificando a necessidade de efeito aleatório no termo quadrático
 out3<-lme(pbf~time+time0+time2, random= ~time+time0|ID,data=dados)
 anova(out2,out3)
#
#==========================================================================
# Modelos Errados!!!!!!
#============================================================================
#     1- Simetria Composta
 out1<-lme(pbf~time+time0+time2, random= ~1|ID,data=dados) #simetria composta
 summary(out1) 
 plot(out1)
 plot(Variogram(out1,form=~time,resType="n", robust=T))
#
#     2- Uma única reta
 out0<-lme(pbf~time, random= ~time|ID,data=dados)
 summary(out0)
 names(out0)
 plot(out0)
# Resíduos não-transformados
 res = out0$residuals[,1]
 pred = out0$fitted[,1]
 length(res)
#
 x = seq(min(res)-1,max(res)+1,by=1)
 hist(res,prob=T,xlab="Resíduos",ylab="Densidade",main="")
 lines(x,dnorm(x,mean=mean(res),sd=sd(res)))
#
 qqnorm(res,xlab="Quantis da Normal Padrão",ylab="Quantis dos Resíduos",main="")
 qqline(res)
#
 plot(pred,res,xlab="Valores preditos",ylab="Resíduos")
 abline(h=0)
 lines(lowess(pred,res),type="l",lwd=2,col="red")
#
 plot(time,res,xlab="Tempo em relação à menarca",ylab="Resíduos")
 abline(h=0)
 lines(lowess(time,res),type="l",lwd=2,col="red")
#
 plot(Variogram(out0,form=~time,resType="n", robust=T))
#======================================================================
#  Modelo Marginal com Estrutura Espacial
#=======================================================================
#   Exponencial 
 modelo1<-gls(pbf~time+time0+time2,correlation=corExp(form=~time|ID))
 summary(modelo1)
 plot(modelo1)
 plot(Variogram(modelo1,form=~time,resType="n", robust=T))
#
#=========================================================================
# Gráfico de Observados vs Preditos
#=========================================================================
 plot(out2,pbf~fitted(.),adj=-0.3) # predito vs obs.
 out22<-data.frame(fitted,dados$ID,dados$time)
 datalong1 = as.ldframe(out22,ID,time)
 par(mfrow=c(1,2))
 plot(datalong,pbf)
 lines(lowess(time,pbf),col="red",lwd=2)
 lines(lowess(time,fitted),col="blue",lwd=2)
 abline(v=0,lty=2)
#
 summary(out2)
#============================================================================