#### Ejercicio 10.1, base de datos melon del paquete agricolae, cuadrados latinos, codigo escrito por Sofia Ochoa.
R

library(agricolae)
data(melon)

#### 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 melon como la variable predictora fija y el rendimiento como la variable de respuesta

m3<-lmer(yield~variety*(1|row)+(1|col),data=melon)

summary(m3)
anova(m3)

plot(fitted(m3), resid(m3))
qqnorm(resid(m3))


# ¿Como sabemos que en realidad es un cuadrado latino? pista: usar la funcion xtabs

xtabs(~row+col+variety,data=melon)
xtabs(~row+col,data=melon)
#existe un caso para cada fila y columna, por que lo que si es cuadros latinos

####que te dicen las desviaciones estandar de las variables aleatorias, realiza seleccion de modelos usando el AIC

#el valor de col es muy bajo y por lo mismo la variación entre individuos de columnas diferentes es bajo
#, por lo que se propone un modelo considerando las filas como la variable aleatoria.


m4<-lmer(yield~variety+(1|row),data=melon)

#y para comparar se utiliza una anova, para obtener un cociente dee razo n de verosimilitud

anova(m3,m4)
library(bbmle)
AICtab(m3,m4,weights=TRUE)

#### Valida el modelo y determina si existen diferencias el rendimiento de acuerdo a la variedad del melon y si es el caso

summary(m4)
anova(m4)

plot(fitted(m4), resid(m4))
qqnorm(resid(m4))
#transformacion a log de res

#### indica los valores predichos para cada variedad y cuales variedades difieren entre si
efectos3<- allEffects(m4, xlevels=list(variety=levels(variety)))
efectos3
summary(efectos3)

########
#Graficas
library(plotrix)
#vectores de valores promedios y minimos y maximos
P<-c(47.25,42.25,49.25,35.00)
LO<-c(37.94098,32.94098,39.94098,25.69098)
UP<-c(56.55902, 51.55902, 58.55902, 44.30902)

plotCI(x=P,ui=UP,li=LO,xlab="Variedad",ylab="Valores predichos")
#grado

#modelo nulo
m5<-lmer(yield~1+(1|row),data=melon)
anova(m4,m5)

pboot <- function(m5,m4) {
s <- simulate(m5)
L0 <- logLik(refit(m5,s))
L1 <- logLik(refit(m4,s))
2*(L1-L0)
}
obsdevM <- c(2*(logLik(m5)-logLik(m4)))
psnM <- replicate(1000, pboot(m5,m4))

pchisq(abs(obsdevM), 3, lower.tail=FALSE)

mean(psnM > abs(obsdevM))

#el valor explica la sobreposicion de los intervalos de confianza