#############################################################################################
### Modelos que utilizan variables aleatorias para modelar la falta de independencia      ###
### en las observaciones. Modelos jerarquicos, multinivel, anidados, mixtos.              ###
### Roberto Munguia-Steyer                                   2011/11/08 y 2011/11/10      ###
#############################################################################################

R

##  Cargando paquetes
#install.packages(c("MEMSS", "RLRsim", "glmmML","languageR","mlmRev","SASmixed", "foreign"))
library(nlme)
#library(MEMSS)
library(lattice)
library(ggplot2)
library(bbmle)
library(effects)
library(multcomp)
library(faraway)
library(foreign)

##  a) El modelo más sencillo que lidia con la falta de independencia entre observaciones.
data(sleep)
#help(sleep)

str(sleep)
head(sleep)


sn <- sleep
names(sn) <- c("suext", "droga", "ind")

#graficando
xyplot(suext~ droga|ind, data=sn, type="o", ylab= "Horas extras de sueño", xlab= "Droga inductora de sueño")

#Grafica mas elegante
#con ggplot
sn$drogac <- as.numeric(sn$droga)
sugraf <- ggplot(sn, aes(drogac,suext))+ scale_y_continuous("Horas extras de sueño") + scale_x_continuous("Droga")
sugraf <- sugraf + geom_point() + geom_line() + facet_wrap(~ind, ncol = 5)

#Analisis de prueba de t pareada

m0 <- t.test(suext~droga, paired=T, data=sn)
m0

#Modelo mixto, individuo como variable aleatoria 

m1 <- lme(suext~droga, random=~1|ind, data=sn)
anova(m1)
summary(m1) #La desviación estandar de los efectos aleatorios indica la variacion existente entre los distintos individuos.

#Validacion del modelo
plot(m1)

er <- resid(m1, type="normalized")
boxplot(er~droga, data=sn)

#Los modelos mixtos con la funcion lme, permiten incorporar estructuras de correlacion y varianza
#probemos modelar varianzas diferentes para cada droga

m2 <- update(m1, weights=varIdent(form=~1|droga))

logLik(m1); logLik(m2)
AICtab(m1,m2, base=T, weights=TRUE) #El poder de explicacion ganado por anhadir un parametro adicional es muy pobre, quedemonos con el modelo mas parsimonioso

##  b) Modelos mixtos con un factor de mas de dos niveles, ejemplo tomado de Oehlert A First Course in
##  Design and Analysis of Experiments, Analisis de varianza con bloques


mealybug <- data.frame(trt = gl(3, 1, 15, labels=c("water","spores","oil")),
                       block = gl(5,3),
                       avediff = c(-7.5, 1.5,  7.5, 11.5, 19.5, 32.5,
                                    9.5, 1.5, 15,    4.5,  2,   16,
                                    3.5, 5,   11))
#Observemos la base de datos
mealybug

cc <- mealybug
names(cc) <- c("trat", "bloque", "difins")

#Inspeccion visual
p <- ggplot(cc, aes(trat, difins, group = bloque, shape = trat))
p+ geom_point()+facet_wrap(~bloque, ncol=5)

#Ajuste del modelo, planta es la variable aleatoria, nuestro bloque

m3 <- lme(difins~trat, random=~1|bloque, data=cc)
summary(m3)
anova(m3)

#Validacion del modelo, con pocas observaciones se complica la inspeccion
plot(m3)
qqnorm(m3)

#¿El tratamiento con esporas difiere del que utiliza un aceite?
cm <-glht(m3, linfct = mcp(trat= "Tukey"))
iccm<- confint(cm)

par(mar=c(5,8,3,2), las=1)
plot(iccm)

#valores predichos y sus intervalos de confianza
efcc<- allEffects(m3, xlevels=list(trat=levels(cc$trat)))
efcc

#Valores predichos con los intervalos de confianza al 95 %
summary(efcc)

# Ejercicio 8.1
# en un campo de cultivo utilizaron cinco metodos de irrigacion en 8 diferentes bloques(los cuales podian diferir en tipo de suelo y luminosidad)
# Utiliza un modelo mixto para determinar si existen diferencias en el peso de los frutos en función del tipo de irrigacion
# Los datos se encuentran en la base de datos ("METHOD.sav"), la cual se encuentra en formato del programa SPSS
# Utiliza el paquete foreign para importar la base de datos a R con los siguientes comandos 
bds <- read.spss("METHOD.sav")
bds <- data.frame(bds)

#1. Realiza la inspeccion visual y ajuste del modelo
#2. Valida el modelo
#3. Existe diferencia en el numero de frutos producidos con los distintos tipos de irrigacion
#4. Determina cuales son los grupos que difieren entre si

str(bds)
names(bds) <- c("bloque", "irrig", "pesfrut")
levels(bds$irrig)
table(bds$irrig)

levels(bds$irrig)[5]<- "irrigdesc"


xyplot(pesfrut~irrig|bloque, data=bds)

girri <- ggplot(bds, aes(irrig, pesfrut, group = bloque, shape = irrig))
girri+ geom_point()+facet_wrap(~bloque, ncol=4)

#Ajuste del modelo

m4 <- lme(pesfrut~irrig, random=~1|bloque, data=bds)
anova(m4)
summary(m4)

#Validacion del modelo
e2 <- resid(m4, type="normalized")
boxplot(e2~irrig, data=bds)
qqnorm(m4)
plot(fitted(m4), e2)

m5<- update(m4, weights=varIdent(form=~1|irrig))
summary(m5)


logLik(m4); logLik(m5)
AICtab(m4,m5, base=T, weights=TRUE)

#Que niveles difieren entre si

cm2 <-glht(m4, linfct = mcp(irrig= "Tukey"))
iccm2<- confint(cm2)

summary(cm2, test = univariate())
summary(cm2, test = adjusted("bonferroni"))

par(mar=c(5,8,3,2), las=1)
plot(iccm2)



## c) Los modelos pueden presentar solo la variable aleatoria, anova de efectos aleatorios

data(pulp)
#help(pulp)
names(pulp)
summary(pulp)

#Variabilidad en el grado de eficiencia de los trabajadores, que tan blancas les salen las hojas de papel
stripchart(bright~operator,data=pulp,method="jitter",vertical=TRUE, pch=19, ylab= "Nivel de blancura", xlab="Trabajador")

am <- lme(bright ~ 1, data = pulp, random = ~ 1 | operator)
summary(am) #Mayor variacion dentro que entre los operadores

#Validacion del modelo
par(las=1, mfrow=c(1,2))
qqnorm(resid(am),main="")
plot(fitted(am),resid(am),xlab="Predichos",ylab="Residuos")

## d) Analisis de varianza anidado, base de datos de Crawley (2007), pag. 475
## Glicogeno en el higado de ratas

rat <- read.table("rats.txt", sep="\t", header=T)
str(rat)
names(rat)


names(rat) <- c("glico", "trat", "rata", "higado")
rat$trat <- as.factor(rat$trat)
rat$higado <- as.factor(rat$higado)
rat$rata <- as.factor(rat$rata)

#Son seis ratas, corregir los codigos
rat$rata <- rep(1:6, each=6)
rat


#Inspeccion visual
#con ggplot
grat <- ggplot(rat, aes(trat, glico, shape = higado))
grat2 <- grat  + geom_point() + facet_wrap(~rata, ncol=3)  + scale_y_continuous("Glicogeno")
grat2

#Ajuste del modelo

m4 <- lme(glico ~ trat, random=~1|rata/higado, data=rat)
summary(m4)
anova(m4)

#Analisis de componentes de la varianza, en que nivel se encuentra la variacion
VarCorr(m4)
vars<- c(36.0843,14.1617,21.1678)
100*vars/sum(vars)

#Varianza encontrada entre ratas 50.52 %, 19.83 % dentro de los pedazos de higado y 29 % en las diferentes lecturas dentro de los pedazos

#Analisis de residuos
par(mfrow=c(1,2))
qqnorm(resid(m4))
plot(fitted(m4), resid(m4))

#Con el paquete lme4
sessionInfo()

detach("package:effects")
detach("package:nlme")

sessionInfo()

#Cargando el paquete lme4
library(lme4)
rat$rathig <- paste(rat$rata, rat$higado, sep="") #Niveles unicos de higado dentro de rata

m4b <- lmer(glico ~ trat + (1|rata) + (1|rathig), data=rat)
summary(m4b)
anova(m4b) #donde estan los valores de probabilidad
#ver http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-are-p_002dvalues-not-displayed-when-using-lmer_0028_0029_003f

#Desvinculando lme4 y cargando nlme

detach("package:lme4")
library(nlme)

## e) Modelos mixtos con intercepto y pendiente aleatorias, modelando estructuras de correlacion y varianza.
data(BodyWeight)
#help(BodyWeight)

pc <- BodyWeight

str(pc)
head(pc)
names(pc) <- c("peso","tiempo","rata","dieta")
ratet <- paste("r", 1:16, sep="")

levels(pc$rata) <-ratet

xyplot(peso~ tiempo|rata, groups=dieta,data=pc, type="o")
 
#con ggplot
p <- ggplot(pc, aes(tiempo, peso, group = dieta, shape = dieta))
pg <- p + geom_line(colour = "darkgrey") + geom_point() + facet_wrap(~dieta+rata, ncol=8)+ scale_y_continuous("Peso") +scale_x_continuous("Tiempo")
pg

#Con intercepto aleatorio
mr1 <- lme(peso~tiempo*dieta, random=~1|rata, data=pc)
anova(mr1)
summary(mr1)


#Adicionalmente con pendiente aleatorias

mr2 <- lme(peso~tiempo*dieta, random=~tiempo|rata,data=pc)

AIC(mr1);AIC(mr2)



#Validacion del modelo
qqnorm(mr2)
plot(mr2)


e2 <-resid(mr2, type="normalized")

#Inspeccion de la varianza por grupos
boxplot(e2~pc$dieta)

#Existe incremento en la varianza en funcion del tiempo transcurrido?
plot(mr2, resid(.)~tiempo, abline=0)# falta de independencia, tambien autocorrelacion temporal? 

xyplot(resid(mr2)~pc$tiempo|pc$rata)# Residuos estandarizados por tiempo e individuo

#Heterogeneidad de varianza por grupos
mr3 <- lme(peso~tiempo*dieta, random=~tiempo|rata, weights=varIdent(form=~1|dieta), data=pc)


#+ Autocorrelacion temporal
plot(ACF(mr2),alpha=0.05)

mr4 <- update(mr3, correlation=corAR1(form=~tiempo))

AICtab(mr1,mr2,mr3,mr4, base=T, weights=T)

anova(mr4)
summary(mr4)

# Podemos observar los contrastes entre las pendientes de los niveles 1-2 y 1-3, y la 2-3
anova(mr4, L = c(dieta2 =1, dieta3 =-1))

# Si quisieramos realizar selección de modelos, necesitamos realizar la optimizacion de los modelos con ML en vez de REML
mr4a <- update(mr4, method="ML")
mr5 <- update(mr4, .~.-dieta:tiempo, method= "ML")

anova(mr4a, mr5) #nos quedamos con el modelo mr4, la interacción es significativa.

#Validacion del modelo
qqnorm(mr4)
plot(mr4) 

#Ejercicio 8.2
#Efectue el equivalente de un analisis de medidas repetidas con la base de datos Orthodont ubicada en el paquete nlme
# 1. Efectue la inspeccion visual, ajuste y validacion del modelo con diferentes estructuras de error: intercepto y 
# pendiente aleatoria
# 2. Realice la seleccion del modelos, empleando el metodo de optimizacion por maxima verosimilitud
# Son diferentes las tasas de crecimiento entre hombres y mujeres