R para datos de biologging

class: center, middle, inverse, title-slide

# R para datos de biologging
## Modelos mixtos en R
### Rocío Joo
### Nov.-Dic. 2020

---

# Objetivos de esta unidad

* Implementar modelos mixtos en R e interpretar resultados
 
---

# Recordando teoría sobre modelos lineales, generalizados y mixtos

Empezaremos con un trabajo en grupo aleatorio sobre teoría de los modelos.

[30 minutos] para discutir y llenar el cuestionario cuyo link les enviaré.

Veamos las respuestas y conversemos. Luego colgaré las diapositivas.

---

# Explorando datos

Utilizaremos datos y algunos códigos de [ourcodingclub](https://ourcodingclub.github.io/tutorials/mixed-models/).

Descargamos primero sus datos de dragones.

```r
library(ggplot2) # plots
library(lme4) # lmer
library(MuMIn) # dredge

path_data <- "./data/"
path_plots <- "./img/"
```

```r
download.file(url = "https://raw.githubusercontent.com/ourcodingclub/CC-Linear-mixed-models/master/dragons.RData",
              destfile = paste0(path_data,"dragons.RData"))
```

---

# Explorando datos

```r
load(paste0(path_data,"dragons.RData"))
```

Los dragones tomaron un test de inteligencia (*testScore*) 
en 3 lugares (*sites*), ubicados en 
diferentes niveles de montañas (*mountainRange*), 
También se tomaron medidas corporales (*bodyLength*)

```r
head(dragons,5)
```

```
##   testScore bodyLength mountainRange  X site
## 1 16.147309   165.5485      Bavarian NA    a
## 2 33.886183   167.5593      Bavarian NA    a
## 3  6.038333   165.8830      Bavarian NA    a
## 4 18.838821   167.6855      Bavarian NA    a
## 5 33.862328   169.9597      Bavarian NA    a
```

```r
str(dragons)
```

```
## 'data.frame':	480 obs. of  5 variables:
##  $ testScore    : num  16.15 33.89 6.04 18.84 33.86 ...
##  $ bodyLength   : num  166 168 166 168 170 ...
##  $ mountainRange: Factor w/ 8 levels "Bavarian","Central",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ X            : logi  NA NA NA NA NA NA ...
##  $ site         : Factor w/ 3 levels "a","b","c": 1 1 1 1 1 1 1 1 1 1 ...
```

---

# Regresión lineal

* Se piensa que hay una relación lineal entre la longitud corporal y 
el nivel de inteligencia. 
* Se hará una regresión lineal

```r
reg_lineal <- lm(formula = testScore ~ bodyLength, data = dragons)
summary(reg_lineal)
```

```
## 
## Call:
## lm(formula = testScore ~ bodyLength, data = dragons)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -56.962 -16.411  -0.783  15.193  55.200 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -61.31783   12.06694  -5.081 5.38e-07 ***
## bodyLength    0.55487    0.05975   9.287  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.2 on 478 degrees of freedom
## Multiple R-squared:  0.1529,	Adjusted R-squared:  0.1511 
## F-statistic: 86.25 on 1 and 478 DF,  p-value: < 2.2e-16
```

Pendiente: cambio en la variable respuesta por unidad de la variable independiente.

---

# Supuestos

* Hay una relación lineal entre variable dependiente e independiente(s).

```r
ggplot(data = dragons, mapping = aes(x = bodyLength, y = testScore)) +
  geom_point() +
  geom_smooth(method = "lm") # esto hace la regresión otra vez 
```

```
## `geom_smooth()` using formula 'y ~ x'
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-6-1.png)

```r
# la línea se podría suplantar con resultados de reg_lineal como
# geom_line(aes(y = reg_lineal$fitted.values))
```

---

# Supuestos

* Los residuales son normales con media 0 y varianza constante

```r
plot(reg_lineal, which = 1) 
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-7-1.png)

* La media no es realmente cero
* La varianza no parece constante

---

# Supuestos

* Los residuales son normales con media 0 y varianza constante

```r
dragons$res <- residuals(reg_lineal)
ggplot(data = dragons, aes(x = res)) +
  geom_histogram() +
  theme_bw()
```

```
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-8-1.png)
 * No parece muy normal, pero tampoco extremadamente lejos

---

# Supuestos

* Los residuales son normales con media 0 y varianza constante

```r
plot(reg_lineal, which = 2) 
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-9-1.png)

* No parece estar tan alejado de la normalidad.

---

# Supuestos

* Combinando información de residuales y ajuste

```r
dragons$predicted <- predict(reg_lineal) # valores predichos del modelo

ggplot(dragons, aes(x = bodyLength, y = testScore)) +
  geom_smooth(method = "lm", se = FALSE, color = "lightgrey") + # línea de regresión
  geom_segment(aes(xend = bodyLength, yend = predicted), alpha = .2) + # distancia entre observado y predicho
  geom_point(aes(color = res)) +  # colorear de acuerdo al tamaño del residual
  scale_color_viridis_c(begin = 0.2, end = 0.8, option = "plasma") +  # nueva escala de colores 
  # <
  geom_point(aes(y = predicted), shape = 1, cex = 0.5) + # dibujar los puntos predichos encima
  theme_bw()
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-10-1.png)

---

# Ejemplos problemáticos

Los tomamos prestados de [Applied Statistics with R](https://daviddalpiaz.github.io/appliedstats/model-diagnostics.html#unusual-observations)

Varianza no constante

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-12-1.png)
![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-13-1.png)

¿Qué hacer aquí?

---

# Ejemplos problemáticos

Los tomamos prestados de [Applied Statistics with R](https://daviddalpiaz.github.io/appliedstats/model-diagnostics.html#unusual-observations)

No linearidad

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-15-1.png)
![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-16-1.png)

¿Qué hacer aquí?

---

# Ejemplos problemáticos

Los tomamos prestados de [Applied Statistics with R](https://daviddalpiaz.github.io/appliedstats/model-diagnostics.html#unusual-observations)

No normalidad

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-17-1.png)

¿Qué hacer aquí?

---

# Outliers

Adaptamos ejemplos de [Applied Statistics with R](https://daviddalpiaz.github.io/appliedstats/model-diagnostics.html#unusual-observations).

Empezamos generando datos y modelándolos con una regresión lineal

```r
set.seed(42)
ex_data  = data.frame(x = 1:10,
                      y = 10:1 + rnorm(n = 10))
ex_model = lm(y ~ x, data = ex_data)
```

---

# Outliers

Añadamos un punto (x, y) y reajustemos una regresión lineal

```r
point_1 = c(5.4, 11) # nuevo punto
ex_data_1 = rbind(point_1, ex_data) # añadiéndolo
model_1 = lm(y ~ x, data = ex_data_1) # nuevo modelo

ggplot(data = ex_data_1, aes(x = x, y = y)) +
  geom_point(color = "#999999") + # plotting all points
  geom_point(aes(x = x[1], y = y[1]), cex = 4, shape = 1) + # circunferencia para el nuevo punto
  geom_smooth(data = ex_data, method = "lm", se = FALSE, aes(color = "inicial")) + # línea de regresión del modelo inicial
  geom_smooth(method = "lm", se = FALSE, aes(color = "nuevo"), linetype = 2) + # línea de regresión del nuevo modelo
  scale_color_manual(name = "modelos", labels = c("inicial", "nuevo"), values = c("#56B4E9", "#E69F00")) +
theme_bw()
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-19-1.png)
El punto "extraño" no genera mayores cambios en el modelo.

---

# Outliers

Otro punto (x, y):

```r
point_2 = c(18, -5.7) # nuevo punto
ex_data_2 = rbind(point_2, ex_data) # añadiéndolo
model_2 = lm(y ~ x, data = ex_data_2) # nuevo modelo

ggplot(data = ex_data_2, aes(x = x, y = y)) +
  geom_point(color = "#999999") + # plotting all points
  geom_point(aes(x = x[1], y = y[1]), cex = 4, shape = 1) + # circunferencia para el nuevo punto
  geom_smooth(data = ex_data, method = "lm", se = FALSE, aes(color = "inicial")) + # línea de regresión del modelo inicial
  geom_smooth(method = "lm", se = FALSE, aes(color = "nuevo"), linetype = 2) + # línea de regresión del nuevo modelo
  scale_color_manual(name = "modelos", labels = c("inicial", "nuevo"), values = c("#56B4E9", "#E69F00")) +
theme_bw()
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-20-1.png)
El nuevo punto no influencia el modelo.

---

# Outliers

Otro punto (x, y):

```r
point_3 = c(14, 5.1) # nuevo punto
ex_data_3 = rbind(point_3, ex_data) # añadiéndolo
model_3 = lm(y ~ x, data = ex_data_3) # nuevo modelo

ggplot(data = ex_data_3, aes(x = x, y = y)) +
  geom_point(color = "#999999") + # plotting all points
  geom_point(aes(x = x[1], y = y[1]), cex = 4, shape = 1) + # circunferencia para el nuevo punto
  geom_smooth(data = ex_data, method = "lm", se = FALSE, aes(color = "inicial")) + # línea de regresión del modelo inicial
  geom_smooth(method = "lm", se = FALSE, aes(color = "nuevo"), linetype = 2) + # línea de regresión del nuevo modelo
  scale_color_manual(name = "modelos", labels = c("inicial", "nuevo"), values = c("#56B4E9", "#E69F00")) +
theme_bw()
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-21-1.png)
El nuevo modelo parece bastante distinto

---

# Supuestos

* ¿Por qué es importante chequearlos?

Porque en ellos se basan los modelos. Importante siempre. 
Necesario si el modelo es para explicar.

* ¿Por qué no usar tests?

Se puede, pero no hay que confiar a ciegas en un p-valor arbitrario. Ver referencias.

---

# Desafío

* En grupos aleatorios, tratar de identificar si hay algún punto que esté influyendo mucho
sobre la regresión ajustada a los datos de dragones.

* Dentro del GD: CURSO R-biologging/R-Rocío/04-Modelos-Mixtos/Actividades/Dragones/, 
  * Crear una carpeta ahí con sus iniciales, e.g. "CI-CZ"
  * Un google doc con sus nombres
  * Un script R
  * Archivo(s) de gráfico

[25 minutos]

---

# Modelos mixtos

* Otro escenario: se piensa que el nivel de inteligencia estaría explicado por 
la longitud corporal, pero que el resultado podría estar afectado por el nivel de 
montaña en el que se tomó cada muestra.

* ¿Qué corresponde hacer aquí?

* Intentaremos explicar (buena) parte de la variación del nivel de inteligencia
por la longitud corporal (efecto fijo).

* Parte de la variación residual estaría asociada al nivel de montaña (efecto 
aleatorio).

* Intentaríamos modelar esta relación entre inteligencia y longitud corporal
controlando la variación debida a nivel de montaña.

```r
mod_mix <- lmer(formula = testScore ~ bodyLength + (1|mountainRange), 
                data = dragons)
```

---

```r
summary(mod_mix)
```

```
## Linear mixed model fit by REML ['lmerMod']
## Formula: testScore ~ bodyLength + (1 | mountainRange)
##    Data: dragons
## 
## REML criterion at convergence: 3991.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4815 -0.6513  0.0066  0.6685  2.9583 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  mountainRange (Intercept) 339.7    18.43   
##  Residual                  223.8    14.96   
## Number of obs: 480, groups:  mountainRange, 8
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 43.70938   17.13489   2.551
## bodyLength   0.03316    0.07865   0.422
## 
## Correlation of Fixed Effects:
##            (Intr)
## bodyLength -0.924
```

---

# Modelos mixtos

```r
dragons$predicted_modmix <- predict(mod_mix) # valores predichos del modelo

ggplot(dragons, aes(x = bodyLength, y = testScore, color = mountainRange)) +
  geom_point(alpha = 0.5) +
  geom_line(aes(y = predicted_modmix), size = 1.5) +
  scale_color_brewer(palette = "Set2") +
  theme_bw()
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-24-1.png)

Efectivamente no se observa una relación fuerte de testScore con bodyLength 
después de haber controlado por mountainRange.

---

# Modelos mixtos

Residuales vs. valores predichos

```r
plot(mod_mix, which = 1)
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-25-1.png)

No se observan problemas de linealidad, media constante o varianza constante.

---

# Modelos mixtos

```r
qqnorm(resid(mod_mix)) 
qqline(resid(mod_mix))  
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-26-1.png)

O de normalidad. Quizás algún outlier.

¿Cómo interpretar estos resultados? Discutamos.

---

# Modelos mixtos

Un ejemplo más real que el de los dragones, de [Clay et al. (2020)](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/1365-2656.13267).

```r
path_file <- paste0(path_data, "BothSites_forwindexperienced.csv")
```

```r
download.file(url = "https://raw.githubusercontent.com/tommyclay/alba-wind-behaviour/master/Data_inputs/BothSites_forwindexperienced.csv", destfile = path_file)
```

```r
dat <- read.csv(file = path_file, na.strings = "NA")
head(dat)
```

```
##   ID TripID   Site Sex Year   WindSp  WindDir
## 1  9    269 Crozet   F 2010 9.100572 157.8553
## 2  9    269 Crozet   F 2010 9.100572 157.8553
## 3  9    269 Crozet   F 2010 9.475041 157.7442
## 4  9    269 Crozet   F 2010 9.475041 157.7442
## 5  9    269 Crozet   F 2010 9.475041 157.7442
## 6  9    269 Crozet   F 2010 9.475041 157.7442
```

---

# Modelos mixtos

```r
summary(dat)
```

```
##        ID            TripID                 Site        Sex       
##  Min.   :  1.0   Min.   :  1.0   Crozet       :368698   F:199930  
##  1st Qu.: 63.0   1st Qu.: 97.0   South Georgia: 25454   M:194222  
##  Median :126.0   Median :196.0                                    
##  Mean   :130.2   Mean   :196.3                                    
##  3rd Qu.:200.0   3rd Qu.:299.0                                    
##  Max.   :259.0   Max.   :385.0                                    
##       Year          WindSp           WindDir       
##  Min.   :2010   Min.   : 0.0451   Min.   :-180.00  
##  1st Qu.:2013   1st Qu.: 6.4381   1st Qu.:  13.33  
##  Median :2013   Median : 8.7915   Median :  80.39  
##  Mean   :2014   Mean   : 8.8056   Mean   :  54.00  
##  3rd Qu.:2016   3rd Qu.:11.1239   3rd Qu.: 122.06  
##  Max.   :2016   Max.   :23.6461   Max.   : 180.00
```

El resumen ayuda a ver fácilmente las clases de datos y rangos de valores.

* Hay una línea por punto GPS. 
* No hay lon, lat o tiempo porque no son relevantes aquí.

---

# Modelos mixtos

Transformando a factor `ID`, `TripID` y `Year`

```r
var_fac <- c("ID", "TripID", "Year") # variables a transformar
dat[var_fac] <- lapply(dat[,var_fac], factor) # aplicar función factor a esas variables (columnas)

str(dat)
```

```
## 'data.frame':	394152 obs. of  7 variables:
##  $ ID     : Factor w/ 259 levels "1","2","3","4",..: 9 9 9 9 9 9 9 9 9 9 ...
##  $ TripID : Factor w/ 385 levels "1","2","3","4",..: 269 269 269 269 269 269 269 269 269 269 ...
##  $ Site   : Factor w/ 2 levels "Crozet","South Georgia": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Sex    : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Year   : Factor w/ 7 levels "2010","2011",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ WindSp : num  9.1 9.1 9.48 9.48 9.48 ...
##  $ WindDir: num  158 158 158 158 158 ...
```

---

# Modelos mixtos

De [Clay et al. (2020)](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/1365-2656.13267):

* *We first assessed whether **wind speeds** experienced by birds
varied by **population** and **sex***

```r
WindSp ~ Site + Sex
```

(Site = population)

* *Wind **speed** was modelled as the response variable with a Gaussian error 
distribution, and the factors **population**, **sex**, their two-way 
**interaction** and **year** included as covariates.*

```r
lmer(WindSp ~ Site*Sex + Year)
```

* *A **random effect** of **trip identity** nested within **individual identity** 
was included to account for variation in the number and duration of trips per 
individual respectively*

```r
lmer(WindSp ~ Site*Sex + Year + (1|ID/TripID))
```

---

# Modelos mixtos

```r
mod_lin <- lmer(WindSp ~ Site*Sex + Year + (1|ID/TripID), data = dat, REML = TRUE)
```

* *We performed multi-model inference on the full set of 
predictor combinations and assessed the best-supported model as 
that with the lowest Akaike Information Criterion (AIC).*

* *As small differences in AIC are not considered to be meaningful (Burnham &
Anderson, 2002), if multiple models were within two AIC units, the
best model was deemed to be that which had the fewest number
of parameters (i.e. the most parsimonious) (Harrison et al., 2018).*

```r
mod_lin_mv <- lmer(WindSp ~ Site*Sex + Year + (1|ID/TripID), 
                   data = dat, REML = FALSE)
# para comparar modelos es mejor hacerlo estimando por máxima verosimilitud
options(na.action = "na.fail") # sin esto nos da error "global model na action argument not set" 
# - así nos aseguramos de que la siguiente función nos arroje resultados si no hay
# valores faltantes. 
m_set <- dredge(mod_lin_mv) # Genera una tabla de selección de modelos con combinaciones de términos fijos
```

---

# Modelos mixtos

```r
m_set
```

```
## Global model call: lmer(formula = WindSp ~ Site * Sex + Year + (1 | ID/TripID), 
##     data = dat, REML = FALSE)
## ---
## Model selection table 
##    (Int) Sex Sit Yer Sex:Sit df   logLik    AICc delta weight
## 16 8.364   +   +   +       + 13 -1018031 2036088  0.00   0.91
## 8  8.413   +   +   +         12 -1018035 2036093  4.62   0.09
## 6  8.421   +       +         11 -1018044 2036109 20.74   0.00
## 7  8.702       +   +         11 -1018047 2036116 27.78   0.00
## 5  8.702           +         10 -1018055 2036130 41.48   0.00
## 12 8.474   +   +           +  7 -1018061 2036136 47.95   0.00
## 2  8.523   +                  5 -1018064 2036138 49.81   0.00
## 4  8.539   +   +              6 -1018064 2036140 51.33   0.00
## 1  8.906                      4 -1018074 2036157 68.35   0.00
## 3  8.921       +              5 -1018074 2036158 70.10   0.00
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | ID/TripID'
```

El mejor modelo es el modelo completo.

---

# Modelos mixtos

```r
plot(mod_lin, which = 1)
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-38-1.png)
No se observan problemas de linealidad, media constante o varianza constante.

---

# Modelos mixtos

```r
qqnorm(resid(mod_lin)) 
qqline(resid(mod_lin))  
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-39-1.png)

No tan alejado de la normalidad.

---

# Modelos mixtos

```r
summary(mod_lin)
```

```
## Linear mixed model fit by REML ['lmerMod']
## Formula: WindSp ~ Site * Sex + Year + (1 | ID/TripID)
##    Data: dat
## 
## REML criterion at convergence: 2036073
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4755 -0.6795  0.0071  0.6653  4.6705 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  TripID:ID (Intercept)  1.5963  1.2634  
##  ID        (Intercept)  0.1943  0.4408  
##  Residual              10.2058  3.1947  
## Number of obs: 394152, groups:  TripID:ID, 385; ID, 259
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)              8.3655     0.4116  20.322
## SiteSouth Georgia       -0.8036     0.4259  -1.887
## SexM                     0.9017     0.1535   5.873
## Year2011                -1.0823     0.4551  -2.378
## Year2012                 1.3849     0.4834   2.865
## Year2013                 0.2590     0.4267   0.607
## Year2014                 0.2139     0.5103   0.419
## Year2015                -0.4569     0.5407  -0.845
## Year2016                 0.1356     0.4309   0.315
## SiteSouth Georgia:SexM  -1.1939     0.4684  -2.549
## 
## Correlation of Fixed Effects:
##             (Intr) StSthG SexM   Yr2011 Yr2012 Yr2013 Yr2014 Yr2015 Yr2016
## SiteSothGrg -0.024                                                        
## SexM        -0.137  0.173                                                 
## Year2011    -0.879 -0.012 -0.062                                          
## Year2012    -0.831 -0.332 -0.036  0.759                                   
## Year2013    -0.939 -0.009 -0.044  0.861  0.808                            
## Year2014    -0.786 -0.010 -0.037  0.722  0.678  0.766                     
## Year2015    -0.731 -0.013 -0.074  0.681  0.634  0.720  0.606              
## Year2016    -0.930 -0.009 -0.051  0.853  0.801  0.906  0.759  0.712       
## StStGrg:SxM  0.045 -0.602 -0.328  0.020  0.012  0.014  0.012  0.024  0.017
```

---

# Modelos mixtos

A ver, el mejor está compuesto por año, población y sexo como covariables 
(e interacción población y sexo) en efectos fijos.

¿Qué tipo de gráfico nos ayudaría a interpretar las relaciones modeladas?

* Haremos una predicción de la velocidad del viento para cada combinación de
valores de sexo, población y año, las 3 covariables fijas.

```r
pred_df <- expand.grid(Sex = levels(dat$Sex), Site = levels(dat$Site), 
                       Year = levels(dat$Year)) 
prediccion <- predict(mod_lin, newdata = pred_df, re.form = NA) # re.form No efectos aleatorios
pred_df$fit <- prediccion
```

---

# Modelos mixtos

* Para mantenerlo simple, hacemos un boxplot de la velocidad del viento en función de población y sexo.

```r
ggplot(data = pred_df, aes(x = Site, y = fit, color = Sex)) +
  geom_boxplot() +
  scale_color_brewer(palette = "Dark2") +
  ylab("Velocidad del viento predicha") +
  theme_bw()
```

![](04-Modelos-Mixtos_files/figure-html/unnamed-chunk-42-1.png)

---