#
#==============================================================
#       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
#  Sigma=matrix(c(225,112.5,112.5,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 
#=========================================================
  dados<-read.table("testet.txt",header=T)
  names(dados)    # formato largo
  dim(dados)
  dados[1:15,]
  plot(dados$x1,dados$x2,xlab="y1",ylab="y2")
  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)
  t.test(dados$x1-dados$x2)
#   Usando o verdadeiro valor da Var(mean(x1)-mean(x2))
 numt<-mean(dados$x1-dados$x2)
 dentsq<-(1/n)*(cov(dados)[1,1]+cov(dados)[2,2]-2*cov(dados)[1,2])
 t<-numt/sqrt(dentsq)
 t
#      Ignorando covari?ncia
 dentsq_er<-(1/n)*(cov(dados)[1,1]+cov(dados)[2,2])
 t_er<-numt/sqrt(dentsq_er)
 t_er
#==================================================================
#  Usando modelo de Regressão
#
#================================================================
#         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)
  dim(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")
# 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 Covariancia - 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)
#
# Simetria composta e mesma vari?ncia (2 componentes de vari?ncia)
#   estimando os componentes de variância  
  out1<-gls(values ~ grupo, dados1, correlation = corSymm(form = ~1 |ident))
  summary(out1)
  intervals(out1) # IC para os componentes de vari?ncia
#
# Simetria composta e mesma vari?ncia (2 componentes de vari?ncia)
#   fixando a correlação = 0.8  
  cs <- corCompSymm(0.8, form = ~ 1 | ident,fixed=T)
  out2<-gls(values ~ grupo, dados1, correlation = cs)
  summary(out2)
  intervals(out2) # IC para os componentes de vari?ncia
#  
#==========================================================================
#       Extensao 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(nlme)
  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(dados3)
#
#========================================================
#         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,]
  dim(dados4)
  names(dados4)
#
#  Descritiva dos Dados
  require(ggplot2)
  names(dados4)[3] <- "Momento" 
  names(dados4)[2] <- "Press?o" 
  p<-ggplot(dados4,aes(Momento, Press?o,group=ident)) + geom_line()
  p 
# 
  names(dados4)[2] <- "values" 
  interaction.plot(dados4$Momento,rep(1,300),dados4$values,xlab="momento",
  ylab="press?o (mmHg)",legend=F, main="grafico de medias")
# ANOVA usual: um fator em blocos (somente v?lida sob esfericidade)
  anova<-aov(values~factor(grupo)+factor(ident),data=dados4)           
  summary(anova)
#
  anova_random.lme<-lme(values ~ factor(grupo), data=dados4, random= ~1|ident, m="REML")
  summary(anova_random.lme) 
#  
# 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
#
 cvm<-cbind(x1,x2,x3)
 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 de Esfericidade - ANOVA com Corre??o de GL
#====================================================================
#
# Teste para Esfericidade: mauchly.test(stats)
# Veja tamb?m: anova.mlm 
# http://www.psych.upenn.edu/~baron/rpsych/rpsych.html#htoc60
# 
# 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
#
#===============================================================================
#  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$Momento, dados4$ident)
#================================================================================