# Capítulo 10. Introducción a los modelos mixtos

#  10.3. Entendiendo los modelos mixtos
#  10.3.1. Un ejemplo con bentos marino

library(ADER)
data(RIKZ)
RIKZ$Richness <- apply(RIKZ[, 2:76] > 0, 1, sum)
RIKZ$Richness <- rowSums(RIKZ[, 2:76] > 0)
RIKZ$Beach <- factor(RIKZ$Beach)

#  10.3.2. Una aproximación al análisis de los datos con modelos
#              lineales

tmp <- lm(Richness ~ NAP, data = RIKZ)
plot(RIKZ$NAP, RIKZ$Richness, ylab="Riqueza de bentos",
  xlab="Nivel intermareal (m)")
abline(tmp, lwd=2)

tmp1<-lm(Richness ~ Beach + NAP, data = RIKZ)

library(visreg)
visreg(tmp1, "NAP", by="Beach", band=F, overlay=T, legend=F,
  pch=16, points.par=list(cex=1), xlab="Nivel intermareal (m)",
  ylab="Riqueza de bentos", ylim=c(0,25))
abline(v=0, lty=3, lwd=2)

tmp2<-lm(Richness ~ NAP + Beach:NAP, data = RIKZ)
visreg(tmp2, "NAP", by="Beach", band=F, overlay=T, legend=F,
  pch=16, points.par=list(cex=1), xlab="Nivel intermareal (m)",
  ylab="Riqueza de bentos", ylim=c(0,25))
abline(v=0, lty=3, lwd=2)

tmp3<-lm(Richness ~ Beach*NAP, data = RIKZ)
  visreg(tmp3, "NAP", by="Beach", band=F, overlay=T, legend=F,
  pch=16, points.par=list(cex=1), xlab="Nivel intermareal (m)",
  ylab="Riqueza de bentos", ylim=c(0,25))
abline(v=0, lty=3, lwd=2)

summary(tmp)$r.squared
summary(tmp1)$r.squared
summary(tmp2)$r.squared
summary(tmp3)$r.squared

anova(tmp3)

#  10.3.3. Reanalizando los datos con modelos lineales mixtos

#  A) El paquete nlme

library(nlme)
lme1 <- lme(Richness ~ NAP, data = RIKZ, random = ~ 1|Beach)
summary(lme1)
anova(lme1)

#  B) El paquete lme4

library(lme4)
library(car)
lmer1 <- lmer(Richness ~ NAP + (1|Beach), data = RIKZ)
summary(lmer1)
Anova(lmer1, test.statistic="F")

#  C) Otros efectos aleatorios en modelos lineales mixtos

library(lme4)
library(car)
lmer3 <- lmer(Richness ~ NAP +(NAP|Beach), data = RIKZ)
summary(lmer3)
Anova(lmer3, test.statistic="F")
coef(lmer3)

library(lattice)
dotplot(ranef(lmer3, condVar=TRUE))

library(lme4)
lmer3_1 <- lmer(Richness ~ NAP +(NAP||Beach), data = RIKZ)
summary(lmer3_1)

#  10.4. Inferencia con modelos mixtos

#  10.4.1. Test de la F

library(lme4)
lmer1 <- lmer(Richness ~ NAP + (1|Beach), data = RIKZ)
anova(lmer1)

library(lmerTest)
lmer1.skr <- lmer(Richness ~ NAP + (1|Beach), data = RIKZ)
anova(lmer1.skr)

library(car)
Anova(lmer1, test.statistic="F")

#  10.4.2. Test de hipótesis para modelos anidados (likelihood ra-
#              tio test)

lme1 <- lme(Richness ~ NAP, data = RIKZ, random = ~ 1|Beach)
lme2 <- lme(Richness ~ NAP, data = RIKZ, random = ~ 1 + NAP|Beach)
anova(lme1, lme2)

lme2 <- lme(Richness ~ NAP, data = RIKZ, random = ~ 1|Beach)
tmp <- gls(Richness ~ NAP, data = RIKZ)
anova(lme2, tmp)

library(pbkrtest)
glmer1 <- glmer(Richness ~ NAP +(1|Beach), data = RIKZ,
  family=poisson)
glmer3 <- glmer(Richness ~ NAP +(1+ NAP|Beach), data = RIKZ,
  family=poisson)
set.seed(2)
PBmodcomp(glmer1, glmer3, nsim=99)

#  10.4.3. Comparación de modelos con AIC
#  A) Generación automática de modelos

fm1 <- lme(Richness ~ NAP + temperature + salinity, data = RIKZ,
  random = ~ 1|Beach, method="ML")

library(MuMIn)
options(na.action = "na.fail")
dd <- dredge(fm1)
dd
subset(dd, delta < 2)


#  10.5. Análisis de los residuos

Res <- residuals(lmer3, type="response")
Fit <- fitted(lmer3)

par(mfrow=c(2,2))
plot(Res ~ Fit, xlab="Fitted values", ylab="Residuals")
abline(h=0)
plot(Res ~ RIKZ$NAP, xlab="NAP", ylab="Residuals")
abline(h=0)
hist(Res, main="", xlab="Residuals")
qqnorm(Res, main="")
qqline(Res)

#  10.6. Estimando componentes de la varianza en modelos
#          mixtos

library(MuMIn)
r.squaredGLMM(lmer3)

lmer4 <- lmer(Richness ~ NAP + salinity + penetrability +
  (1|Beach), data = RIKZ)

library(partR2)
partR2(lmer4, partvars=c("NAP", "salinity", "penetrability"),
  nboot=20)

#  10.7. Modelos lineales generalizados mixtos (GLMM)

glmer3 <- glmer(Richness ~ NAP + (NAP|Beach), data = RIKZ,
  family=poisson)
summary(glmer3)
Res <- residuals(glmer3, type="pearson")
Fit <- fitted(glmer3)

par(mfrow=c(2,2))
plot(Res ~ Fit, xlab="Fitted values", ylab="Residuals")
abline(h=0)
plot(Res ~ RIKZ$NAP, xlab="NAP", ylab="Residuals")
abline(h=0)
hist(Res, main="", xlab="Residuals")
qqnorm(Res, main="")
qqline(Res)

library(MASS)
glmer3.nb <- glmer.nb(Richness ~ NAP + (NAP|Beach), data = RIKZ)

#  10.8. Otras funciones para ajustar modelos mixtos en R

library(lmerTest)
mod.lmeT <- lmer(Richness ~ NAP + (NAP|Beach), data = RIKZ)
summary(mod.lmeT)
anova(mod.lmeT, type=2)
ranova(mod.lmeT)

library(afex)
mod.afex <- mixed(Richness ~ NAP + (NAP|Beach), data = RIKZ, type=2)
mod.afex
summary(mod.afex)

library(glmmTMB)
mod.glmmTMB <- glmmTMB(Richness ~ NAP + (NAP|Beach), data = RIKZ,
  family=poisson, REML=TRUE)
summary(mod.glmmTMB)

library(phyr)
data(comm_a)
data(envi)
data(traits)
data(phylotree)
comm <- comm_a
comm$site <- row.names(comm)

library(data.table)
dat <- data.frame(melt(setDT(comm), id.vars = "site",
  variable.name = "sp", value.name="freq"))
dat <- merge(dat, envi, by="site")
dat <- merge(dat, traits, by="sp")

mod.phyr <- pglmm(freq ~ 1 + shade + (1|sp__) + (1|site) +
  (1|sp__@site), data = dat, family = "gaussian",
  REML = FALSE, cov_ranef = list(sp = phylotree))
summary(mod.phyr)

library(brms)
mod.brm <- brm(Richness ~ NAP + (NAP|Beach), data = RIKZ,
  family=poisson, chains = 2, cores = 2)
summary(mod.brm)

#  10.9. Evitando problemas de ajuste

lmer3 <- lmer(Richness ~ NAP +(NAP|Beach), data = RIKZ)
glmer3 <- glmer(Richness ~ NAP + (NAP|Beach), data = RIKZ,
  family=poisson)
lmer3.s <- lmer(Richness ~ scale(NAP) +(scale(NAP)|Beach),
  data = RIKZ)
lmer3 <- lmer(Richness ~ NAP +(NAP|Beach), data = RIKZ)

strict_tol <- lmerControl(optCtrl=list(xtol_abs=1e-10,
  ftol_abs=1e-10))
lmer3.tol <- update(lmer3, control=strict_tol)

lmer3.all <- allFit(lmer3)
ss <- summary(lmer3.all)
ss$msgs
ss$sdcor
ss$fixef

lmer3 <- lmer(Richness ~ NAP +(NAP|Beach), data = RIKZ)
lmer3.corr <- lmer(Richness ~ NAP +(NAP||Beach), data = RIKZ)