Tarea para entregar: caso dieta de Tyto Alba (problema 2)
Biometría 2, Módulo 2 - FCEyN, UBA, 2019
Juan Manuel Berros
October 16, 2019
Problema 2 – Dieta de Tyto alba
La lechuza de campanarios (Tyto alba) es una rapaz que se alimenta principalmente de pequeños mamíferos. Sin embargo, puede incorporar en menores proporciones aves, reptiles y anfibios. Un grupo de investigadores quiere estudiar cómo varía el consumo de aves por parte de T. alba entre distintas localidades de la provincia de Entre Ríos. Para ello analizan las egagrópilas de lechuzas (cúmulo de pelos, plumas, huesos y demás restos no digeribles que estas rapaces regurgitan en los nidos) colectadas en 41 localidades de la provincia, 12 colectadas en la región fitogeográfica de Delta e Islas del Paraná, 11 en la región Pampa y 18 en el Espinal. Además los investigadores tienen información sobre el porcentaje de suelo cubierto por construcciones (% urbano) en cada sitio de muestreo del ambiente.
El objetivo del estudio es determinar si el porcentaje de aves consumidas por T. alba varía entre las regiones fitogeográficas y el nivel de urbanización. Base de datos
Tyto.txt.
Estadística descriptiva
Leemos el dataset y agregamos una nueva variable: la proporción de presas que son aves en cada observación.
tyto <-
read_table2("tyto.txt") %>%
mutate(
Region = factor(Region),
proporcion = N_aves/N_presas
)
tyto ID Region N_presas N_aves porc_urbano proporcion
Min. : 1 Delta :12 Min. :102.0 Min. : 0.00 Min. :0.00000 Min. :0.00000
1st Qu.:11 Espinal:18 1st Qu.:214.0 1st Qu.: 5.00 1st Qu.:0.00000 1st Qu.:0.01623
Median :21 Pampa :11 Median :316.0 Median : 9.00 Median :0.00000 Median :0.02961
Mean :21 Mean :332.8 Mean :13.51 Mean :0.07976 Mean :0.04080
3rd Qu.:31 3rd Qu.:418.0 3rd Qu.:20.00 3rd Qu.:0.12000 3rd Qu.:0.05908
Max. :41 Max. :762.0 Max. :70.00 Max. :0.40000 Max. :0.14019 ¿Cuál es la variable respuesta?
La variable de respuesta es N_aves, el número de presas que son aves, inferido a partir de las egagrópilas de lechuzas. Ese número cobra sentido en relación a una cantidad total de N_presas.
¿Qué valores toma y cómo es su distribución?
N_aves es una variable cuantitativa discreta que toma valores no negativos acotados entre 0 y el total de aves (N_presas) en cada sitio. Puede modelarse con la distribución binomial, como el número de éxitos en un número total de ensayos.
¿Cuáles son los parámetros de la misma?
La distribución binomial toma dos parámetros, el total de ensayos y la cantidad de éxitos, usualmente representados como \((n, \pi)\).
Describa el tipo de variables explicatorias involucradas.
Tenemos dos variables explicatorias:
porc_urbano, el porcentaje de suelo cubierto con construcciones, una medida de cuán urbanizada está el área. Se trata de una variable cuantitativa continua que toma valores entre 0 y 1 (es una proporción).Region, la región de recolección de datos, una variable cualitativa o factor con tres niveles:Delta,PampayEspinal.
Explore si la proporción de aves consumidas se asocia con ambiente muestreado.
fig_box <-
ggplot(tyto) +
aes(x = Region, y = proporcion, color = Region) +
stat_summary(
fun.data = mean_sdl,
fun.args = list(mult = 1), # mult es de cuántos SD dibuja el rango
position = position_nudge(.25)
) +
geom_boxplot(
width = .2,
position = position_nudge(-.25)
) +
geom_jitter(
shape = 1,
size = 2,
width = .05
) +
labs(
subtitle = "Cuantiles, datos y media +- desvío estándar"
) +
guides(color = F)
fig_scat <-
ggplot(tyto) +
aes(x = porc_urbano, y = proporcion, color = Region) +
geom_point() +
geom_label_repel(
data = filter(tyto, ID == 39),
mapping = aes(label = ID)
)
# Esto funciona con library(patchwork)
fig_box + fig_scat + plot_layout(widths = c(6, 3)) +
plot_annotation(title = "Fig. 1 Distribución de N_aves/N_presas en cada localidad",
tag_levels = "A")
La figura 1 A sugiere que hay menos proporción de aves consumidas en la región Espinal y una proporción intermedia en la región Pampa.
En la figura 1 B se marca el dato de ID 39, que luego será relevante en el análisis.
Estadísticos por región:
Construya un modelo lineal generalizado para predecir la proporción de aves consumidas en función del ambiente y el nivel de urbanización.
El modelo elegido (binomial) tiene la función de enlace logit o log odds. Planteamos primero el modelo completo, con interacciones:
\[ \log \left( \frac{\pi_i}{1 - \pi_i} \right) = \log \text{odds}_i = \eta_i = \beta_D + \beta_E \text{ESPINAL}_i + \beta_P \text{PAMPA}_i + \beta_{D:U} U_i + \beta_{E:U} \text{ESPINAL}_i U_i + \beta_{P:U} \text{PAMPA}_i U_i \tag{1} \]
donde las variables ESPINAL y PAMPA son dummies del factor Region y la variable \(U\) es el porcentaje de urbanización.
m1 <- glm(proporcion ~ Region * porc_urbano, family = binomial, weights = N_presas, data = tyto)
# summary(m1)Explicite los 3 componentes del GLM.
El componente aleatorio es la variable de respuesta de distribución binomial, el número de N_aves en el total de N_presas.
El predictor lineal es el descripto en la ecuación \((1)\).
La función de enlace es el logaritmo de los odds.
¿De qué maneras puede introducir la proporción de aves consumidas en el modelo?
Creamos una nueva variable, N_aves/N_presas.
Seleccione el mejor modelo, teniendo en cuenta los supuestos y la ocurrencia de sobre/sub-dispersión.
Calculamos a continuación el estadístico de sobredispersión:
\[ \hat{\phi} = \frac{ \sum^{n}_{i=1} e^2_{i, \text{Pearson}} }{ \text{g.l. residuos }} \]
calcular_sobredispersion <- function(modelo) {
residuos_cuadrados <- resid(modelo, type = "pearson")^2
sum(residuos_cuadrados/modelo$df.residual) %>% round(2)
}
calcular_sobredispersion(m1)[1] 6.42Siguiendo el criterio de Zuur et al. (2009), diremos que hay sobredispersión pues el estadístico \(\hat{\phi} > 1.5\).
Antes de cambiar el modelo, comprobemos si hay outliers culpables de la sobredispersión observada.
tyto %<>%
mutate(
residuos_pearson = residuals(m1, type = "pearson"),
predichos_VR = predict(m1, type = "response"),
predichos_link = predict(m1, type = "link"),
cooks_D = cooks.distance(m1)
)fig_resid <-
ggplot(tyto) +
aes(x = predichos_VR, y = residuos_pearson, color = Region) +
geom_hline(yintercept = 0) +
geom_point() +
geom_label_repel(
data = filter(tyto, cooks_D == max(cooks_D)),
mapping = aes(label = ID)
)
fig_cook <-
ggplot(tyto) +
aes(x = ID, y = cooks_D, color = Region) +
geom_segment(aes(xend = ID, yend = 0)) +
geom_point()
fig_resid + fig_cook +
plot_annotation(title = "Fig. 2 Regresión con el modelo completo, binomial",
tag_levels = "A")
En base a la distancia de Cook, notamos que el dato de ID \(39\) es muy influyente en la regresión, si bien su residuo no está alejado del cero. Realizaremos el ajuste nuevamente, sin ese dato, para comprobar si la sobredispersión se resuelve de ese modo.
tyto_sin_outlier <- filter(tyto, ID != 39)
m2 <- glm(proporcion ~ Region * porc_urbano, family = binomial, weights = N_presas, data = tyto_sin_outlier)
calcular_sobredispersion(m2)[1] 6.28Como vemos, el problema de sobredispersión no desaparece. Además, volviendo al scatterplot de la Fig. 1 A, notamos que el dato sólo es “atípico” porque es el único dato de la región Delta con un porcentaje urbano mayor a cero, pero no parece ser un error y su remoción no se justificaría. Mas aún, su proporción de aves no es extrema. Todo indica que debemos conservarlo.
Para resolver la sobredispersión, intentaremos con la técnica OLRE (observation level random effects).
tyto %<>% mutate(ID = factor(ID))
m3 <- glmmTMB::glmmTMB(proporcion ~ Region * porc_urbano + (1|ID),
weights = N_presas, family = binomial, data = tyto)
calcular_sobredispersion(m3)[1] 0Como vemos, toda la varianza extra es “absorbida” por el efecto aleatorio artificial que agregamos con la variable ID. Conservamos este modelo. Nos interesa ahora saber si hay interacción entre la Region y el porc_urbano, es decir, si el porcentaje de urbanización afecta de manera diferente a la proporción de aves en las presas según la región considerada. Formalmente, consiste en testear la hipótesis
\[ H_0: \beta_{D:U} = \beta_{E:U} = \beta_{P:U} \]
donde la hipótesis alternativa es que al menos uno de esos coeficientes difiere del resto.
Single term deletions
Model:
proporcion ~ Region * porc_urbano + (1 | ID)
Df AIC LRT Pr(>Chi)
<none> 285.75
Region:porc_urbano 2 282.91 1.1615 0.5595El comando drop1 quita variables respetando el principio de marginalidad, de modo que comienza con las interacciones. Vemos que al quitar el efecto de interacción entre ambas variables, la verosimilitud del modelo no cambia significativamente. Es decir, introducir coeficientes de interacción no mejora el ajuste significativamente, de modo que podemos plantear un modelo más simple, aditivo:
\[ \eta_i = \beta_D + \beta_E \text{ESPINAL}_i + \beta_P \text{PAMPA}_i + \beta_U U_i \tag{2} \]
donde hay un solo coeficiente asociado a la variable porc_urbano, \(\beta_U\).
Estime la magnitud del efecto de la/s variables predictoras seleccionadas.
A continuación, planteo el modelo de \((2)\) con una “distribución” cuasi binomial (dejamos el estudio del modelo con OLRE para un trabajo futuro).
m4 <- glm(proporcion ~ Region + porc_urbano,
family = quasibinomial, weights = N_presas, data = tyto)
summary(m4)
Call:
glm(formula = proporcion ~ Region + porc_urbano, family = quasibinomial,
data = tyto, weights = N_presas)
Deviance Residuals:
Min 1Q Median 3Q Max
-4.3096 -2.0803 -0.3761 1.6889 3.7720
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.56988 0.15113 -17.005 < 2e-16 ***
RegionEspinal -1.23203 0.29432 -4.186 0.000168 ***
RegionPampa -0.74154 0.31545 -2.351 0.024177 *
porc_urbano -0.09567 1.16531 -0.082 0.935013
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasibinomial family taken to be 6.225518)
Null deviance: 390.32 on 40 degrees of freedom
Residual deviance: 246.11 on 37 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 5Como vemos, el coeficiente asociado a porc_urbano no es significativamente distinto de cero. Es decir, alcanzaría con la Region para predecir la proporción de aves consumidas. Lo mismo se evidencia al notar que el intervalo de confianza para el coeficiente asociado a porc_urbano incluye el cero:
2.5 % 97.5 %
(Intercept) -2.88 -2.29
RegionEspinal -1.83 -0.67
RegionPampa -1.38 -0.14
porc_urbano -2.48 2.12Por ende, esperamos que drop1 nos sugiera remover esta variable. Dado que estamos en un modelo cuasi binomial, usamos el test \(F\) para la comparación de modelos:
Single term deletions
Model:
proporcion ~ Region + porc_urbano
Df Deviance F value Pr(>F)
<none> 246.11
Region 2 369.13 9.2472 0.0005535 ***
porc_urbano 1 246.15 0.0063 0.9370256
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Efectivamente, podemos quitar la variable porc_urbano. Planteamos entonces el modelo que sólo usa la región como covariable:
\[ \eta_i = \beta_D + \beta_E \text{ESPINAL}_i + \beta_P \text{PAMPA}_i \tag{3} \]
Call:
glm(formula = proporcion ~ Region, family = quasibinomial, data = tyto,
weights = N_presas)
Deviance Residuals:
Min 1Q Median 3Q Max
-4.3769 -2.0656 -0.3743 1.7165 3.7179
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.5704 0.1488 -17.269 <2e-16 ***
RegionEspinal -1.2399 0.2747 -4.513 6e-05 ***
RegionPampa -0.7563 0.2563 -2.951 0.0054 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasibinomial family taken to be 6.050867)
Null deviance: 390.32 on 40 degrees of freedom
Residual deviance: 246.15 on 38 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 5Comprobamos que la Region es significativa:
Single term deletions
Model:
proporcion ~ Region
Df Deviance F value Pr(>F)
<none> 246.15
Region 2 390.32 11.128 0.000157 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Por ende, nos quedamos con el modelo resultante, con los siguientes coeficientes \(\beta_j\) y sus intervalos de confianza 95%:
coeficientes <-
confint(m5) %>%
as_tibble(rownames = "Region") %>%
mutate(
Estimado = coefficients(m5)
) %>%
select(Region, `2.5 %`, Estimado, `97.5 %`) %>%
mutate(Region = str_replace(Region, "\\(Intercept\\)", "Delta"),
Region = str_replace(Region, "Region", ""))Waiting for profiling to be done...El modelo elegido entonces es:
\[ \eta_i = \text{logit}_i = -2.57 + (-1.24) \cdot \text{ESPINAL}_i + (-0.76) \cdot \text{PAMPA}_i \tag{4} \] De aquí podemos calcular el logit o log odds de cada región:
\[ \text{logit}_D = \eta_D = \beta_D = -2.57 \\ \text{logit}_E = \eta_E = \beta_D + \beta_E = -3.81 \\ \text{logit}_P = \eta_P = \beta_D + \beta_P = -3.33 \]
Alternativamente, emmeans lo hace por nosotros y nos da los IC de nivel asintótico 95%:
Region emmean SE df asymp.LCL asymp.UCL
Delta -2.57 0.149 Inf -2.86 -2.28
Espinal -3.81 0.231 Inf -4.26 -3.36
Pampa -3.33 0.209 Inf -3.74 -2.92
Results are given on the logit (not the response) scale.
Confidence level used: 0.95 A continuación obtenemos los odds estimados para cada localidad,
\[ \text{odds}_D = e^{\beta_D} = e^{-2.57} = 0.08 \\ \text{odds}_E = e^{\beta_D} e^{\beta_E} = e^{-3.81} = 0.02 \\ \text{odds}_P = e^{\beta_D} e^{\beta_P} = e^{-3.33} = 0.04 \] con sus IC de confianza 95%:
odds_IC <-
logits_IC %>%
as_tibble() %>%
mutate_if(is.double, ~ round(exp(.), 2)) %>%
select(-SE, -df)
odds_ICPor otro lado, \(e\) elevado a cada \(\beta_j\) corresponde con la magnitud del efecto u odds ratio de cada región:
\[ \text{OR}_D = e^{\beta_D} = e^{-2.57} = 0.08 \\ \text{OR}_E = e^{\beta_E} = e^{-1.24} = 0.29 \\ \text{OR}_P = e^{\beta_P} = e^{-0.76} = 0.47 \]
A continuación se muestran estos mismos ORs junto a sus IC de nivel 95%:
Calcule (a mano) la probabilidad estimada de consumo de aves para la región Espinal
Finalmente, podemos despejar la probabilidad de que una presa consumida sea un ave en cada región, es decir, cada valor de \(\pi\):
\[ \pi_D = \frac{e^{\beta_D}}{1 + e^{\beta_D}} = 0.071 \\ \pi_E = \frac{e^{\beta_D + \beta_E}}{1 + e^{\beta_D + \beta_E}} = 0.022 \\ \pi_P = \frac{e^{\beta_D + \beta_P}}{1 + e^{\beta_D + \beta_P}} = 0.035 \]
Como se ve, la probabilidad estimada para la región Espinal es 0.022.
Nuevamente, eemeans calcula estas probabilidades por nosotros con sus intervalos de confianza:
Region prob SE df asymp.LCL asymp.UCL
Delta 0.0711 0.00983 Inf 0.0541 0.0929
Espinal 0.0217 0.00489 Inf 0.0139 0.0336
Pampa 0.0347 0.00698 Inf 0.0233 0.0513
Confidence level used: 0.95
Intervals are back-transformed from the logit scale Sintetice en un gráfico las estimaciones del modelo seleccionado. Concluya en forma general e Informe sus resultados.
contrast estimate SE df z.ratio p.value
Delta - Espinal 1.240 0.275 Inf 4.513 <.0001
Delta - Pampa 0.756 0.256 Inf 2.951 0.0089
Espinal - Pampa -0.484 0.311 Inf -1.554 0.2659
Results are given on the log odds ratio (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates contrast odds.ratio SE df z.ratio p.value
Delta / Espinal 3.455 0.949 Inf 4.513 <.0001
Delta / Pampa 2.130 0.546 Inf 2.951 0.0089
Espinal / Pampa 0.617 0.192 Inf -1.554 0.2659
P value adjustment: tukey method for comparing a family of 3 estimates
Tests are performed on the log odds ratio scale Una manera es graficar las probabilidades estimadas por región:
fig_link <-
plot(contrastes, comparison = T) +
geom_hline(yintercept = 0) +
geom_text(aes(label = round(the.emmean, 3)), nudge_y = .2) +
labs(
subtitle = "Log odds de consumir un ave",
x = "Media Marginal Estimada (E.M.Mean)"
)
fig_vr <-
plot(contrastes_VR, comparison = T) +
geom_text(aes(label = round(the.emmean, 3)), nudge_y = .2) +
labs(subtitle = "Probabilidad de consumir un ave")
fig_link + fig_vr + plot_annotation(title = "Fig. 3 Estimaciones", tag_level = "A")
En la figura 3 B se ve que en la región Delta hay mayor probabilidad de que las presas consumidas sean aves con respecto a las regiones Pampa y Espinal. Esto se corresponde a que en la figura 3 A se vea un mayor logit en el Delta respecto de las otras regiones, con IC que no se solapan.
De los contrastes, se sigue que en el Delta los odds aumentan un 245% respecto del Espinal y un 130% respecto de la Pampa. En la escala del log odds, las diferencias son otras pero los p valores son los mismos, pues es en esa escala que se realizan los tests.
Por otro lado, los odds no cambian significativamente entre Pampa y Espinal, o al menos no tenemos el suficiente poder estadístico para detectar esa diferencia.
Otro modo de visualizar lo que hizo nuestro modelo se ve en la figura 4:
tyto %<>%
mutate(
residuos_pearson = residuals(m5, type = "pearson"),
predichos_VR = predict(m5, type = "response"),
predichos_link = predict(m5, type = "link"),
cooks_D = cooks.distance(m5)
)
ggplot(tyto) +
aes(x = porc_urbano, y = proporcion, color = Region) +
geom_segment(aes(yend = predichos_VR, xend = porc_urbano), linetype = "dotted") +
geom_line(aes(y = predichos_VR), size = 1) +
geom_point(aes(y = predichos_VR), shape = 4, color = "Black") +
geom_point() +
labs(
title = "Fig. 4 La predicción es constante por región e ignora el porcentaje urbano"
)
Las líneas horizontales representan la probabilidad estimada de que las presas sean aves. Este valor es una constante que no depende del porc_urbano y es el mismo para todas las observaciones de cada región. Los valores estimados coinciden exactamente con el cociente entre el total de N_aves y el total de N_presas observados en cada región, como se ve en el siguiente chunk:
tyto %>%
group_by(Region) %>%
summarise(pi = sum(N_aves)/sum(N_presas)) %>%
mutate(pi = round(pi, 4))