#
 #=========================================================
 #           Leitura e Preliminares dos Dados 
 #=========================================================
 #
 # Fonte: Applied Longitudinal Analysis - 2004 FLW
 # http://www.biostat.harvard.edu/~fitzmaur/ala/
 #
 # Subamostra (N=100) de Dados de Níveis de chumbo no sangue  
 # no estudo de tratamento de crianças expostas a chumbo.
 #
 # Source: Data courtesy of Dr. George G. Rhoads (Chair, TLC 
 # Steering Committee).
 #
 # Referência: Treatment of Lead-exposed Children (TLC) Trial Group. (2000).
 # Safety and Efficacy of Succimer in Toddlers with Blood Lead Levels of 
 # 20-44 µg/dL.  Pediatric Research, 48, 593-599.
 #
 # Descrição: 
 # 
 # O estudo envolvendo o Tratamento de Crianças Expostas ao chumbo (TLC)
 # foi clínico aleatorizado para placebo e um agente succimer em crianças
 # com níveis de chumbo no sangue entre 20-44 micrograms/dL.
 # Estes dados consistem de quatro medidas repetidas de níveis de chumbo
 # no sangue obtidos na baseline (semana 0), semana 1, semana 4 e semana 6
 # em 100 crianças aleatoriamente alocadas entre tratamento e placebo.
 #
 # Lista de variáveis no arquivo de dados (chumbo.txt): 
 #
 # Cada linha contém as seguintes 6 variáveis (formato largo):
 #
 # 1- ID, 
 # 2- Groupo, 
 # 3- Week0, 
 # 4- Week1, 
 # 5- Week4, 
 # 6- Week6. 
 #
 # Resumo: estudo clínico aleatório, balanceado: 50 indivíduos em 
 # cada grupo avaliados em 4 momentos, objetivo (comparar os grupos
 # avaliar a evolução no tempo) e resposta contínua. 
 #
 source("longit.r") #carrega a biblioteca longit.r
 #require(growth)
 dir() # verifica os arquivos no diretório
 dados<-read.table("chumbo.txt",header=T)
 attach(dados)
 names(dados)
 dim(dados)
 dados[1:15,]
 fix(dados) #Planilha de Dados
 # Descritiva dos dados
 summary(Week0) 
 summary(Week1)
 summary(Week4)
 summary(Week6)
 boxplot(Week0,Week1,Week4,Week6,ylab="pb (mg/dl)",xlab="Semana")
 axis(1, 1:4, c(0,1,4,6))
 ### Cuidado pois o Boxplot não considera a estrutura longitudinal ####
 datawide<-data.frame(dados$ID,dados$Groupo,dados$Week0,dados$Week1,
 dados$Week4,dados$Week6)
 cor(datawide[,3:6]) # matriz de correlação dos dados #deveria ser feito com resíduos
 #
 #========================================================
 #         Manipulando Banco de Dados Longitudinais
 #            Banco Longo (uma medida por linha)
 #=========================================================
 #
 # Transformando o banco wide para longo
 datalong<-reshape(data=datawide,direction="long", idvar="ID", 
 varying = list(names(datawide)[3:6]), v.names="chumbo", time= c(0,1,4,6),
 timevar="tempo")
 datalong$Groupo<-datalong$dados.Groupo
 attach(datalong)
 dim(datalong)
 names(datalong)
 datalong[1:15,]
#  
 length(unique(ID)) # qtos individuos no estudo 
 table(tempo) # qtos individuos por tempo
 table(tempo,Groupo) # qtos individuos por tempo e dieta
 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
 # Qtos individuos por tratamento? 
 tam <- function(x) {
        length(unique(x))
        }
 tapply(ID,Groupo, FUN=tam) #mostra qtos indivíduos por dieta
 #
 #==========================================================
 #          Estatísticas Descritivas
 #==========================================================
 #
 tapply(chumbo,list(Groupo), mean) # media por grupo e tempo
 tapply(chumbo,list(Groupo), sd) # dp por grupo e tempo
 tapply(chumbo,list(Groupo,tempo), mean)
 plot(Groupo,chumbo)   
 plot(tempo,chumbo)
 #
 # Os comandos abaixo somente funcionam na com o longit.r
 # Os gráficos de perfis mais abaixo não necessitam do longit.r
 #
 datalong1 <- as.ldframe(datalong, ID, tempo)  
 class(datalong)
 plot(datalong1, chumbo, color = Groupo, cols = c("blue", "red"))
 legend(2, 60, legend = paste("Tipo", unique(datalong1$Groupo)), 
 col = c("blue","red"), lty = 1, bty = "n", cex = 0.8)
 #
 # Gráfico para 6 indivíduos escolhidos aleatoriamente
 #
 plot(datalong1, chumbo, n.full = 6, color = Groupo, cols = c("blue", "red"))
 legend(2,60, legend = paste("Tipo", unique(datalong1$Groupo)), 
 col = c("blue","red"), lty = 1, bty = "n", cex = 0.8)
 #
 # Gráfico de perfis excluindo a baseline
 #
 plot(datalong[which(datalong1$Time!=0),], chumbo, color = Groupo, cols = c("blue", "red"))
 legend(2,60, legend = paste("Tipo", unique(datalong1$Groupo)), 
 col = c("blue","red"), lty = 1, bty = "n", cex = 0.8)
 #
 # Identifica outliers
 #
 outliers<-identify(datalong1, chumbo) 
 datalong[outliers,] # armazenou os outliers
 datalong2<-datalong1[which(datalong1$Subject !=40),] #exclui o ind. 40
 #
 #==================================================================
 #          Graficos de Perfis
 #==================================================================
 #
 teste<-as.numeric(Groupo)
 teste<-teste*2 #(muda de cores (1,2) para (2,4))
 #
 # function(choice,xl=range(time),yl=range(y))
 # {# Plots profiles for subjects in vector choice. 
 m=length(unique(ID))
 n=length(chumbo)/m
 xl=range(tempo)
 yl=range(chumbo)  
 choice=ID  
 j=ID==choice[1]
 plot(tempo[j],chumbo[j],xaxt="n",
 type="l",col = teste[j],ylim=yl,xlim=xl,ylab="chumbo (mg/dl)",xlab="Mês",
 main="Efeito de Tratamento")
 #
 if(length(choice)>1){
 for(it in 2:length(choice)){
 j=ID==choice[it]
 lines(tempo[j],chumbo[j],col = teste[j])
 }
 axis(1, c(0,1,4,6))
 }
 cat("\n")
        }
 #Fazendo um alisamento no grafico de perfis 
 plot(lowess(tempo,chumbo),type="l",xaxt="n",ylab="chumbo",xlab="mês")  
      axis(1, c(0,1,4,6)
 plot(lowess(Groupo,chumbo),type="l",xaxt="n",ylab="chumbo",xlab="Grupo")  
      axis(1, 1:2, c(0,1))
 #
 #
 #==================================================================
 #          Graficos de Perfis por Grupo
 #==================================================================
 #
 attach(datalong)
 placebo<-subset(datalong,Groupo=="P")
 plot(lowess(placebo$tempo,placebo$chumbo),
 xaxt="n", col = "blue",
 type="l",lty=1,
 ylim=range(c(10,30)),ylab="Pb",xlab="mês") 
 trat<-subset(datalong,Groupo=="A")
 lines(lowess(trat$tempo,trat$chumbo),lty=2,type="l",col="red")
 axis(1, c(0,2,4,6))
 legend(1,30, legend = paste("Tipo", unique(datalong1$Groupo)), 
 col = c("blue","red"), lty = 1, bty = "n", cex = 0.8)
 #
 #===============================================================
 #       Explorando os Perfis Médios por Grupo
 #===============================================================
 # Gráfico de Interação
 attach(datalong) 
 interaction.plot(tempo, Groupo, chumbo)
 # Box-plots por Dieta
 boxplot(chumbo ~ tempo + Groupo, notch = T)
 abline(v = 4.5, lty = 2)
 abline(v = 10.5, lty = 2)
 title("Dispersão do chumbo por grupo")
 #===============================================================
 #       Explorando a Estrutura de Variância
 #===============================================================
 #          
 perfil.std(y = datawide[, 3:6], id = datawide$ID, group = datawide$Groupo)
 interaction.plot(tempo, Groupo, chumbo, trace.label = "Groupo",
 xlab = "Mês", ylab = "Variância Médio", fun = var, fixed = T)
 #
 #====================================================================
 #      Explorando a Estrutura de Dependência dos Tempos
 #====================================================================
 # 
 require(car)
 attach(datawide)
 dados.std <- perfil.std(y = datawide[, 3:6], id = ID, group = Groupo, 
 flag = FALSE)
 scatterplot.matrix(dados.std[, 2:5], diagonal = "histogram", 
 smooth = F)
#graphics.off()
#=====================================================================
# Explorando a Estrutura de Variância - Deveria utilizar resíduos
#     Vamos utilizar uma análise para cada grupo
#======================================================================
 cor(dados.std[, 2:5]) # análise geral (sem estratificar por grupo)
#
 datawideP <-subset(datawide,Groupo=="P")
 dados.std <- perfil.std(y = datawideP[, 3:6], id = ID, 
 flag = FALSE)
 scatterplot.matrix(dados.std[, 2:5], diagonal = "histogram", 
 smooth = F)
 cor(dados.std[, 2:5])
#
 datawideA <-subset(datawide,Groupo=="A")
 dados.std <- perfil.std(y = datawideA[, 3:6], id = ID, 
 flag = FALSE)
 scatterplot.matrix(dados.std[, 2:5], diagonal = "histogram", 
 smooth = F)
 cor(dados.std[, 2:5])
#
=================================================================
# Explorando a Estrutura de Variância - Utilizando resíduos 
#   Modelo Maximal --> Mínimos Quadrados Ordinários  
#======================================================================
#
  attach(datalong) 
 fit<-lm(chumbo ~  factor(Groupo)*factor(tempo)) # 8 parâmetros
 res0<-subset(fit$residuals,tempo==0)
 res1<-subset(fit$residuals,tempo==1)
 res4<-subset(fit$residuals,tempo==4)
 res6<-subset(fit$residuals,tempo==6)
 res<-data.frame(res0,res1,res4,res6)
 cor(res) # possível indicação da estrutura de simetria composta
 cov(res) # possível heterocedasticidade
#
#==========================================================================
#   Modelos de Efeitos Fixos com Estrutura de Covariância: simetria
#   composta, AR1, não-estruturada.
#       EMVR
#===========================================================================
#
 require(lme4)
 require(nlme)
## MODELO 1
## Modelo linear de efeito fixo (sem intercepto)
#      mais simples de comparar os grupos em cada tempo
#  
 out0<-lm(chumbo ~ factor(tempo) + factor(tempo):factor(Groupo) -1,data=datalong)
 summary(out0)
 plot(out0)
 out<-gls(chumbo ~ factor(tempo) + factor(tempo):factor(Groupo) -1,
 correlation=corCompSymm (form=~1|ID),data=datalong) 
 summary(out)
 plot(out)
 out1<-gls(chumbo ~ factor(tempo) + factor(tempo):factor(Groupo) -1,
 correlation=corAR1 (form=~1|ID),data=datalong) 
 summary(out1)
 plot(out1)
 out2<-gls(chumbo ~ factor(tempo) + factor(tempo):factor(Groupo) -1,
 correlation=corSymm (form=~1|ID),data=datalong) 
 summary(out2)
 plot(out2)
#
## MODELO 2
## Modelo linear de efeito fixo (com intercepto)
#  
 out0<-lm(chumbo ~ factor(tempo)*factor(Groupo),data=datalong)
 summary(out0)
 out<-gls(chumbo ~ factor(tempo)*factor(Groupo),
 correlation=corCompSymm (form=~1|ID),data=datalong) 
 summary(out)
 out1<-gls(chumbo ~ factor(tempo)*factor(Groupo),
 correlation=corAR1 (form=~1|ID),data=datalong) 
 summary(out1)
 out2<-gls(chumbo ~ factor(tempo)*factor(Groupo),
 correlation=corSymm (form=~1|ID),data=datalong) 
 summary(out2)

##=================================================================
#     GEE - MODELO MAXIMAL: TEMPO COMO FATOR
#==================================================================
#   
 require(geepack)
 dados1<-datalong[order(datalong$ID),]
 dados1[1:15,]
 attach(dados1)
#
#  Modelo 1 - Sem Intercepto
# 
# Mínimos Quadrados Ordinários (Parâmetros: 8 + 1)
 out<-lm(chumbo ~ factor(tempo)+ factor(tempo)*factor(Groupo)-1,data=dados1) 
 summary(out)
# Independência (Parâmetros: 8 + 1)
 out1<-geeglm(chumbo ~ factor(tempo)+ factor(tempo)*factor(Groupo)-1, 
 id = ID, data = dados1, corstr = "independence")
 summary(out1)
 # AR(1) # Não é válido pois não é igualmente espaçado
 out2<-geeglm(chumbo ~ factor(tempo)+ factor(tempo)*factor(Groupo)-1, 
 id = ID, data = dados1, corstr="corar1")
 summary(out2)
# Simetria Composta (Parâmetros: 8 + 2)
 out3<-geeglm(chumbo ~ factor(tempo)+ factor(tempo)*factor(Groupo)-1, 
 id = ID, data = dados1, corstr = "exchangeable")
 summary(out3)
# Nao Estruturada (Parametros: 8 + 7)
 out4<-geeglm(chumbo ~ factor(tempo)+ factor(tempo)*factor(Groupo)-1, 
 id = ID, data = dados1, corstr = "unstructured")
 summary(out4)
# Avaliando os resíduos
 fit<-out4
 res0<-subset(fit$residuals,tempo==0)
 res1<-subset(fit$residuals,tempo==1)
 res4<-subset(fit$residuals,tempo==4)
 res6<-subset(fit$residuals,tempo==6)
 res<-data.frame(res0,res1,res4,res6)
 cor(res) # possível indicação da estrutura de simetria composta
 cov(res) # possível heterocedasticidade
#
#   Modelo 2 (com intercepto)
#
# Mínimos Quadrados Ordinários (Parâmetros: 8 + 1)
 out<-lm(chumbo ~ factor(tempo)*factor(Groupo),data=dados1) 
 summary(out)
# Independência (Parâmetros: 8 + 1)
 out1<-geeglm(chumbo ~ factor(tempo)*factor(Groupo), id = ID, data = dados1, 
 corstr = "independence")
 summary(out1)
 # AR(1) # Não é válido pois não é igualmente espaçado
 out2<-geeglm(chumbo ~ factor(tempo)*factor(Groupo), id = ID, data = dados1, 
 corstr="corar1")
 summary(out2)
# Simetria Composta (Parâmetros: 8 + 2)
 out3<-geeglm(chumbo ~ factor(tempo)*factor(Groupo), id = ID, data = dados1, 
 corstr = "exchangeable")
 summary(out3)
# Nao Estruturada (Parametros: 8 + 7)
 out4<-geeglm(chumbo ~ factor(tempo)*factor(Groupo), id = ID, data = dados1, 
 corstr = "unstructured")
 summary(out4)
# Avaliando os resíduos
 fit<-out4
 res0<-subset(fit$residuals,tempo==0)
 res1<-subset(fit$residuals,tempo==1)
 res4<-subset(fit$residuals,tempo==4)
 res6<-subset(fit$residuals,tempo==6)
 res<-data.frame(res0,res1,res4,res6)
 cor(res) # possível indicação da estrutura de simetria composta
 cov(res) # possível heterocedasticidade
#
# Mínimos Quadrados Ordinários (Parâmetros: 8 + 1)
 out<-lm(chumbo ~ factor(tempo)*factor(Groupo),data=dados1) 
 summary(out)
# Independência (Parâmetros: 8 + 1)
 out1<-geeglm(chumbo ~ factor(tempo)*factor(Groupo), id = ID, data = dados1, 
 corstr = "independence")
 summary(out1)
 # AR(1) # Não é válido pois não é igualmente espaçado
 out2<-geeglm(chumbo ~ factor(tempo)*factor(Groupo), id = ID, data = dados1, 
 corstr="corar1")
 summary(out2)
# Simetria Composta (Parâmetros: 8 + 2)
 out3<-geeglm(chumbo ~ factor(tempo)*factor(Groupo), id = ID, data = dados1, 
 corstr = "exchangeable")
 summary(out3)
# Nao Estruturada (Parametros: 8 + 7)
 out4<-geeglm(chumbo ~ factor(tempo)*factor(Groupo), id = ID, data = dados1, 
 corstr = "unstructured",tol = 0.00001, maxiter = 50)
 summary(out4)
# Avaliando os resíduos
 fit<-out4
 res0<-subset(fit$residuals,tempo==0)
 res1<-subset(fit$residuals,tempo==1)
 res4<-subset(fit$residuals,tempo==4)
 res6<-subset(fit$residuals,tempo==6)
 res<-data.frame(res0,res1,res4,res6)
 cor(res) # possível indicação da estrutura de simetria composta
 cov(res) # possível heterocedasticidade
#

#=================================================================
#     GEE - MODELO: TEMPO contínuo (linear e quadrático) e 
#                   termo de Interação tempo*FATOR
#==================================================================
 library(MASS)
 library(geepack)
 library(gee)
 dados1<-datalong[order(ID),]
 dados1[1:15,]
#
# Centrando na média para evitar multicolinearidade
  dados1[,3]<-dados1[,3]-mean(dados1[,3])
  attach(dados1)
#
# Mínimos Quadrados Ordinários (Parâmetros: 5 + 1)
  out<-lm(chumbo ~ tempo + I(tempo^2) + tempo:factor(Groupo) +
  I(tempo^2):factor(Groupo))
  summary(out)
#
# Independência (Parâmetros: 5 + 1)
  out1<-gee(chumbo ~ tempo + I(tempo^2) + 
  tempo:factor(Groupo) +
  I(tempo^2):factor(Groupo), id = ID, data = dados1, 
  corstr="independence")
  summary(out1)
# AR(1)
  out2<-gee(chumbo ~ tempo + I(tempo^2) + tempo:factor(Groupo) +
  I(tempo^2):factor(Groupo), id = ID, data = dados1, 
  corstr="AR-M", Mv=1,tol = 0.00001, maxiter = 50)
  summary(out2)
# Simetria Composta (Parâmetros: 5 + 2)
  out3<-gee(chumbo~tempo + I(tempo^2) + tempo:factor(Groupo) +
  I(tempo^2):factor(Groupo), id = ID, data = dados1, 
  corstr = "exchangeable")
  summary(out3)
# Nao Estruturada (Parametros: 5 + 7)
  out4<-gee(chumbo ~ tempo + I(tempo^2) + tempo:factor(Groupo) +
  I(tempo^2):factor(Groupo), id = ID, data = dados1, 
  corstr="unstructured")
## Aparente problema: um componente da estrutura de correlação estimado
## maior que 1
  summary(out4)
#
#=================================================================
#     Gráfico Preditos vs Observados
#==================================================================
#
 placebo<-subset(datalong,Groupo=="P")
 plot(lowess(placebo$tempo,placebo$chumbo),
 xaxt="n", col = "blue",
 type="l",lty=1,
 ylim=range(c(10,30)),ylab="Pb",xlab="mês") 
 trat<-subset(datalong,Groupo=="A")
 lines(lowess(trat$tempo,trat$chumbo),lty=1,type="l",col="red")
 axis(1, c(0,2,4,6))
 legend(1,30, legend = paste("Tipo", unique(datalong$Groupo)), 
 col = c("blue","red"), lty = 1, bty = "n", cex = 0.8)
#
#  Independência 
 x<-seq(0,6,0.2)
 y1<-25 - (8.2*x) + (1.3*x^2) # escala original
# y1<-18,2 - (0.62*(x-3.75)) + (0.43*(x-3.75)^2) # centrado na média
 lines(lowess(x,y1),lty=2,type="l",col="red")
 y2<-25 - (0.3*x) + (0.01*x^2)
#  y2<-18.2 - (0.75*(x-3.75)) + (0.86*(x-3.75)^2)
 lines(lowess(x,y2),lty=2,type="l",col="blue")
# 
#   AR(1)
 x<-seq(0,6,0.2)
 y1<-25 - (7.1*x) + (1.1*x^2)
 lines(lowess(x,y1),lty=2,type="l",col="red")
 y2<-25 - (0.52*x) + (0.05*x^2)
 lines(lowess(x,y2),lty=2,type="l",col="blue")
#
#  Simetria Composta
 x<-seq(0,6,0.2)
 y1<-25 - (7.8*x) + (1.2*x^2)
 lines(lowess(x,y1),lty=2,type="l",col="red")
 y2<-25 - (0.63*x) + (0.06*x^2)
 lines(lowess(x,y2),lty=2,type="l",col="blue")
#
#  Não Estruturada
 x<-seq(0,6,0.2)
 y1<-24 - (6*x) + (1*x^2)
 lines(lowess(x,y1),lty=2,type="l",col="red")
 y2<-24 - (0.23*x)
 lines(lowess(x,y2),lty=2,type="l",col="blue")
#
#=================================================================
#     GEE - Geração do banco de dados - Linear e quadrático no tempo
# e        AR(1) \rho =0.7 (=1, Simetria Composta)
#         \beta_0 = 5 \beta_1 = 0,25 \beta_2 = 0.5
#==================================================================
 library(MASS)
 library(geepack)
 n=100 #tamanho da amostra
# mu=c(5.25,5.50,5.75,6.00,6.25) #vetor de médias para os tempos lineares
#      \beta_2=0
 mu=c(5.75,7.5,10.25,14,18.75) #vetor de médias para os tempos
 resp=matrix(rep(rep(0,n),4),n,5)
 rho=0.7
 fator=2
 Sigma=fator*rho^(abs(outer(1:5,1:5,FUN="-"))) #matriz de covariância AR-1
 Sigma[1,1]<-2
 Sigma[2,2]<-2
 Sigma[3,3]<-2
 Sigma[4,4]<-2
 Sigma[5,5]<-2
 Sigma
 for(i in 1:n){ #geração das respostas
 resp[i,]=mvrnorm(1,mu, Sigma)
 }
 fitz=as.data.frame(resp)
 colnames(fitz)=c("Y1","Y2","Y3","Y4","Y5")
 D=fitz #banco a ser manipulado
#
# Transformação do banco original em banco longo
# 
 o.id=rep(1:n,5)
 o.resp=c(fitz$Y1,fitz$Y2,fitz$Y3,fitz$Y4,fitz$Y5)
 o.tempo=c(rep(1,n),rep(2,n),rep(3,n),rep(4,n),rep(5,n))
 o.longo=as.data.frame(cbind(o.id,o.resp,o.tempo)) #banco de dados original long
 o.longo<-o.longo[order(o.longo$o.id),]
 colnames(o.longo)=c("ID","Y","tempo")
 o.longo[1:15,]
# Mínimos Quadrados Ordinários (Parâmetros: 5 + 1)
  out<-lm(o.longo$Y ~ o.longo$tempo 
  + I(o.longo$tempo^2)
  )
  summary(out)
# Independência (Parâmetros: 5 + 1)
  out1<-gee(o.longo$Y ~ o.longo$tempo + I(o.longo$tempo^2), id = ID, 
  data = o.longo, corstr="independence")
  summary(out1)
  out1<-geeglm(o.longo$Y ~ o.longo$tempo 
 + I(o.longo$tempo^2)
  , id = ID, data = o.longo, corstr="independence")
  summary(out1)
# AR(1)
  out2<-gee(o.longo$Y ~ o.longo$tempo 
 + I(o.longo$tempo^2)
  , id = ID,  data = o.longo, corstr="AR-M", Mv=1,tol = 0.00001, maxiter = 50)
  summary(out2)
# Simetria Composta (Parâmetros: 5 + 2)
  out3<-gee(o.longo$Y ~ o.longo$tempo 
 + I(o.longo$tempo^2)
  , id = ID,  data = o.longo, corstr = "exchangeable")
  summary(out3)
# Nao Estruturada (Parametros: 5 + 7)
  out4<-gee(o.longo$Y ~ o.longo$tempo 
 + I(o.longo$tempo^2)
  , id = ID,  data = o.longo, corstr="unstructured", tol = 0.00001, 
  maxiter = 50)
  summary(out4)
#
 library(nlme) 
 attach(o.longo)
# gls - AR(1)          
 out<-gls(Y ~ tempo + I(tempo^2), correlation = corAR1(form = ~1 |ID))
# ,method="ML", weights = varIdent(form = ~1| tempo))
 summary(out)
 plot(out)
# 
#==================================================================
#   GLS - EMV
#===================================================================
#
 library(nlme) 
 attach(dados1)
 Groupo<-relevel(Groupo,ref="P")
#
# Não-Estruturada 
#
# Modelo Maximal --> Tempo Categórico          
 out<-gls(chumbo ~ factor(tempo)*factor(Groupo), correlation = corSymm(form = ~1 |ID),
 method="ML", weights = varIdent(form = ~1| tempo))
# Tempo contínuo - componente quadrática
 out<-gls(chumbo ~ tempo + I(tempo^2) + tempo:factor(Groupo) + 
 I(tempo^2):factor(Groupo), correlation = corSymm(form = ~1 |ID),
 method="ML", weights = varIdent(form = ~1| tempo),data=dados1)
 summary(out)
 plot(out)
 # Homocedasticidade
 out0<-gls(chumbo ~ tempo + I(tempo^2) + tempo:factor(Groupo) + 
 I(tempo^2):factor(Groupo), correlation = corSymm(form = ~1 |ID),
 method="ML",data=dados1)
 summary(out0)
 plot(out0)
 anova(out,out0)
#
# Simetria Composta          
# out1<-gls(chumbo ~ factor(tempo)*factor(Groupo), correlation = corCompSymm(form = ~1 |ID),
# method="ML",weights = varIdent(form = ~1|tempo))
 out1<-gls(chumbo ~ tempo + I(tempo^2) + tempo:factor(Groupo) +
  I(tempo^2):factor(Groupo), correlation = corCompSymm(form = ~1 |ID),
 method="ML",weights = varIdent(form = ~1|tempo))
  summary(out1)
 anova(out,out1)
#
#=================================================================
#     Gráfico Preditos vs Observados
#==================================================================
#
#  par(mfrow = c(2, 1))
 placebo<-subset(datalong,Groupo=="P")
 plot(lowess(placebo$tempo,placebo$chumbo),
 xaxt="n", col = "blue",
 type="l",lty=1,
 ylim=range(c(10,30)),ylab="Pb",xlab="mês") 
 trat<-subset(datalong,Groupo=="A")
 lines(lowess(trat$tempo,trat$chumbo),lty=1,type="l",col="red")
 axis(1, c(0,2,4,6))
 legend(1,30, legend = paste("Tipo", unique(datalong$Groupo)), 
 col = c("blue","red"), lty = 1, bty = "n", cex = 0.8)
#
# Simetria Composta 
 x<-seq(0,6,0.2)
 y1<-26.4 - (7.7*x) + (1.2*x^2) # escala original
 lines(lowess(x,y1),lty=2,type="l",col="red")
 y2<-26.4 - (0.96*x) + (0.1*x^2)
 lines(lowess(x,y2),lty=2,type="l",col="blue")
# 
#  Não-Estruturada
 x<-seq(0,6,0.2)
 y1<-25.9 - (6.9*x) + (1.13*x^2) # escala original
 lines(lowess(x,y1),lty=2,type="l",col="red")
 y2<-25 - (0.8*x) + (0.08*x^2)
 lines(lowess(x,y2),lty=2,type="l",col="blue")
# 
#=====================================================================