9 Aplicación rat pup
En este capítulo se mostrará un ejemplo de un modelo mixto con dos niveles usando la base de datos RatPupWeight del paquete nlme. Este ejemplo corresponde al ejemplo de la sección 3.2 de West (2014).
Figure 9.1: Camada de ratones, tomada de https://www.howcast.com/videos/509444-how-to-care-for-a-pregnant-rat-pet-rats
A continuación la base de datos a utilizar.
## Grouped Data: weight ~ 1 | Litter
## weight sex Litter Lsize Treatment
## 1 6.60 Male 1 12 Control
## 2 7.40 Male 1 12 Control
## 3 7.15 Male 1 12 Control
## 4 7.24 Male 1 12 Control
## 5 7.10 Male 1 12 Control
## 6 6.04 Male 1 12 ControlEsta base de datos sobre crecimiento contiene la información sobre altura (weigth), sexo (sex), camada (litter), tamaño de la camada (Lsize) y tratamiento (Treatment) al que se sometió un rata antes de tener ratones.
La estrategia es crear el modelo de forma incremental. Vamos a comenzar la construcción a partir del modelo mod31 que usa como efectos fijos el tratamiento, el sexo, la interacción entre tratamiento y sexo, el tamaño de la camada. Adicionalmente, este modelo incluye intercepto aleatorio debido a la camada.
mod31 <- lme(weight ~ Treatment * sex + Lsize, random= ~ 1 | Litter,
data=RatPupWeight, method = "REML")Vamos a usar la función anova.lme para estudiar la importancia de cada efecto fijo en el modelo.
## numDF denDF F-value p-value
## (Intercept) 1 292 9093.772 <.0001
## Treatment 2 23 5.082 0.0149
## sex 1 292 52.602 <.0001
## Lsize 1 23 47.374 <.0001
## Treatment:sex 2 292 0.466 0.6282El otro modelo inicial mod31a es un modelo lineal sin intercepto aleatorio y se ajusta con la función gls y no con lm, esto se hace para poder comparar los modelos mod31 y mod31a.
Para comparar los modelos mod31 y mod31a se puede usar nuevamente la función anova.lme. En este caso la hipótesis nula es \(H_0: \sigma_{b0}=0\) frente a \(H_1: \sigma_{b0} > 0\).
## Model df AIC BIC logLik Test L.Ratio p-value
## mod31 1 9 421.3015 455.0747 -201.6508
## mod31a 2 8 508.7072 538.7277 -246.3536 1 vs 2 89.40562 <.0001De la salida anterior se tiene que el estadístico de la prueba LRT (likelihood ratio test) es 89.40562 con un valor-P inferior a 0.0001, es quiere decir que hay evidencias en contra de \(H_0\), lo que significa que incluir \(b_0\) mejora el modelo.
mod32a <- lme(weight ~ Treatment * sex + Lsize, random= ~ 1 | Litter,
data=RatPupWeight, method = "REML",
weights = varIdent(form = ~1 | Treatment))
summary(mod32a)## Linear mixed-effects model fit by REML
## Data: RatPupWeight
## AIC BIC logLik
## 384.0819 425.3602 -181.041
##
## Random effects:
## Formula: ~1 | Litter
## (Intercept) Residual
## StdDev: 0.3134846 0.5147948
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Treatment
## Parameter estimates:
## Control High Low
## 1.0000000 0.6394383 0.5649830
## Fixed effects: weight ~ Treatment * sex + Lsize
## Value Std.Error DF t-value p-value
## (Intercept) 7.888771 0.23223913 292 33.96831 0.0000
## Treatment.L -0.638713 0.13587695 23 -4.70067 0.0001
## Treatment.Q 0.011965 0.11635236 23 0.10283 0.9190
## sexFemale -0.351238 0.04682791 292 -7.50062 0.0000
## Lsize -0.130007 0.01848708 23 -7.03233 0.0000
## Treatment.L:sexFemale 0.066939 0.09135485 292 0.73274 0.4643
## Treatment.Q:sexFemale -0.023417 0.06929365 292 -0.33794 0.7357
## Correlation:
## (Intr) Trtm.L Trtm.Q sexFml Lsize Tr.L:F
## Treatment.L -0.325
## Treatment.Q -0.139 0.154
## sexFemale -0.126 -0.008 -0.103
## Lsize -0.954 0.376 0.190 0.031
## Treatment.L:sexFemale -0.020 -0.304 -0.012 -0.034 0.015
## Treatment.Q:sexFemale -0.021 -0.023 -0.300 0.445 -0.017 -0.029
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.88670114 -0.52493419 0.02123518 0.57307286 2.56409983
##
## Number of Observations: 322
## Number of Groups: 27