#
#======================================================
#                 Geração dos Dados
#=======================================================
# Acompanhamento do comprimento de 120 meninas nascidas em
# uma certa maternidade nos primeiros seis meses de vida
#
# Medidas realizadas em quatro instantes: 0,2,4 e 6 meses
#
# Objetivo: descrever o comportamento do comprimento ao longo
# do tempo
#
# Modelos: 
# Resposta: normal
# Média: comprimento médio (idade) = comprimento (idade -1) + 
#                                    3 * idade - 0.5 * idade^2
# Comprimento médio (0) = 50 cm
# Variância (t) = 100 cm^2, para todo t
#
  x0 <- round(rnorm(120, mean=50, sd=10)) # gera 120 normais
# x2
  x2<-c()
  for(it in 1:120){
  x2[it]<-round(rnorm(1,x0[it]+5,10))
  }
# x4 
  x4=range(120)
  for(it in 1:120){
  x4[it]<-round(rnorm(1,x2[it]+5,10))
  }
# x6 
  x6=range(120)
  for(it in 1:120){
  x6[it]<-round(rnorm(1,x4[it]+5,10))
  }
#
  dados<-data.frame(x0,x2,x4,x6) # banco de dados
  attach(dados)
  names(dados)
  dim(dados)
  dados[1:15,]
  write.table(dados, file="simula4.txt", sep=" ", row.names=F)
#           
#=========================================================
#            Leitura e Preliminares dos Dados 
#=========================================================
  dados<-read.table("altura.txt",header=T)
  attach(dados)
  names(dados)
  dim(dados)
  dados[1:15,]
#
# Descritiva dos dados
  summary(x0) 
  summary(x2)
  summary(x4)
  summary(x6)   
  boxplot(x0,x2,x4,x6,ylab="comprimento (cm)")
### Cuidado pois o Boxplot não considera a estrutura longitudinal ####
  cor(dados) # matriz de correlação dos dados
  cov(dados) # matriz de covariancia dos dados
#
#========================================================
#         Mudando de Banco Largo pra Banco Longo
#========================================================
  dados1<-stack(dados) #empilhando o vetor, R gera duas colunas
  attach(dados1)
  ident<-rep(1:120,4)
  grupo<-c(rep(0,120),rep(2,120),rep(4,120),rep(6,120))
  dados1<-data.frame(ident,values,grupo)
  dados1[1:15,]
#
#========================================================
#         Manipulando Banco de Dados Longitudinais
#            Banco Longo (uma medida por linha)
#=========================================================
#  
  length(unique(ident)) # qtos individuos no estudo 
  as.data.frame(table(ident)) #mostra qtas vezes cada individuo aparece
  max(table(ident)) #numero maximo de medidas por individuo
  min(table(ident)) #numero minimo de medidas por individuo
#
#==========================================================
#          Estatísticas Descritivas
#==========================================================
#
  tapply(values,list(grupo), mean)  # media por grupo e tempo
  tapply(values,list(grupo), median)  # media por grupo e tempo
  tapply(values,list(grupo), sd) # dp por grupo e tempo
#
#==================================================================
#          Graficos de Perfis
#==================================================================
#
# function(choice,xl=range(time),yl=range(y))
# {# Plots profiles for subjects in vector choice. 
  m=length(unique(ident))
  n=length(values)/m
  xl=range(grupo)
  yl=range(values)  
  choice=ident  
  j=ident==choice[1]
  plot(grupo[j],values[j],type="l",ylim=yl,xlim=xl,ylab="Comprimento (cm)",xlab="Mês",main="Comp. Récem-nascidos")
  if(length(choice)>1){
  for(it in 2:length(choice)){
  j=ident==choice[it]
  lines(grupo[j],values[j])}
#axis(1, 0:6, c(0,2,4,6))
  }
  cat("\n")
#       }
#Fazendo um alisamento no grafico de perfis 
  lines(lowess(grupo,values),col="2",type="l",lwd=3)  
#
#====================================================================
#  Teste Não-Paramétrico de Friedman - Comparação de Grupos em Blocos
#====================================================================
# Objetivo: comparar grupos (4 instantes - meses)
# Bloco: cada criança
  help.search("friedman")
  friedman.test(values, grupo, ident)
#
#====================================================================
#  Teste de Esfericidade - ANOVA com Correção de GL
#====================================================================
#
# Teste para Esfericidade: mauchly.test(stats)
# Veja também: anova.mlm 
# http//gribblelab.org/2009/03/09/repeated-measures-anova-using-r/
# 
# ANOVA usual (somente válida sob esfericidade)
  anova<-aov(values~factor(grupo)+factor(ident))           
  summary(anova)
#
# Teste de Esfericidade e MANOVA
#
# O primeiro passo é organizar os dados na forma de matriz em que
# cada indivíduo é uma linha
  cvm<- with(dados1, cbind(values[grupo==0],values[grupo==2],
  values[grupo==4],values[grupo==6]))
# O segundo passo é definir um modelo linear multivariado
#
  mlm<-lm(cvm~1)
  mlm
# O terceiro passo é definir a variável de desenho
  rfactor<-factor(c(1,2,3,4))
  rfactor
# Agora carregar a biblioteca car
  library(car)
  mlm1<-Anova(mlm, idata = data.frame(rfactor), idesign=~rfactor, 
  type="III")
#summary(mlm1)
  summary(mlm1, multivariate=F) # exclui os testes multivariados
#
# Pos-hoc Comparações Múltiplas
#
# não existe comparações múltiplas disponível na Anova
# pode-se fazer comparações dois-a-dois usando o teste-t pra 
# diferenças tomando o devido cuidado com o nível de significância
#
           contraste1 <- cvm[,2]-cvm[,1]
           contraste2 <- cvm[,3]-cvm[,2]
           contraste3 <- cvm[,4]-cvm[,3]
           t.test(contraste1)
           t.test(contraste2)
           t.test(contraste3)
           #
           #===============================================================
           #       Explorando os Perfis Médios por Grupo
           #===============================================================
           # Box-plots 
           boxplot(values ~ grupo, notch = T)
           abline(v = 4.5, lty = 2)
           abline(v = 10.5, lty = 2)
           title("Comportamento da Altura no Tempo")
           #
           #====================================================================
           #      Explorando a Estrutura de Dependência dos Tempos
           #====================================================================
           # 
           source("longit.r") #carrega a biblioteca longit.r
           require(car)
           attach(dados)
           ID<-rep(1:120,1)
           dados.std <- perfil.std(y = dados[1:4],id=ID, flag = FALSE)
           scatterplot.matrix(dados.std[, 2:5], diagonal = "histogram", 
           smooth = F)
           #graphics.off()
           cor(dados.std[, 2:5])
           #
           #==========================================================================
           #      Modelando a Matriz de Covariância - GEE
           #============================================================================
           #
           # Modelos Marginais sao ajustados utilizando a biblioteca gee
           #
           library(gee)
           dados2<-dados1[order(ident),]
           attach(dados2)
           # Minimos Quadrados Ordinarios 
           x1<-rep(c(0,1,0,0),120)
           x2<-rep(c(0,0,1,0),120)
           x3<-rep(c(0,0,0,1),120)
           out<-lm(values ~ x1+x2+x3)
           out 
           # Nao Estruturada (dez parametros)
           out<-gee(values ~ factor(grupo), id = ident, data = dados2, corstr = "unstructured")
           summary(out)
           # Autoregressiva (quatro parametros)
           out<-gee(values ~ factor(grupo), id = ident, data = dados2, corstr = "AR-M",Mv=1)
           summary(out)
           # Esferica (quatro Parametros)
           out<-gee(values ~ factor(grupo), id = ident, data = dados2, corstr = "exchangeable")
           summary(out)
           # Independencia (tres parametro)
           out<-gee(values ~ factor(grupo), id = ident, data = dados2, corstr = "independence")
           summary(out)
           #
           #==========================================================================
           #      Modelando a Matriz de Covariância - GEEPACK
           #============================================================================
           #
           # Modelos Marginais sao ajustados utilizando a biblioteca geepack
           #
           library(geepack)
           dados2<-dados1[order(ident),]
           attach(dados2)
           # Nao Estruturada (dez parametros)
           out<-geeglm(values ~ factor(grupo), id = ident, data = dados2, corstr = "unstructured")
           summary(out)
           # Autoregressiva (quatro parametros)
           out<-geeglm(values ~ factor(grupo), id = ident, data = dados2, corstr = "AR1")
           summary(out)
           # Esferica (quatro Parametros)
           out<-geeglm(values ~ factor(grupo), id = ident, data = dados2, corstr = "exchangeable")
           summary(out)
           # Independencia (tres parametro)
           out<-geeglm(values ~ factor(grupo), id = ident, data = dados2, corstr = "independence")
           summary(out)
           #
           #============================================================================
           #      Modelando a Matriz de Covariância - GLS (library -nlme)
           #==========================================================================
           #
           # Modelos Marginais ajustados pela biblioteca nlme
           #
           library(nlme)
           attach(dados1)
           # Não Estruturada (10 componentes de variância)
           out<-gls(values ~ factor(grupo), dados1, correlation = corSymm(form = ~1 |ident),
           weights = varIdent(form = ~1| grupo))
           summary(out)
           # Simetria Composta (5 componentes de variância)
           out1<-gls(values~factor(grupo),correlation=corCompSymm(form=~1|ident),
           weights = varIdent(form = ~1| grupo))
           out1
           summary(out1)
           # AR(1) (5 componentes de Variância)
           out2<-gls(values~factor(grupo),correlation=corAR1(form=~1|ident),
           weights = varIdent(form = ~1| grupo))
           out2
           summary(out2)
           anova(out,out1,out2)
           # Não Estruturada - Homocedastica (7 componentes de variância)
           intervals(out) # IC para os componentes de variância
           out3<-gls(values ~ factor(grupo), dados1, 
           correlation = corSymm(form = ~1 |ident))
           out3
           anova(out,out3)
           # Não Estruturada (10 componentes de variância - MV)
           out4<-gls(values ~ factor(grupo), dados1, correlation = corSymm(form = ~1 |ident),
           method="ML", weights = varIdent(form = ~1| grupo))
           summary(out4)
           #=================================================================================