Práctico 04: Modelos Mixtos

Modelo lineal mixto.

Los datos del archivo floral_size.csv corresponden a un estudio donde se pretende determinar si el tamaño de las flores varía entre poblaciones. En cada población se muestrearon 30 individuos y de cada individuo se seleccionaron 3 flores (ocasionalmente 2).

library(emmeans)
library(nlme)
library(multcomp)
library(ggplot2)

Csize_mF <- read.csv("floral_size.csv", header = TRUE, stringsAsFactors = TRUE)

## Modelo con nlme (normal, único efecto random posible)
Csize.field1 <- lme(log(CSF) ~ PopF, random = ~ 1 | indivF, data = Csize_mF, method = "REML") 
summary(Csize.field1)


## probar efectos fijos
Csize.field0 <- lme(log(CSF) ~ 1, random = ~ 1 | indivF, data = Csize_mF, method = "ML") 
Csize.field2 <- lme(log(CSF) ~ PopF, random = ~ 1 | indivF, data = Csize_mF, method = "ML") 
anova(Csize.field0, Csize.field2)

# modelo final
Csize.field3 <- lme(log(CSF) ~ PopF, random = ~ 1 | indivF, data = Csize_mF, method = "REML") 
summary(Csize.field3)

# diagnóstico
plot(Csize.field3)
qqnorm(resid(Csize.field3))

# pair-wise differences
CSgrafico_f <- emmeans(Csize.field3, list(pairwise ~ PopF), adjust = "tukey")
CSgrafico_f

# Floral shape in field - plot 
CS_campito <- plot(CSgrafico_f)
letrasF <- cld(CSgrafico_f, Letters = letters)
CS_F <- CS_campito + 
  geom_text(data = letrasF, aes(letrasF$emmean, letrasF$PopF, label = letrasF$.group),
            position = position_nudge(x = 0.05), 
            size = 4) +
  labs(x= "log(Centroid size)", y = "Natural populations") + 
  xlim(3.6,4.25) +
  theme_bw()

CS_F

Modelo linal generalizado mixto.

Modelamos el número de gorriones de las marismas (Ammospiza caudacuta) en función de la cobertura de juncos. La variable respuesta es banded (número de gorriones), por lo que probaremos una distribución de Poisson. Se utilizará sitio como variable aleatoria.

library(lme4)
library(ggplot2)

banded <- read.table("bandedSP.txt", header = TRUE)

# exploración
ggplot(banded, aes(x=Juncus,y=Banded))+geom_point()+facet_wrap(.~Site)

# diferentes modelos random (sólo ML)
m1 <- glmer(Banded ~ Juncus + (1|Site), data = banded, 
              family = poisson)
m2 <- glmer(Banded ~ Juncus + (0+Juncus|Site) + (1|Site), data = banded, 
            family = poisson)
# NOTA: la forma (Juncus|Site) no ajusta.

AIC(m1, m2)
BIC(m1, m2)

# diferentes modelos fijos
m3 <- glmer(Banded ~ 1 + (1|Site), data = banded, 
            family = poisson)
m4 <- glmer(Banded ~ Juncus + (1|Site), data = banded, 
            family = poisson)

anova(m3, m4)

summary(m4)

## Debemos chequear la sobredispersión siempre en el modelo más complejo

overdisp_fun <- function(m) {
  rdf <- df.residual(m)
  rp <- residuals(m,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}

overdisp_fun(m2)

# diagnóstico
plot(m4)
qqnorm(resid(m4))

# si hubiera sobredispersión podríamos usar quasipoisson con
library(MASS)

q1 <- glmmPQL(Banded ~ Juncus, random = ~1|Site, data = banded, 
            family = poisson)

summary(q1)

Ejercicios

La base datos Contraception obtenida a partir del paquete mlRev corresponde a un conjuntos de datos del Center for Multilevel Modelling de la Universidad de Bristol y proceden de la Encuesta de Fecundidad de Bangladesh de 1989 (Huq & Cleland, 1990). Los datos son una submuestra de respuestas de 1934 mujeres agrupadas en 60 distritos y poseen varias variables. La idea es ajustar un modelo para ver si existe una relación en el uso de anticonceptivos (uso) y las variables: edad de las personas (age), si las mujeres ya han tenido hijos (ch) y el lugar donde viven (urban). La variable District sería una variable aleatoria. Ajuste el modelo más adecuado. Nota: Para introducir un término cuadrático use poly(x,2).

Pista:

g1 <- ggplot(cont, aes(x = age, y = uso)) + 
  geom_point(alpha = 0.25, position = position_jitter(height = 0.1)) + 
  stat_smooth(method = "gam", method.args = list(family=binomial)) +
  theme_bw() + 
  facet_wrap(.~urban, ncol = 2)
g1