306-apli-rat-pup.Rmd

# Aplicación rat pup {#apli-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 @West2014.

```{r echo=FALSE, out.width="50%", fig.align='center', fig.cap="Camada de ratones, tomada de https://www.howcast.com/videos/509444-how-to-care-for-a-pregnant-rat-pet-rats"}
knitr::include_graphics("images/rat_pups.jpg")
```

A continuación la base de datos a utilizar.

```{r}
library(nlme)
head(RatPupWeight)
```

Esta 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. 

```{r}
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.

```{r}
anova(mod31)
```

El 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`.

```{r}
mod31a <- gls(weight ~ Treatment * sex + Lsize,
              data=RatPupWeight, method = "REML")
```

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$.

```{r}
anova(mod31, mod31a)
```

De 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.

```{r}
mod32a <- lme(weight ~ Treatment * sex + Lsize, random= ~ 1 | Litter,
             data=RatPupWeight, method = "REML",
             weights = varIdent(form = ~1 | Treatment))

summary(mod32a)
```