#############################################################################################
### 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
#install.packages(c("languageR", "agricolae", "agridat"))
#library(nlme)
#library(agridat)
library(lme4)
library(ggplot2)
library(languageR)
library(lattice)
library(faraway)
library(languageR)
library(MEMSS)
library(mlmRev)



##Valores de probabilidad, bootstrap parametrico Faraway Extending the linear model pp 180
## Los valores de probabilidad obtenidas a traves de los LRTs no siempre son confiables
## Estadistico de Chi que no siempre sigue la distribucion esperada por esa probabilidad

# a) Produccion de penicilina a partir de un licor de maiz, bastante variable, mezcla solo alcanza para cuatro corridas

data(penicillin)
summary(penicillin)
##m1.lme <- lme(yield ~treat, random=~1|blend, penicillin)

m1 <- lmer(yield~treat + (1|blend), penicillin)
summary(m1)
anova(m1)

#Si el estadistico F realmente se encuentra acorde a la distribucion de F
pf(1.24, df1=4, df2=16, lower.tail=FALSE)


#Ajustamos el modelo global y el restringido, con el metodo de optimizacion de maxima verosimilitud
m0 <- lmer(yield~1 + (1|blend), penicillin, REML=FALSE)
m1a <- lmer(yield~treat + (1|blend), penicillin, REML=FALSE)

#Prueba de razon de verosimilitud, estadistico de chi cuadrada
anova(m0, m1a)

#Valores de probabilidad tambien se pueden obtener asi
pchisq(4.05, df=3, lower.tail=FALSE)

#Razon de verosimilitud, cotejala con el valor obtenida con el LRT previo
obsdev <- c(2*(logLik(m1a)-logLik(m0)))
obsdev

#Bootstrap parametrico
lrstat <- numeric(1000)
for(i in 1:1000){
ryield <- unlist(simulate(m0))
modnsim <- lmer(ryield ~ 1 + (1|blend), penicillin, REML=FALSE)
modsim <- lmer(ryield ~ treat + (1|blend), penicillin,REML=FALSE)
lrstat[i] <- 2*(logLik(modsim)-logLik(modnsim))
}

#Que tanto coinciden los valores teoricos a los generados por el bootstrap

summary(lrstat)

#Histograma de razones de verosimilitud simuladas
hist(lrstat)
abline(v=obsdev, col="red")

plot(qchisq((1:1000)/1001,3),sort(lrstat),xlab=expression(chi[3]^2), ylab="Cocientes de razon de verosimilitud simulados") #por que estos valores en qchisq?
abline(0,1)

mean(lrstat > obsdev) #Valor de probabilidad mucho mas cercano al valor de F, 
# el quid es si conocemos los grados de libertad del denominador, disenhos balanceados.

# b)
# Regresemos a nuestro ejemplo sobre las horas de suenho adicionales, base de datos de las horas extra de suenho
# m1 <- lme(suext~droga, random=~1|ind, data=sn)

data(sleep)
#help(sleep)

str(sleep)
head(sleep)


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

sn$drogac <- factor(sn$droga)
levels(sn$drogac) <- c("p", "s")

#Grafica 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)
sugraf

#Effects parametrization
m2 <- lmer(suext~ drogac + (1|ind),  data = sn)
summary(m2)
anova(m2)

#Podemos estimar los valores de probabilidad de diferentes maneras

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

mn <- lmer(suext~ 1 +(1|ind), data=sn)
pboot <- function(mn,m2) {
s <- simulate(mn)
L0 <- logLik(refit(mn,s))
L1 <- logLik(refit(m2,s))
2*(L1-L0)
}
obsdev2 <- c(2*(logLik(mn)-logLik(m2)))
psn <- replicate(1000, pboot(mn,m2))

#Prueba de razon de verosimilitud
anova(mn,m2)

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

#Valor obtenido mediante un bootstrap parametrico
mean(psn>abs(obsdev2))

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

## Ejercicio 9.1 Realiza un bootstrap parametrico con la base de datos
## de los coccidos empleada en la practica previa.

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))


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

#Modelo nulo
m0 <- lmer(difins~1 + (1|bloque), cc, REML=FALSE)
#Modelo con la variable predictora
m1a <- lmer(difins~trat + (1|bloque), cc, REML=FALSE)

#Prueba del cociente de verosimilitud
anova(m0, m1a)

#Bootstrap parametrico
lrstat <- numeric(1000)
for(i in 1:1000){
rdifins <- unlist(simulate(m0))
modnsim <- lmer(rdifins ~ 1 + (1|bloque), cc, REML=FALSE)
modsim <- lmer(rdifins ~ trat + (1|bloque), cc,REML=FALSE)
lrstat[i] <- 2*(logLik(modsim)-logLik(modnsim))
}

obsdev2 <- c(2*(logLik(m0)-logLik(m1a)))
obsdev2 <- abs(obsdev2)
obsdev3 <- obsdev2
obsdev3

summary(lrstat)


#Valores de la razon de verosimilitud vs quantiles de la distribucion de la Chi cuadrada
hist(lrstat)
abline(v=obsdev3, col="red")
plot(qchisq((1:1000)/1001,2),sort(lrstat),xlab=expression(chi[2]^2), ylab="Cocientes de razon de verosimilitud simulados") 
abline(0,1)
#Valor de probabilidad
mean(lrstat > obsdev3)