#
#======================================================
#       Amostra pareada: pré e pós tratamento
#                  100 pacientes
#                 Geração dos Dados
#    Research Question: Existe diferença entre o pré e o pós?
#=======================================================
# Grupo 1: pacientes com hipertensao Sistolica N(150,15)
# Pressao Sistolica Normal ate 120mmHG
# Grupo 2: mesmos pacientes pós tratamento N(110,15)
# \rho = 0.8 (coeficiente de correlacao)
# cov=0.8*15*15=180
  require(MASS)
  n<-100
  mu=c(150,110) #vetor de médias pré e pós tratamento
  Sigma=matrix(c(225,180,180,225),2,2) #matriz de covariância
  resp<-round(mvrnorm(n,mu,Sigma))
  x1<-resp[,1]
  x2<-resp[,2]
  dados<-data.frame(x1,x2) # banco de dados
  names(dados)
  dim(dados)
  dados[1:15,]
#
  write.table(dados, file="testet.txt", sep=" ", row.names=F)
#           
#=========================================================
#           Leitura e Preliminares dos Dados 
#=========================================================
  require(epicalc)
  dados<-read.table("testet.txt",header=T)
  names(dados)    # formato largo
  dim(dados)
  dados[1:15,]
  plot(dados$x1,dados$x2)
  cor(dados$x1,dados$x2)
#
# Descritiva dos dados 
  summary(dados$x1)
  summary(dados$x2)
  boxplot(dados$x1,dados$x2,ylab="pressão sistólica (mmHg)")
  axis(1, 1:2, c(1,2))
### Cuidado pois o Boxplot não considera o pareamento ####
  cor(dados)
  cov(dados)
#
#========================================================
#         Mudando de Banco Largo pra Banco Longo
#========================================================
  dados1<-stack(dados) #empilhando o vetor, R gera duas colunas
#  attach(dados1)
  ident<-rep(1:100,2)
  grupo<-c(rep(1,100),rep(2,100))
  dados1<-data.frame(ident,dados1[1],grupo)
  dados1[1:15,]
  names(dados1)
#
#========================================================
#         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(dados1$values,list(dados1$grupo), mean) # media por grupo 
  tapply(dados1$values,list(grupo),median) # mediana por grupo 
  tapply(dados1$values,list(grupo), sd) # dp por grupo 
  tapply(dados1$values,list(grupo),summary)
#
#==================================================================
#          Graficos de Perfis
#==================================================================
#
# function(choice,xl=range(time),yl=range(y))
# {# Plots profiles for subjects in vector choice. 
# 
   m=length(unique(ident)) # somente 10 perfis m=10
   n=length(dados1$values)/m
   xl=range(grupo)
   yl=range(dados1$values)  
   choice=ident # somente 10 perfis choice=ident[1:20]
   j=ident==choice[1]
   plot(grupo[j],dados1$values[j],xaxt = "n", type="l",ylim=yl,xlim=xl,ylab="Pressão",
   xlab="Momento",main="Tratamento - Pressão")
   if(length(choice)>1){
   for(it in 3:length(choice)){
   j=ident==choice[it]
   lines(grupo[j],dados1$values[j])}
   axis(1, 1:3, c(1,2,3))
   }
# Fazendo um alisamento no grafico de perfis 
  lines(lowess(dados1$grupo,dados1$values),col="red",lwd=2)  
# Fazendo um gráfico de médias
  mean<-c(mean(dados$x1),mean(dados$x2))
  lines(grupo[j],mean,type="l",col="blue",lwd=2)
#====================================================================
#       Teste t - Pareado
#====================================================================
  dif<-dados$x1-dados$x2  
  boxplot(dif)
# verificacao de normalidade           
  shapiro.test(dif)
  plot(qqnorm(dif),ylab="quantis amostrais",xlab="quantis normal")
  shapiro.qqnorm(dif)
# teste t para diferencas           
  t.test(dif)
  t.test(values~grupo, paired=T,data=dados1)
#help.search("wilcoxon") teste não-paramétrico de Wilcoxon
  wilcox.test(dif)
#
#============================================================================
#       ************MODELO MARGINAL************************** 
#      Modelando a Matriz de Covariância - GLS (library -nlme)
#==========================================================================
  library(nlme)
# Ignorando a correlação entre as medidas antes e depois
  out0<-lm(values ~ grupo, data=dados1)
  summary(out0)
  plot(out0)
#
# Não Estruturada (3 componentes de variância)
  out<-gls(values ~ grupo, dados1, correlation = corSymm(form = ~1 |ident),
  weights = varIdent(form = ~1| grupo))
  summary(out)
  intervals(out) # IC para os componentes de variância
## Obs. o Std. Error muda muito (de 2.07 para 0.951)
# Simetria Composta (3 componentes de variância)= Não Estruturada
  out1<-gls(values~ grupo,dados,correlation=corCompSymm(form=~1|ident),
  weights = varIdent(form = ~1| grupo))
  summary(out1)
# AR(1) (3 componentes de Variância)= Não Estruturada
  out2<-gls(values~ grupo,dados,correlation=corAR1(form=~1|ident),
  weights = varIdent(form = ~1| grupo))
  summary(out2)
# Não existe diferença entre os ajustes pois somente temos dois tempos
# Não Estruturada - Homocedastica (2 componentes de variância)
  out3<-gls(values ~ grupo, dados,correlation = corSymm(form = ~1 |ident))
  summary(out3)
  anova(out,out3)
# Assim como foi gerado, aceitamos a homocedasticidade
#=================================================================================
#       ************MODELO MARGINAL************************** 
#      Modelando a Matriz de Covariância - GEE
#============================================================================
# Ignorando a correlação entre as medidas antes e depois
  out0<-lm(values ~ grupo, data=dados1)
  summary(out0)
#
# Modelos Marginais sao ajustados utilizando a biblioteca gee
#
 library(geepack)
 dados2<-dados1[order(ident),]
# Esferica (quatro Parametros)
  out<-geeglm(values ~ grupo, id = ident, data = dados2, corstr = "exchangeable")
  summary(out)
# Independencia (tres parametro)
  out<-geeglm(values ~ grupo, id = ident, data = dados2, corstr = "independence")
  summary(out)
#
# Obs. que os resultados são iguais pois usa a variância robusta.
# Como o GEE calcula os dois componentes de variância
# Utilizando os resíduos - Somente um passo com os resíduos de MQO
  names(out)
  out$residuals
  x<-rep(c(1,2),100) # separando os resíduos do mesmo indivíduo
  res<-data.frame(x,out$residuals)
  sd(out$residual)^2 # variância dos erros
#
  res1<-subset(res,x==1) # resíduo do antes
  res2<-subset(res,x==2) # resíduo do depois
  cor(res1[,2],res2[,2]) # correlação
#
#==========================================================================
#       Extensão para um terceiro tempo
#                  100 pacientes
#     Research Question: Avaliação temporal da PS
#                 Geração dos Dados
#==========================================================================
# Tempo 0: pacientes com hipertensao Sistolica N(150,15)
# Pressao Sistolica Normal ate 120mmHG
# Tempo 30: mesmos pacientes 30 dias após tratamento N(110,15)
#
# Tempo 60: mesmos pacientes 60 dias após tratamento N(110,15)
# \rho = 0.8 (coeficiente de correlacao)
# cov=0.8*15*15=180
  require(MASS)
  n<-100
  mu=c(150,110,110) #vetor de médias tempos:0,30,60
  Sigma=matrix(c(225,180,180,180,225,180,180,180,225),3,3) #matriz de covariância
  resp<-round(mvrnorm(n,mu,Sigma))
  x1<-resp[,1]
  x2<-resp[,2]
  x3<-resp[,3]
  dados3<-data.frame(x1,x2,x3) # banco de dados
  names(dados3)
  dim(dados3)
  dados3[1:15,]
  cor(dados3)
  cov(dados)
#
#====================================================================
#        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/
# 
# Teste de Esfericidade e MANOVA
#
  library(car)
#
# O primeiro passo é organizar os dados na forma de matriz em que
# cada indivíduo é uma linha
  cvm<-cbind(x1,x2,x3)
# 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))
  rfactor
# Agora carregar a biblioteca car
  mlm1<-Anova(mlm, idata = data.frame(rfactor), idesign=~rfactor, 
  type="III")
#summary(mlm1)
  summary(mlm1, multivariate=F) # exclui os testes multivariados
#
#========================================================
#         Mudando de Banco Largo pra Banco Longo
#========================================================
  values<-stack(dados3) #empilhando o vetor, R gera duas colunas
  ident<-rep(1:100,3)
  grupo<-c(rep(1,100),rep(2,100),rep(3,100))
  dados4<-data.frame(ident,values[1],grupo)
  dados4[1:15,]
#
  interaction.plot(dados4$grupo,rep(1,300),dados4$values)
# ANOVA usual (somente válida sob esfericidade)
  anova<-aov(values~factor(grupo)+factor(ident),data=dados4)           
  summary(anova)
#
# 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[,1]-cvm[,3]
 t.test(contraste1)
 t.test(contraste2)
 t.test(contraste3)
#
#===============================================================================
#  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(dados4$values, dados4$grupo, dados4$ident)
#================================================================================
#  Comparação de Dois Grupos ao longo do tempo
#       Split-plot Design
#   Tutorial: http://statistics.ats.ucla.edu/stat/r/seminars/Repeated_Measures/repeated_measures.htm
#===============================================================================
  demo1 <- read.csv("http://www.ats.ucla.edu/stat/data/demo1.csv")
# Convert variables to factor
  demo1 <- within(demo1, {
  group <- factor(group)
  time <- factor(time)
  id <- factor(id)
  })
  par(cex = .6)
#
  with(demo1, interaction.plot(time, group, pulse,
  ylim = c(5, 20), lty= c(1, 12), lwd = 3,
  ylab = "mean of pulse", xlab = "time", trace.label = "group"))
#
  demo1.aov <- aov(pulse ~ group * time + Error(id), data = demo1)
  summary(demo1.aov)
#
#
 demo4 <- read.csv("http://www.ats.ucla.edu/stat/data/demo4.csv")
## Convert variables to factor
  demo4 <- within(demo4, {
  group <- factor(group)
  time <- factor(time)
  id <- factor(id)
  })
  par(cex = .6)
#
  with(demo4, interaction.plot(time, group, pulse,
  ylim = c(10, 60), lty = c(1, 12), lwd = 3,
  ylab = "média bat. cardíaco", xlab = "tempo", trace.label = "grupo"))
#
  demo4.aov <- aov(pulse ~ group * time + Error(id), data = demo4)
  summary(demo4.aov)
########################################################################