#
 #=====================================================================
 #                 Dados - Colírio - Crossover
 #======================================================================
 # 
 # O estudo em questão envolve pacientes com hipertensão ocular 
 # que nunca fizeram uso de colírios. O objetivo do estudo é comparar 
 # o efeito de dois colírios disponíveis no mercado: Betoptic e Timoptol.
 # 
 # O estudo foi conduzido de forma longitudinal. As seguintes variáveis 
 # foram medidas ao nível do  paciente: Sexo e Idade (em anos). 
 # E ao nível do episódio foram medidas as seguintes variáveis:
 #	           IOP: Pressão Intra Ocular (mm Hg);
 #
 # O estudo foi conduzido da seguinte forma para cada paciente: 
 # inicialmente a  variável IOP ao nível do episódio é medida, 
 # a seguir um colírio é escolhido de forma aleatória e submetido ao 
 # paciente pelo período de dois meses. Novamente a mesma variável 
 # ao nível do episódio é medida no paciente. Dois meses são 
 # concedidos ao paciente para limpar o efeito do colírio (período de 
 # “washout”) e em seguida é aplicado o outro colírio pelo mesmo período. 
 # Novamente ao final de dois meses a IOP é medida. 
 #  Trinta e dois pacientes foram utilizados no estudo. E a ordem foi
 #   0 - BETOP e a seguir TIMO
 #   1 - TIMO e a seguir BETOP
 #           
 #=======================================================================
 #           Leitura e Preliminares dos Dados 
 #========================================================================
 dados<-read.table("colirio_page.txt",dec=",",h=T) 
 names(dados)
 dim(dados)
 dados[1:15,]
 attach(dados) # o risco é todo seu....
#
#                  Descritiva dos dados
# 
 summary(dados$IOP_sem)
 summary(dados$IOP_betop)
 summary(dados$IOP_timo)
 boxplot(dados$IOP_sem, dados$IOP_betop, dados$IOP_timo,ylab="pressão ocular (mmHg)")
 axis(1, 1:3, c("sem","betop","timol"))
### Cuidado pois o Boxplot não considera o pareamento ####
#       Correlações
 plot(dados$IOP_betop, dados$IOP_timo,xlab="betop (mmHG)",ylab="timo (mmHg)")
 abline(0,1)
 cor(dados$IOP_betop, dados$IOP_timo)
#
#====================================================================
#       Suposição: Não Existe Efeito da Ordem
#       Teste t - Pareado
#====================================================================
 dif1<-dados$IOP_sem - dados$IOP_betop  
 dif2<-dados$IOP_sem - dados$IOP_timo
 dif3<-dados$IOP_betop - dados$IOP_timo
 dif<-dif3
 summary(dif)
 boxplot(dif)
# verificacao de normalidade           
  shapiro.test(dif)
  plot(qqnorm(dif),ylab="quantis")
# teste t para diferencas           
  t.test(dif)
#help.search("wilcoxon")
  wilcox.test(dif)
#=====================================================================
#   Mudando de Banco Largo pra Banco Longo - Ignorando a ordem
#=====================================================================
 resp<-data.frame(dados$IOP_sem,dados$IOP_betop,dados$IOP_timo) 
 resp<-stack(resp) #empilhando o vetor, R gera duas colunas
 IOP<-resp$values
 SEXO<-as.numeric(dados$SEXO)
 sexo<-stack(data.frame(SEXO,SEXO,SEXO))
 idade<-stack(data.frame(dados$IDADE,dados$IDADE,dados$IDADE))
 ordem<-stack(data.frame(dados$Ordem,dados$Ordem,dados$Ordem))
 ident<-rep(1:32,3)
 sexo<-sexo$values
 idade<-idade$values
 ordem<-ordem$values
 ind1 = NULL
 for (i in 1:length(resp$ind)){
 if (resp$ind[i]=="dados.IOP_timo"){ind1[i]=3}
 else if (resp$ind[i]=="dados.IOP_betop"){ind1[i]=2}
 else {ind1[i]=1}
}
 dados1<-data.frame(ident,sexo,idade,IOP,ind1)
 dados1[1:35,]
 boxplot(dados1$IOP~dados1$ind1)
 attach(dados1)
#
#====================================================================
#         Manipulando Banco de Dados Longitudinais
#            Banco Longo (uma medida por linha)
#=====================================================================
#  
 length(unique(dados1$ident)) # qtos individuos no estudo 
 as.data.frame(table(dados1$ident)) #mostra qtas vezes cada individuo aparece
 max(table(dados1$ident)) #numero maximo de medidas por individuo
 min(table(dados1$ident)) #numero minimo de medidas por individuo
#
#======================================================================
#          Estatísticas Descritivas
#=======================================================================
#
 grupo<-dados1$ind1
 tapply(IOP,list(grupo), mean) # media por grupo e tempo
 tapply(IOP,list(grupo), sd) # dp por grupo e tempo
#
#==================================================================
#          Graficos de Perfis
#==================================================================
#
  values<-dados1$IOP
# function(choice,xl=range(time),yl=range(y))
# {# Plots profiles for subjects in vector choice. 
   m=length(unique(dados1$ident))
   n=length(values)/m
   xl=range(grupo)
   yl=range(values)  
   choice=dados1$ident  
   j=ident==choice[1]
   plot(grupo[j],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 2:length(choice)){
   j=ident==choice[it]
   lines(grupo[j],values[j])}
   axis(1, 1:3, c("sem","betop","timol"))
   }
   cat("\n")
#       }
#Fazendo um alisamento no grafico de perfis 
#  plot(lowess(grupo,values),type="l",ylab="Altura",xlab="Idade")  
   lines(lowess(grupo,values),col = "red",cex=2,lwd=2.5)  
#
#==========================================================================
#      Testando o Efeito da Ordem
#============================================================================
#
  tapply(values,list(grupo,ordem), mean)
  tapply(values,list(grupo,ordem), median)
  ordem0<-subset(dados1,ordem==0)
  summary(ordem0$values)
  ordem1<-subset(dados1,ordem==1)
  summary(ordem1$values)
#  dados2<-data.frame(values,grupo)
#  med0<-subset(dados2,ordem==0)
  plot(lowess(ordem0$ind1,ordem0$values),col = "red",cex=1.2,
  xaxt = "n", type="l",xlim=xl,ylab="Pressão",
    xlab="Momento",main="Tratamento - Pressão", ylim=c(10,24))
#  med1<-subset(dados2,ordem==1)
  lines(lowess(ordem1$ind1,ordem1$values),type="l",col="blue",cex=1.2)
  legend(13,140,lty=c(1,2),c("betop-timol","timol-betop"),bty="n")
  axis(1, 1:3, c("sem","betop","timol"))
#
  boxplot(ordem0$values[ordem0$ind1==1],ordem1$values[ordem1$ind1==1],
  ordem0$values[ordem0$ind1==2],ordem1$values[ordem1$ind1==2],
  ordem0$values[ordem0$ind1==3],ordem1$values[ordem1$ind1==3],xaxt = "n")
  axis(1, 1:6, c("sem_0","sem_1","betop_0","betop_1","timol_0","timol_1"))
#
# Teste t para ordens
#
  dif0<-ordem0$values[ordem0$ind1==2]- ordem0$values[ordem0$ind1==3]  
  dif1<-ordem1$values[ordem1$ind1==2]- ordem1$values[ordem1$ind1==3]  
  boxplot(dif0,dif1)
  summary(dif0)
  summary(dif1)
  sd(dif0)
  sd(dif1)
  t.test(dif0,dif1)
#
#######################################################################
###             Material do Andrei ###
##########################################################################
#
#======================================================================================
# Verificando a Hipótese de não influência da Ordem e de Aleatoriedade na Escolha
#======================================================================================
#
# Primeiro teste para influência da ordem
#
dados_aux<-as.data.frame(dados)
#
dados0<-as.data.frame(dados_aux[dados_aux$Ordem==0,])
dados1<-as.data.frame(dados_aux[dados_aux$Ordem==1,])
#
dif0_betop<-dados0$IOP_betop-dados0$IOP_sem
dif0_timo<-dados0$IOP_timo-dados0$IOP_sem
#
dif1_betop<-dados1$IOP_betop-dados1$IOP_sem
dif1_timo<-dados1$IOP_timo-dados1$IOP_sem
#
c(mean(dif0_betop<-dados0$IOP_betop-dados0$IOP_sem),
mean(dif1_betop<-dados1$IOP_betop-dados1$IOP_sem),
mean(dif0_timo<-dados0$IOP_timo-dados0$IOP_sem),
mean(dif1_timo<-dados1$IOP_timo-dados1$IOP_sem))
#
t.test(x=dif1_betop,y=dif0_betop,alternative="two.sided",mu=0,paired=FALSE) #não pareado
t.test(x=dif1_timo,y=dif0_timo,alternative="two.sided",mu=0,paired=FALSE) #não pareado
#
# Segunda teste para influência da ordem
#
primeiro<-c(dados0$IOP_betop,dados1$IOP_timo)
segundo<-c(dados0$IOP_timo,dados1$IOP_betop)
sem<-c(dados0$IOP_sem,dados1$IOP_sem)
#
t.test(x=primeiro,y=sem,alternative="less",mu=0,paired=TRUE) #pareado primeiro
t.test(x=segundo,y=sem,alternative="less",mu=0,paired=TRUE) #pareado segundo
#
t.test(x=primeiro,y=segundo,alternative="two.sided",mu=0,paired=TRUE) #pareado diferenças (two sided)
t.test(x=primeiro,y=segundo,alternative="greater",mu=0,paired=TRUE) #pareado diferenças (greater)
#
# Teste para os medicamentos
#
t.test(x=IOP_betop,y=IOP_sem,alternative="less",mu=0,paired=TRUE) #pareado (Betop é eficaz?)
t.test(x=IOP_timo,y=IOP_sem,alternative="less",mu=0,paired=TRUE) #pareado (Timo é eficaz?)
t.test(x=IOP_betop,y=IOP_timo,alternative="two.sided",mu=0,paired=TRUE) #pareado (existe diferença entre eles?)
#
# Teste para Escolha Aleatória
#
plot(x=c(1,2,3),y=c(mean(dados1$IOP_sem),mean(dados1$IOP_betop),mean(dados1$IOP_timo)),type="l",col="red",
ylim=c(0,22),xlab="Sem, Betop e Timo",ylab="pressão ocular (mmHg)")
lines(x=c(1,2,3),y=c(mean(dados0$IOP_sem),mean(dados0$IOP_betop),mean(dados0$IOP_timo)),col="black")
#
valor<-mean(dados1$IOP_sem)-mean(dados0$IOP_sem)
#
vector<-c()
npath<-5000
for (i in 1:npath) {
	k0<-1
	k1<-1
	tipo0<-c()
	tipo1<-c()
	for (j in 1:32) {
		aux<-rbinom(1,1,0.5)
		if (aux==0) {
			tipo0[k0]<-sem[j]
			k0<-k0+1
		}
		if (aux==1) {
			tipo1[k1]<-sem[j]
			k1<-k1+1
		}
	}
	aux2<-mean(tipo1)-mean(tipo0)
	vector[i]<-aux2
} # simulação de monte carlo com 5.000 paths
#
mean(vector)
#
hist(vector)
#
y<-ecdf(vector)
#
plot(y)
#
2*pnorm(q=valor,mean=mean(vector),sd=sd(vector),lower.tail=FALSE)
#
#  Valor empirico
#
(sum(ifelse(vector>valor,1,0))+sum(ifelse(vector<(-valor),1,0)))/npath
#
#    ####### Fim Material Andrei ###########################################
#
#==========================================================================
#      Modelando a Matriz de Covariância
#============================================================================
#
           # Modelos Marginais sao ajustados utilizando a biblioteca gee
           #
           library(gee)
           dados2<-dados1[order(ident),]
           attach(dados2)
           # Nao Estruturada (quatro parametros)
           out<-gee(values ~ grupo, id = ident, data = dados2, corstr = "unstructured")
           summary(out)
           # Autoregressiva (quatro parametros)
           out<-gee(values ~ grupo, id = ident, data = dados2, corstr = "AR-M",Mv=1)
           summary(out)
           # Esferica (quatro Parametros)
           out<-gee(values ~ grupo, id = ident, data = dados2, corstr = "exchangeable")
           summary(out)
           # Independencia (tres parametro)
           out<-gee(values ~ grupo, id = ident, data = dados2, corstr = "independence")
           summary(out)
           #
           #============================================================================