#
#=========================================================
#           Leitura e Preliminares dos Dados 
#=========================================================
#
# Fonte: Applied Longitudinal Analysis - 2004 FLW
# http://www.biostat.harvard.edu/~fitzmaur/ala/
#
# Subamostra (N=100) de Dados de Niveis de chumbo no sangue  
# no estudo de tratamento de criancas expostas a chumbo.
#
# Source: Data courtesy of Dr. George G. Rhoads (Chair, TLC 
# Steering Committee).
#
# Referencia: 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.
#
# Descricao: 
# 
# O estudo envolvendo o Tratamento de Criancas Expostas ao chumbo (TLC)
# foi clinico aleatorizado para placebo e um agente succimer em criancas
# com niveis de chumbo no sangue entre 20-44 micrograms/dL. Valores de referencia
# de normalidade para chumbo no sangue fica entre 1 e 14 mg/dl.
# Estes dados consistem de quatro medidas repetidas de niveis de chumbo
# no sangue obtidos na baseline (semana 0), semana 1, semana 4 e semana 6
# em 100 criancas aleatoriamente alocadas entre tratamento e placebo.
#
# Lista de variaveis no arquivo de dados (chumbo.txt): 
#
# Cada linha contem as seguintes 6 variaveis (formato largo):
#
# 1- ID, 
# 2- Groupo, 
# 3- Week0, 
# 4- Week1, 
# 5- Week4, 
# 6- Week6. 
#
# Resumo: estudo clinico aleatorio, balanceado: 50 individuos em 
# cada grupo avaliados em 4 momentos, objetivo (comparar os grupos
# avaliar a evolucao no tempo) e resposta continua.
# Objetivos: (1) comparar os dois tratamentos; (2) avaliar o comportamento
# temporal.
#
 rm(list=ls(all=TRUE)) # limpar a console
 require(lme4)
 library(plyr)
 library(nlme)
 library(ggplot2)
 dir() # verifica os arquivos no diretorio
 dados<-read.table("chumbo.txt",header=T)
 names(dados)
 dim(dados)
 dados[1:15,]
#
#  Observe que os dados estao no formato largo, um individuo por linha e as
#   medidas repetidas apresentadas por colunas
# Descritiva dos dados geral
 summary(dados$Week0) 
 summary(dados$Week1)
 summary(dados$Week4)
 summary(dados$Week6)
#
 boxplot(dados$Week0,dados$Week1,dados$Week4,dados$Week6,ylab="pb (mg/dl)",xlab="Semana")
 axis(1, 1:4, c(0,1,4,6))
### Cuidado pois o Boxplot nao considera a estrutura longitudinal por individuo####
#
# Descritiva por grupo
 summary(subset(dados, Groupo=="P")[,3:6]) 
 summary(subset(dados, Groupo=="A")[,3:6]) 
#
#  A medida muda muito pouco no grupo Placebo e no grupo tratado existe uma
#   grande reducaoo da Linha base para a semana 1.
#
 datawide<-data.frame(dados$ID,dados$Groupo,dados$Week0,dados$Week1,
 dados$Week4,dados$Week6)
#
# Vamos precisar mudar o formato do banco para longo para continuar as nossas analises. 
#  
#===================================================================================
#         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
 dim(datalong)
 names(datalong)
 datalong[1:15,]
#  
 length(unique(datalong$ID)) # qtos individuos no estudo 
 table(datalong$tempo) # qtos individuos por tempo
 table(datalong$tempo,datalong$Groupo) # qtos individuos por tempo e dieta
 as.data.frame(table(datalong$ID)) #mostra qtas vezes cada individuo aparece
 max(table(datalong$ID)) #numero maximo de medidas por individuo
 min(table(datalong$ID)) #numero minimo de medidas por individuo
 # Qtos individuos por tratamento? 
 tam <- function(x) {
        length(unique(x))
        }
 tapply(datalong$ID,datalong$Groupo, FUN=tam) #mostra qtos indiv?duos por dieta
#
#=====================================================================================
#          Estatisticas Descritivas (ja realizadas no formato largo)
#=====================================================================================
#
 tapply(datalong$chumbo,list(datalong$Groupo), mean) # media por grupo e tempo
 tapply(datalong$chumbo,list(datalong$Groupo), sd) # dp por grupo e tempo
 tapply(datalong$chumbo,list(datalong$Groupo,datalong$tempo), mean)
 plot(datalong$Groupo,datalong$chumbo)   
 plot(datalong$tempo,datalong$chumbo)
#
#=====================================================================================
#          Graficos de Perfis
#=====================================================================================
#
 p1=ggplot(datalong, aes(x=tempo, y=chumbo, color=Groupo)) +	theme_bw() + 
   geom_line(aes(group=ID)) + theme(legend.position="top") + 
   scale_x_continuous(breaks=unique(datalong$tempo))
 p1
#
#==================================================================
#          Graficos de alisamento por Grupo
#==================================================================
#
 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(datalong$Groupo)), 
 col = c("blue","red"), lty = 1, bty = "n", cex = 0.8)
#
#   Ou se voce quiser usar o ggplot
 p3=ggplot(datalong, aes(x=tempo, y=chumbo, group = Groupo, shape = Groupo)) + theme_bw() +
         stat_summary(fun.y="mean",geom="line", size=1.1, aes(linetype = Groupo)) +
         theme(legend.position="top") + scale_x_continuous(breaks=unique(datalong$tempo)) 
 p3
# 
#   Outra forma alternativa
 interaction.plot(datalong$tempo, datalong$Groupo, datalong$chumbo)
#
#  Outra forma de mostrar as diferencas entre os grupos ao longo do tempo
#
 p2=ggplot(datalong, aes(x=factor(tempo), y=chumbo, fill=Groupo)) + geom_boxplot(notch=TRUE) + 
   theme_bw() + theme(legend.position="top") + stat_summary(fun.y="mean", geom="point", size=2,
    position=position_dodge(width=0.75), color="white", show.legend=FALSE)
 p2
#
#  Boxplot em uma forma alternativa.
#   
 boxplot(datalong$chumbo ~ datalong$tempo + datalong$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 Variancia
#=======================================================================
#
# Avaliacao da correlacao geral e por grupo no formato largo
#
 round(cor(datawide[,3:6]),3) # matriz de correlacao dos dados #deveria ser feito com res?duos
 round(cor(subset(datawide, dados.Groupo=="P")[,3:6]),3) 
 round(cor(subset(datawide, dados.Groupo=="A")[,3:6]),3) 
#
# Existe uma indicacao de similar correlacao em cada grupo, e maior no grupo Placebo.
#
 interaction.plot(datalong$tempo, datalong$Groupo, datalong$chumbo, trace.label = "Grupo",
 xlab = "M?s", ylab = "Variancia Media", fun = var, fixed = T)
#
#========================================================================
# Explorando a Estrutura de Variancia - Utilizando residuos 
#   Modelo Maximal --> Minimos Quadrados Ordinarios  
#========================================================================
#
 fit<-lm(chumbo ~  factor(Groupo)*factor(tempo),data=datalong) # 8 par?metros
 res0<-subset(fit$residuals,datalong$tempo==0)
 res1<-subset(fit$residuals,datalong$tempo==1)
 res4<-subset(fit$residuals,datalong$tempo==4)
 res6<-subset(fit$residuals,datalong$tempo==6)
 res<-data.frame(res0,res1,res4,res6)
 cor(res) # possivel indicacao da estrutura de simetria composta
 cov(res) # possivel heterocedasticidade
#
#==========================================================================
#   Modelos de Efeitos Fixos com Estrutura de Covariancia: simetria
#   composta, nao-estruturada.
#       EMVR
#===========================================================================
#
## MODELO 1
## Modelo linear de efeito fixo (com intercepto) e Homocedasticidade
# 
 datalong$Groupo<-relevel(datalong$Groupo, ref="P") 
 out0<-lm(chumbo ~ factor(tempo)*factor(Groupo),data=datalong) # independente
 summary(out0)
 out1<-gls(chumbo ~ factor(tempo)*factor(Groupo),
 correlation=corCompSymm (form=~1|ID),data=datalong) # simetria composta
 summary(out1)
#  AR(1) nao e indicado por nao ser igualmente espacado
#  out2<-gls(chumbo ~ factor(tempo)*factor(Groupo),
# correlation=corAR1 (form=~1|ID),data=datalong)  # AR(1)
# summary(out2)
 out3<-gls(chumbo ~ factor(tempo)*factor(Groupo),
 correlation=corSymm (form=~1|ID),data=datalong) # Nao estruturada
 summary(out3)
#
#  Estrutura heterocedastica e Nao Estruturada
 out5<-gls(chumbo ~ factor(tempo)*factor(Groupo),correlation=corSymm (form=~1|ID),
 weights = varIdent(form = ~1 | factor(tempo)),data=datalong)
 summary(out5)
 plot(out5)
 anova(out5,out3)  # heterocedastico e significativo 
#
#  Estrutura heterocedastica e Simetria composta
 out4<-gls(chumbo ~ factor(tempo)*factor(Groupo),correlation=corCompSymm (form=~1|ID),
 weights = varIdent(form = ~1 | factor(tempo)),data=datalong)
 summary(out)
 plot(out4)
 anova(out5,out4) # teste significativo --> out5
#
## MODELO 2
## 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)
 out1<-gls(chumbo ~ factor(tempo) + factor(tempo):factor(Groupo) -1,
 correlation=corCompSymm (form=~1|ID),data=datalong) 
 summary(out1)
 plot(out1)
 out3<-gls(chumbo ~ factor(tempo) + factor(tempo):factor(Groupo) -1,
 correlation=corSymm (form=~1|ID),data=datalong) 
 summary(out3)
 plot(out3)
##  Estrutura heterocedastica e Nao Estruturada
 out5<-gls(chumbo ~ factor(tempo)*factor(Groupo)-1,correlation=corSymm (form=~1|ID),
           weights = varIdent(form = ~1 | factor(tempo)),data=datalong)
 summary(out5)
 plot(out5)
#  Estrutura heterocedastica e Simetria composta
 out4<-gls(chumbo ~ factor(tempo)*factor(Groupo) -1,correlation=corCompSymm (form=~1|ID),
 weights = varIdent(form = ~1 | factor(tempo)),data=datalong)
 summary(out4)
 plot(out4)
#
 anova(out5,out4)
#
 ##  Estrutura heterocedastica e Nao Estruturada sem efeito principal de grupo.
 datalong$a<-model.matrix(~factor(datalong$tempo))[,-c(1)]
 datalong$inter<-datalong$a*model.matrix(~factor(datalong$Groupo))[,-c(1)]
 out6<-gls(chumbo ~ factor(tempo) + inter,correlation=corSymm (form=~1|ID),
           weights = varIdent(form = ~1 | factor(tempo)),data=datalong)
 summary(out6)
 names(out6)
 out6$varBeta
 plot(out6)
# 
## MODELO 3
## Modelo duas retas com knot em tempo=1
#
# datalong$tempom<- datalong$tempo - mean(datalong$tempo)
 datalong$tempo1<-ifelse(datalong$tempo>1 ,datalong$tempo - 1, 0)
 out7<-gls(chumbo ~ tempo + tempo1 + tempo:Groupo + tempo1:Groupo,
           correlation=corSymm (form=~1|ID),
           weights = varIdent(form = ~1 | factor(tempo)),data=datalong)
 summary(out7)
 plot(out7)
# 
# A comparacao dos dois modelos pelo AIC e BIC vence o de tempos discretos 
## 
##=================================================================
#     GEE - MODELO MAXIMAL: TEMPO COMO FATOR
#==================================================================
#   
 require(geepack)
 dados1<-datalong[order(datalong$ID),]
 dados1[1:15,]
#
#  Modelo 1 - Sem Intercepto
# 
# Minimos Quadrados Ordinarios (Parametros: 8 + 1)
 out<-lm(chumbo ~ factor(tempo)+ factor(tempo)*factor(Groupo)-1,data=dados1) 
 summary(out)
# Independencia (Parametros: 8 + 1)
 out1<-geeglm(chumbo ~ factor(tempo)+ factor(tempo)*factor(Groupo)-1, 
 id = ID, data = dados1, corstr = "independence")
 summary(out1)
 # AR(1) # Nao e valido pois nao e igualmente espacado
# 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 residuos
 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) # possivel indicacao da estrutura de simetria composta
 cov(res) # possivel heterocedasticidade
#
#   Modelo 2 (com intercepto)
#
# Minimos Quadrados Ordinarios (Parametros: 8 + 1)
 out<-lm(chumbo ~ factor(tempo)*factor(Groupo),data=dados1) 
 summary(out)
# Independencia (Parametros: 8 + 1)
 out1<-geeglm(chumbo ~ factor(tempo)*factor(Groupo), id = ID, data = dados1, 
 corstr = "independence")
 summary(out1)
 # AR(1) # Nao e valido pois nao e igualmente espacado
# out2<-geeglm(chumbo ~ factor(tempo)*factor(Groupo), id = ID, data = dados1, 
# corstr="corar1")
# summary(out2)
#
 # Simetria Composta (Parametros: 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 residuos
 fit<-out4
 res0<-subset(fit$residuals,dados1$tempo==0)
 res1<-subset(fit$residuals,dados1$tempo==1)
 res4<-subset(fit$residuals,dados1$tempo==4)
 res6<-subset(fit$residuals,dados1$tempo==6)
 res<-data.frame(res0,res1,res4,res6)
 cor(res) # possivel indicacao da estrutura de simetria composta
 cov(res) # possivel heterocedasticidade
#
#
#=================================================================
#     GEE - MODELO: TEMPO continuo (linear e quadratico) e 
#                   termo de Interacao tempo*FATOR
#==================================================================
 library(MASS)
 library(geepack)
 library(gee)
 dados1<-datalong[order(ID),]
 dados1[1:15,]
#
# Centrando na media para evitar multicolinearidade
  dados1[,3]<-dados1[,3]-mean(dados1[,3])
  attach(dados1)
#
# Minimos Quadrados Ordinarios (Parametros: 5 + 1)
  out<-lm(chumbo ~ tempo + I(tempo^2) + tempo:factor(Groupo) +
  I(tempo^2):factor(Groupo))
  summary(out)
#
# Independencia (Parametros: 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 (Parametros: 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 correlacao estimado
## maior que 1
  summary(out4)
#
#=================================================================
#     Grafico 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)
#
#  Independencia 
 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")
# 
#=====================================================================