#############################################################################################
### Modelos que utilizan variables aleatorias para modelar la falta de independencia      ###
### en las observaciones. Modelos mixtos: Modelos cruzados, lmer y valores de prob.       ###
### Roberto Munguia-Steyer                                   2011/11/15 y 2011/11/17      ###
#############################################################################################

R

##  Cargando paquetes

library(lme4)
library(languageR)
library(lattice)
library(faraway)
library(languageR)
library(MEMSS)
library(mlmRev)
library(SASmixed)
library(effects)
library(agricolae)


##  Modelos mixtos y efectos cruzados : Cuando existe mas de una fuente de variacion que queremos controlar 
##  el disenho RCB es insuficiente, en este caso el disenho adecuado seria el de cuadrados latinos.
##  Un disenho de cuadrados latinos es un caso particular de los modelos mixtos con efectos cruzados


# Produccion de leche en bovinos, bloques vaca y tiempo, Oehlert (2000) ex 13.12 pp 341.
# Columns are cow, period, treatment (1--roughage, 2--limited grain, 3--full
# grain), and milk production in pounds per six week period.

bd <- read.table("pleche.txt", header=T, sep="\t")

names(bd) <- c("vaca", "periodo", "trat", "leche")
bd$vaca <- factor(bd$vaca)
bd$periodo <- factor(bd$periodo)
bd$trat <- factor(bd$trat)
levels(bd$trat) <- c("fibra", "gramlim", "gramtot")

xtabs(~vaca+periodo, data=bd)

#Por que en este periodo es una variable aleatoria?
m1 <- lmer(leche ~ trat + (1|vaca) + (1|periodo), data=bd)

summary(m1)
anova(m1)

#Validacion del modelo
plot(fitted(m1), resid(m1))
qqnorm(resid(m1))


#Calculo del valor de probabilidad usando MCMC y bootstraps parametricos

#Bootstrap parametrico, modelo reducido y global
#Numero de simulaciones

#Modelo nulo
mn <- lmer(leche~ 1 + (1|vaca) + (1|periodo), data=bd)

pboot <- function(mn,m1) {
s <- simulate(mn)
L0 <- logLik(refit(mn,s))
L1 <- logLik(refit(m1,s))
2*(L1-L0)
}
obsdev2 <- c(2*(logLik(mn)-logLik(m1)))
psn <- replicate(1000, pboot(mn,m1))

#Prueba de razon de verosimilitud, asumiendo que el estadistico tiene una dist de ji cuadrada
anova(mn,m1)

#Valor de probabilidad, usando una prueba de ji cuadrada
pchisq(abs(obsdev2), 1, lower.tail=FALSE)

#Valor obtenido mediante el bootstrap parametrico, en este caso la correccion resulta trivial
mean(psn>abs(obsdev2))

#Cadenas de markov montecarlo (MCMC)
mcmc = pvals.fnc(m1, nsim=10000, withMCMC=TRUE)
#aovlmer.fnc(m2b, mcmc$mcmc, c("drogacp", "drogacs")), categorico
str(mcmc)

#Cuales son los niveles del tratamiento que difieren entre si
mcmc$fixed
aovlmer.fnc(m1, mcmc$mcmc, c("tratgramlim", "tratgramtot"))

efectos <- allEffects(m1, xlevels=list(trat=levels(trat)))
efectos
summary(efectos)


# b) Data abrasion, perdida de peso por uso, cuatro materiales diferentes (A,B,C,D), 
#la maquina procesa cuatro muestras a la vez, efecto de posicion
help(abrasion)
str(abrasion)

xtabs(~run+position, data=abrasion)

names(abrasion) <- c("turno", "posicion", "material", "desgaste")

m2 <- lmer(desgaste ~ material + (1|posicion) + (1|turno), abrasion)

#Validacion del modelo
plot(fitted(m2), resid(m2))
qqnorm(resid(m2))

#Efecto del material
anova(m2)

summary(m2)

##Modelo nulo
mn2 <- lmer(desgaste~ 1 + (1|posicion) + (1|turno), data=abrasion)

#Bootstrap parametrico
pboot <- function(mn2,m2) {
s <- simulate(mn2)
L0 <- logLik(refit(mn2,s))
L1 <- logLik(refit(m2,s))
2*(L1-L0)
}
obsdev4 <- c(2*(logLik(mn2)-logLik(m2)))
psn2 <- replicate(1000, pboot(mn2,m2))

#Prueba de razon de verosimilitud, asumiendo que el estadistico tiene una dist de ji cuadrada
anova(mn2,m2)

#Valor de probabilidad, usando una prueba de ji cuadrada
pchisq(abs(obsdev4), 3, lower.tail=FALSE)

#Valor obtenido mediante el bootstrap parametrico, en este caso la correccion resulta trivial por el bajo valor de p
mean(psn2 > abs(obsdev4))

#Comparaciones entre niveles del factor
efectos2 <- allEffects(m2, xlevels=list(material=levels(material)))
efectos2
summary(efectos2)

# Ejercico 10.1, base de datos melon del paquete agricolae, cuadrados latinos
# Realiza un modelo  mixto que cuente como  variables aleatorias la hilera y la columna en la que se encuentraba cada melon, 
# considerando la variedad del frijol como la variable predictora fija y el rendimiento como la variable de respuesta
# ¿Como sabemos que en realidad es un cuadrado latino? pista: usar la funcion xtabs
# que te dicen las desviaciones estandar de las variables aleatorias, realiza seleccion de modelos usando el AIC
# Valida el modelo y determina si existen diferencias el rendimiento de acuerdo a la variedad del melon y si es el caso
# indica los valores predichos para cada variedad y cuales variedades difieren entre si