Modelos de regresión

Fecha de publicación

12 de noviembre de 2023

Objetivos del manual

Familiarizarse con el formato para reportes dinámicos Rmarkdown/quarto
Ser capaz de documentar el manejo y análisis de datos en R usando Rmarkdown/quarto

1 Modelos lineales como marco unificador

Tradicionalmente, los modelos estadísticos se han enseñado como herramientas desconectadas sin una relación clara entre ellas. Sólo hace falta darle un vistazo a un libro clásico de estadística para biólogos (Sokal & Rolhf):

Sin embargo, la mayoría de esos modelos estadísticos comunes no son más que casos especiales de modelos lineales. Por lo tanto, aprenderlos como tales puede simplificar mucho las cosas. Este enfoque tiene varias ventajas:

En primer lugar, todo se reduce a y=a⋅x+b, lo que simplifica en gran medida el aprendizaje
Esto también significa que no hay necesidad de aprender sobre los parámetros, las hipótesis y la interpretación de los resultados para cada caso especial.
Los modelos lineales se han ampliado para tener en cuenta distribuciones y estructuras de datos complejas (por ejemplo, modelos mixtos, modelos lineales generalizados, modelos cero-inflados, etc.) proporcionando una plataforma más flexible.
Los modelos lineales se aplican en todos los paradigmas estadísticos (por ejemplo, frecuentista, bayesiano)

2 Objetivo del manual

- Comprender la inferencia estadística a través de una única herramienta de modelaje (modelos lineales en sentido amplio)

Familiarizarse con la construcción de modelos lineales
Extender los modelos lineales a diferentes estructuras de datos

Paquetes a utilizar en este manual:

Código

# instalar/cargar paquetes

sketchy::load_packages(
  c("ggplot2", 
    "viridis", 
    "lmerTest", 
    "sjPlot")
  )

3 Cómo simular datos

3.1 Generación de números aleatorios en R

La estadística nos permite inferir patrones en los datos. Solemos utilizar conjuntos de datos reales para enseñar estadística. Sin embargo, puede ser circular entender el funcionamiento interno de una herramienta estadística probando su capacidad para inferir un patrón que no estamos seguros de encontrar en los datos (y no tenemos idea del mecanismo que produjo ese patrón). Las simulaciones nos permiten crear escenarios controlados en los que conocemos con seguridad los patrones presentes en los datos y los procesos subyacentes que los han generado.

R ofrece algunas funciones básicas para la simulación de datos. Las más utilizadas son las funciones generadoras de números aleatorios. Los nombres de estas funciones comienzan con r (r____()). Por ejemplo, runif():

Código

# simular variable uniforme
unif_var <- runif(n = 100, min = 0, max = 10)

El resultado es un vector numérico de longitud 100 (n = 100):

Código

# imprimir variable
unif_var

  [1] 9.889093 3.977455 1.156978 0.697487 2.437494 7.920104 3.400624 9.720625
  [9] 1.658555 4.591037 1.717481 2.314771 7.728119 0.963015 4.534478 0.847007
 [17] 5.606659 0.087046 9.857371 3.165848 6.394489 2.952232 9.967037 9.060213
 [25] 9.887391 0.656457 6.270388 4.904750 9.710244 3.622208 6.799935 2.637199
 [33] 1.857143 1.851432 3.792967 8.470244 4.980761 7.905856 8.384639 4.569039
 [41] 7.994758 3.819431 7.597012 4.367756 9.042177 3.195349 0.825691 8.162891
 [49] 8.984762 9.664964 5.730689 7.200795 7.740586 6.277608 7.229893 3.868313
 [57] 1.627908 1.872283 3.912495 2.739012 1.919177 5.043918 7.638404 6.936689
 [65] 5.440542 6.590872 4.687284 4.818055 3.370636 4.245263 2.870151 6.011915
 [73] 8.407423 6.208370 1.345516 5.677224 4.434263 4.379754 6.236172 9.326533
 [81] 8.884926 8.785406 2.421769 7.414538 3.876563 0.789517 0.948356 7.621427
 [89] 3.478940 4.167667 3.440162 0.084109 9.115750 1.822054 7.228034 5.719633
 [97] 5.400364 3.549474 8.240918 1.861368

Podemos explorar el resultado graficando un histograma:

Código

# crear histograma
ggplot(data = data.frame(unif_var), mapping = aes(x = unif_var)) + geom_histogram()

Muestra una distribución uniforme que va de 0 a 10.

También podemos simular números aleatorios procedentes de una distribución normal utilizando rnorm():

Código

# crear una variable normal
norm_var <- rnorm(n = 1000, mean = 2, sd = 1)

# graficar histograma
ggplot(data = data.frame(norm_var), mapping = aes(x = norm_var)) + geom_histogram()

Tenga en cuenta que todas las funciones generadoras de números aleatorios tienen el argumento ‘n’, que determina la longitud del vector generado (es decir, el número de números aleatorios), además de algunos argumentos adicionales relacionados con parámetros específicos de la distribución.

Las variables continuas (es decir, los vectores numéricos) pueden convertirse en variables discretas (es decir, números enteros) simplemente redondeándolas:

Código

v1 <- rnorm(n = 5, mean = 10, sd = 3)

v1

[1] 15.4269  8.5450 15.1387  8.9203  9.3693

Código

round(x = v1, digits = 0)

[1] 15  9 15  9  9

Ejercicio 1

¿Qué hacen las funciones rbinom() y rexp()?
Ejecutela y haga histogramas de sus resultados
¿Qué hacen los argumentos ‘mean’ y ‘sd’ en rnorm()? Juege con diferentes valores y comprueba el histograma para hacerse una idea de su efecto en la simulación

3.2 Generación de variables categóricas

La forma más sencilla de generar variables categóricas es utilizar el vector de ejemplo letters' (oLETTERS’) para asignar niveles de categoría. Podemos hacerlo utilizando la función rep(). Por ejemplo, el siguiente código crea un vector categórico (caracteres) con dos niveles, cada uno con 4 observaciones:

Código

rep(letters[1:2], each = 4)

[1] "a" "a" "a" "a" "b" "b" "b" "b"

También podemos replicar este patrón utilizando el argumento ‘times’. Este código replica el vector anterior 2 veces:

Código

rep(letters[1:2], each = 4, times = 2)

 [1] "a" "a" "a" "a" "b" "b" "b" "b" "a" "a" "a" "a" "b" "b" "b" "b"

Otra opción es simular una variable a partir de una distribución binomial y luego convertirla en un factor:

Código

# correr rbinom
binom_var <- rbinom(n = 50, size = 1, prob = 0.5)

binom_var

 [1] 1 0 0 1 1 0 0 0 0 0 0 0 1 1 1 0 0 1 0 1 0 0 0 0 0 0 0 1 0 1 1 1 0 0 0 1 0 0
[39] 1 1 0 0 0 0 0 1 0 1 1 1

Código

# convertir a factor
categ_var <- factor(binom_var, labels = c("a", "b"))

categ_var

 [1] b a a b b a a a a a a a b b b a a b a b a a a a a a a b a b b b a a a b a a
[39] b b a a a a a b a b b b
Levels: a b

3.3 Muestreo aleatorio

La otra herramienta importante de R para jugar con datos simulados es sample(). Esta función permite tomar muestras de tamaños específicos de vectores. Por ejemplo, tomemos el ejemplo del vector ‘letters’:

Código

letters

 [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" "s"
[20] "t" "u" "v" "w" "x" "y" "z"

Podemos tomar una muestra de este vector como es:

Código

# tomar muestra
sample(x = letters, size = 10)

 [1] "i" "r" "q" "a" "c" "m" "y" "z" "u" "v"

El argumento ‘size’ nos permite determinar el tamaño de la muestra. Tenga en cuenta que obtendremos un error si el tamaño es mayor que el propio vector:

Código

sample(x = letters, size = 30)

Error in sample.int(length(x), size, replace, prob): cannot take a sample larger than the population when 'replace = FALSE'

Esto sólo puede hacerse cuando el muestreo es con reemplazo (replacement). El muestreo con reemplazo puede aplicarse estableciendo el argumento replace = TRUE:

Código

sample(x = letters, size = 30, replace = TRUE)

 [1] "j" "i" "h" "i" "k" "u" "n" "i" "w" "k" "c" "t" "m" "p" "i" "h" "d" "w" "c"
[20] "q" "j" "e" "e" "j" "k" "c" "l" "t" "r" "h"

3.4 Iterar un proceso

A menudo, las simulaciones deben repetirse varias veces para descartar resultados espurios debidos al azar o simplemente para probar diferentes parámetros. Las funciones de simulación de datos mencionadas anteriormente pueden ejecutarse varias veces (por ejemplo, iteradas) utilizando la función replicate():

Código

# replicar
repl_rnorm <- replicate(n = 3, expr = rnorm(2), simplify = FALSE)

# ver clase
class(repl_rnorm)

[1] "list"

Código

# imprimir
repl_rnorm

[[1]]
[1]  1.19184 -0.33992

[[2]]
[1]  0.78911 -0.63213

[[3]]
[1] -1.49312 -0.13441

3.5 Hacer que las simulaciones sean reproducibles

El último truco que necesitamos para ejecutar simulaciones en R es la capacidad de reproducir una simulación (es decir, obtener exactamente los mismos datos y resultados simulados). Esto puede ser útil para que otros investigadores puedan ejecutar nuestros análisis exactamente de la misma manera. Esto puede hacerse fácilmente con la función set.seed(). Pruebe a ejecutar el siguiente código. Debería obtener la misma salida:

Código

# definir semilla
set.seed(10)

# crear variable uniforme
runif(n = 2)

[1] 0.50748 0.30677

4 Crear juegos de datos

4.1 Juegos de datos con variables numéricas y categóricas

Ahora que sabemos cómo simular variables continuas y categóricas. Podemos juntarlas para crear conjuntos de datos simulados. Esto se puede hacer utilizando la función data.frame():

Código

# crear variable categorica
grupo <- rep(letters[1:2], each = 3)

# crear variable continuaa
size <- rnorm(n = 6, mean = 5, sd = 1)

# poner juntas en un data frame
df <- data.frame(grupo, size)

# imprimir
df

grupo	size
a	4.8158
a	3.6287
a	4.4008
b	5.2946
b	5.3898
b	3.7919

Por supuesto, podríamos añadir más variables a este cuadro de datos:

Código

# crear variable categorica
grupo <- rep(letters[1:2], each = 3)
individuo <- rep(LETTERS[1:6])

# crear variables continuas
size <- rnorm(n = 6, mean = 5, sd = 1)
weight <- rnorm(n = 6, mean = 100, sd = 10)

# poner todo en un data frame
df <- data.frame(grupo, individuo, size, weight)

# imprimir
df

grupo	individuo	size	weight
a	A	4.6363	109.874
a	B	3.3733	107.414
a	C	4.7435	100.893
b	D	6.1018	90.451
b	E	5.7558	98.049
b	F	4.7618	109.255

Y eso es un juego de datos simulados en su forma más básica. Se parece mucho al tipo de datos con los que trabajamos en biología.

5 Cómo utilizar datos simulados para entender el comportamiento de las herramientas estadísticas

5.1 Prueba de concepto: el Teorema Central del Límite

El Teorema del Límite Central afirma que, si tomamos muestras aleatorias de una población, los promedios de esas muestras seguirán una distribución normal, aunque la población no esté distribuida normalmente. Además, la distribución normal resultante debe tener un promedio cercano al promedio de la población. El teorema es un concepto clave para la estadística inferencial, ya que implica que los métodos estadísticos que funcionan para las distribuciones normales pueden ser aplicables a muchos problemas que implican otros tipos de distribuciones. No obstante, el objetivo aquí es sólo mostrar cómo se pueden utilizar las simulaciones para entender el comportamiento de los métodos estadísticos.

Para comprobar si esas afirmaciones básicas sobre el Teorema del Límite Central son ciertas, podemos utilizar datos simulados en R. Vamos a simular una población de 1000 observaciones con una distribución uniforme:

Código

# simular popublacion uniforme
unif_pop <- runif(1000, min = 0, max = 10)

# ver histograma
ggplot(data = data.frame(unif_pop), mapping = aes(x = unif_pop)) + geom_histogram()

Podemos tomar muestras aleatorias usando sample() así:

Código

sample(x = unif_pop, size = 30)

 [1] 9.28420 1.02626 2.57517 3.32485 6.89990 2.29404 0.33737 8.21366 3.30364
[10] 8.03793 2.59174 7.81770 5.65426 0.63831 2.83470 4.20434 4.76330 4.42193
[19] 6.97830 7.92625 0.68121 3.52323 6.51103 5.38289 7.97210 1.80062 4.21282
[28] 3.33866 8.91223 4.71163

Este proceso puede ser replicado varias veces con replicate():

Código

# replicar
samples <- replicate(n = 100, expr = mean(sample(x = unif_pop, size = 30)))

El código anterior toma 100 muestras con 30 valores cada una. Ahora podemos comprobar la distribución de las muestras:

Código

# ver distribucion/ histograma
ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram()

… asi como el promedio:

Código

mean(samples)

[1] 5.0212

Como era de esperar, las muestras siguen una distribución normal con una media cercana a la media de la población, que es:

Código

mean(unif_pop)

[1] 5.0527

Probemos con una distribución más compleja. Por ejemplo, una distribución bimodal:

Código

# usar semilla
set.seed(123)

# simular variables
norm1 <- rnorm(n = 1000, mean = 10, sd = 3)
norm2 <- rnorm(n = 1000, mean = 20, sd = 3)

# juntar en una sola variable
bimod_pop <- c(norm1, norm2)

# ver histograma
ggplot(data = data.frame(bimod_pop), mapping = aes(x = bimod_pop)) + geom_histogram()

Código

# replicar muestreo
samples <- replicate(200, mean(sample(bimod_pop, 10)))

# ver histograma
ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram()

Código

# ver promedios
mean(samples)

[1] 15.231

Código

mean(bimod_pop)

[1] 15.088

Ejercicio 2

Intenta explorar el Teorema del Límite Central como en el caso anterior, pero esta vez utilizando:
1. Una distribución exponencial (rexp())
2. Una distribución binomial (rbinom())

Para cada distribución: grafique un histograma y compare los promedios de la población y de las muestras

Regresión lineal simple

Las regresiones lineales se basan en la ecuación lineal a = mx + b que aprendimos en el colegio. La representación formal tiene este aspecto:

$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1}$

$\hat{Y}$: variable respuesta
$\beta_{o}$: intercepto
$\beta_{1}$: estimado de la magnitud del efecto de $x_{1}$ en $\hat{Y}$ (también conocido como tamaño de efecto, coeficiente o simplemente “estimado”)
$x_{1}$: variable predictora

El objetivo más común de una regresión lineal es la estimación de los valores de $\beta_{*}$. Esto se consigue encontrando la línea recta que mejor se ajusta y que representa la asociación entre un predictor y la respuesta:

Estos $\beta_{*}$ son el tamaño de efecto estimado del predictor correspondiente (por ejemplo, $\beta_{1}$ es el tamaño de efecto $x_{1}$). Su valor representa el cambio medio en $\hat{Y}$ (en unidades $\hat{Y}$) para una unidad de cambio en $\beta_{*}$. Por lo tanto, la hipótesis nula es que esos $\beta_{*}$ no son diferentes de 0:

$\qquad \mathcal{H}_0: \hat{Y} = \beta_0 + 0 * x_{1}$

lo que equivale a esto:

$\qquad \mathcal{H}_0: \hat{Y} = \beta_0$

Por ejemplo, un modelo de regresión con este resultado:

6 Regresión lineal en R

Para sacar el máximo provecho de los modelos lineales necesitamos sentirnos cómodos con ellos. Empezaremos explorando la función de regresión lineal de R lm(). En R la mayoría de los modelos lineales y sus extensiones comparten formatos comunes de entrada de datos y resultados, lo que facilita su aplicación una vez que entendemos sus fundamentos.

Utilizaremos el juego de datos ‘trees’ que viene por defecto con R. ‘trees’ proporciona medidas del diámetro (etiquetado como ‘Girth’), altura y volumen de 31 cerezos talados:

Código

head(trees)

Girth	Height	Volume
8.3	70	10.3
8.6	65	10.3
8.8	63	10.2
10.5	72	16.4
10.7	81	18.8
10.8	83	19.7

La función básica de R para construir un modelo lineal es lm(). Veamos los componentes básicos de un modelo de regresión utilizando lm():

Podemos correr este modelo para ver el resultado:

Código

reg_mod <- lm(formula = Height ~ Girth, data = trees)

summary(reg_mod)


Call:
lm(formula = Height ~ Girth, data = trees)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.582  -2.769   0.316   2.473   9.946 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   62.031      4.383   14.15  1.5e-14 ***
Girth          1.054      0.322    3.27   0.0028 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.54 on 29 degrees of freedom
Multiple R-squared:  0.27,  Adjusted R-squared:  0.244 
F-statistic: 10.7 on 1 and 29 DF,  p-value: 0.00276

Esto es lo que significan los elementos de la salida:

Llamada (call): la función y los parámetros que se utilizaron para crear el modelo
Residuales (residuals): distribución de los residuales. Los residuales son la diferencia entre lo que predijo el modelo y el valor real de y. Esta es una representación gráfica de los residuales:

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

Coeficientes (coefficients): contiene los tamaños de efecto (‘Estimates’), una medida de su incertidumbre (‘Standar error’), la estadística asociada (‘t value’) y el valor p (‘Pr(>|t|)’). Los tamaño de efecto o estimados se dan como el cambio promedio en y por cada aumento de 1 unidad en x. Para este ejemplo el estimado es de 1,0544 $in / in$. A pesar de que las unidades se anulan (`1,0544 $in / in$ = 1,0544) tenerlas en cuenta sigue siendo biológicamente significativo. Significan que, en promedio, y un aumento de 1 pulgada en la circunferencia que se espera un aumento de 1,0544 pulgadas en la altura.
Error estándar de los residuales (residual standard error): se explica por sí mismo. El error estándar de los residuales
R-cuadrado múltiple (multiple R-squared): el coeficiente de determinación, que pretende ser una medida de lo bien que su modelo se ajusta a los datos
R-cuadrado ajustado (adjusted R-squared): similar al ‘R-cuadrado múltiple’ pero penalizado por el número de parámetros
Estadístico F (F-statistic): estadístico para una prueba global que comprueba si al menos uno de sus coeficientes es distinto de cero
Valor de p: probabilidad de una prueba global que comprueba si al menos uno de sus coeficientes es distinto de cero

Utilizaremos lm() para mostrar la flexibilidad de los modelos de regresión. Los componentes de regresión se añadirán gradualmente para que podamos tomarnos el tiempo de entender cada uno de ellos así como los correspondientes cambios en la salida de la regresión.

6.1 Modelo solo con intercepto

Primero vamos a crear una variable numérica de respuesta:

Código

# definir semilla
set.seed(123)

# numero de observaciones
n <- 50

#  variables aleatorias
y <- rnorm(n = n, mean = 0, sd = 1)

# put it in a data frame
y_data <- data.frame(y)

Esta única variable puede introducirse en un modelo de regresión sólo con intercepto. Para ello, debemos suministrar la fórmula del modelo y los datos a lm():

Código

# run model
y_mod <- lm(formula = y ~ 1, data = y_data)

Lo que equivale a:

$\hat{Y} \sim \beta_{o}$

Podemos obtener el resumen por defecto de los resultados del modelo ejecutando summary() en el objeto de salida ‘y_mod’:

Código

summary(y_mod)


Call:
lm(formula = y ~ 1, data = y_data)

Residuals:
   Min     1Q Median     3Q    Max 
-2.001 -0.594 -0.107  0.664  2.135 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.0344     0.1309    0.26     0.79

Residual standard error: 0.926 on 49 degrees of freedom

Puede ser bastante informativo graficar los tamaños de efecto (aunque en este caso sólo tenemos uno):

Código

ci_df <- data.frame(param = names(y_mod$coefficients), 
                    est = y_mod$coefficients, confint(y_mod))

ggplot(ci_df, aes(x=param, y=est)) + 
  geom_hline(yintercept = 0, color="red", lty = 2) +
  geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + 
  labs(x = "Parámetro", y = "Tamaño de efecto") + 
  coord_flip()

Interpretación del modelo

Para evaluar la importancia de la asociación nos centramos en la tabla de coeficientes:

            Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.034404    0.13094 0.26275  0.79385

En este ejemplo no hay predictores en el modelo, por lo que sólo tenemos una estimación para el intercepto ($\beta_0$)
El modelo nos dice que el intercepto se estima en 0.0344 y que este valor no es significativamente diferente de 0 (valor p = 0.79385)
En este caso el intercepto es simplemente la media de la variable de respuesta

Código

mean(y_data$y)

[1] 0.034404

Rara vez tenemos predicciones sobre el intercepto, por lo que tendemos a ignorar esta estimación.

Caso de estudio

Ruhs, E. C., Martin, L. B., & Downs, C. J. (2020). The impacts of body mass on immune cell concentrations in birds. Proceedings of the Royal Society B, 287(1934), 20200655.

“Encontramos que un modelo sólo con intercepto (intercept-only model) explicaba mejor las concentraciones de linfocitos y eosinófilos en las aves, indicando que las concentraciones de estos tipos de células eran independientes de la masa corporal.”

Ejercicio 3

Cambie el argumento ‘mean’ en la llamada de la función “rnorm()` (línea 130) por un valor distinto de 0 y observe cómo cambian los valores en la tabla de coeficientes
Cambia el argumento sd en la llamada a la función rnorm() (línea 130) por un valor más alto y observe cómo cambian los valores en la tabla de coeficientes

6.2 Añadir un predictor no asociado

Podemos crear 2 variables numéricas no relacionadas así:

Código

# definir semilla
set.seed(123)

# numero de observaciones
n <- 50

#  variables aleatorias
y <-  rnorm(n = n, mean = 0, sd = 1)
x1 <-  rnorm(n = n, mean = 0, sd = 1)

# crear data frame
xy_datos <- data.frame(x1, y)

Estas dos variables pueden introducirse en un modelo de regresión para evaluar la asociación entre ellas:

Código

# construir model
xy_mod <- lm(formula = y ~ x1, data = xy_datos)

# graficar
ggplot(xy_datos, aes(x = x1, y = y)) + 
  geom_smooth(method = "lm", se = FALSE) + 
  geom_point() # graficar points

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

Que es equivalente a esto:

$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1}$

Imprimamos el resumen de este modelo:

Código

summary(xy_mod)


Call:
lm(formula = y ~ x1, data = xy_datos)

Residuals:
   Min     1Q Median     3Q    Max 
-2.004 -0.624 -0.123  0.687  2.106 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.0398     0.1340    0.30     0.77
x1           -0.0367     0.1475   -0.25     0.80

Residual standard error: 0.935 on 48 degrees of freedom
Multiple R-squared:  0.00129,   Adjusted R-squared:  -0.0195 
F-statistic: 0.0618 on 1 and 48 DF,  p-value: 0.805

… y graficar los tamaños de efecto:

Código

ci_df <- data.frame(param = names(xy_mod$coefficients), 
                    est = xy_mod$coefficients, confint(xy_mod))

ggplot(ci_df, aes(x=param, y=est)) + 
  geom_hline(yintercept = 0, color="red", lty = 2) +
  geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + 
  labs(x = "Parámetro", y = "Tamaño de efecto") + 
  coord_flip()

Deberíamos “diagnosticar” la idoneidad del modelo inspeccionando más de cerca la distribución de los residuales. La función plot_model() del paquete sjPlot hace un buen trabajo para crear gráficos de diagnóstico para modelos lineales:

Código

plot_model(xy_mod, type = "diag")

[[1]]

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


[[2]]


[[3]]

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

Interpretación del modelo

Cuadro con coeficientes:

             Estimate Std. Error  t value Pr(>|t|)
(Intercept)  0.039774    0.13396  0.29690  0.76782
x1          -0.036679    0.14750 -0.24867  0.80467

En este ejemplo hemos añadido un predictor al modelo, por lo que hemos obtenido una estimado adicional (y una fila extra, ‘x1’)
El modelo nos dice que la estimación de ‘x1’ es -0.03668 y que no es significativamente diferente de 0 (p = 0.80467)

6.3 Simular un predictor asociado

Podemos utilizar la fórmula del modelo lineal anterior para simular dos variables continuas asociadas así:

Código

# definir semilla
set.seed(123)

# numero de observaciones
n <- 50
b0 <- -4
b1 <- 3
error <- rnorm(n = n, sd = 3)

# variables aleatorias
x1 <-  rnorm(n = n, mean = 0, sd = 1)
y <- b0 + b1 * x1 + error

# crear data frame
xy_datos2 <- data.frame(x1, y)

Note que también hemos añadido un término de error, por lo que la asociación no es perfecta. Vamos a correr el modelo y graficar la asociación entre las dos variables:

Código

# construir model
xy_mod2 <- lm(formula = y ~ x1, data = xy_datos2)

# graficar
ggplot(xy_datos2, aes(x = x1, y = y)) + 
  geom_smooth(method = "lm", se = FALSE) +  
  geom_point() # graficar points

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

La fórmula es la misma que la del modelo anterior:

$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1}$

Este es el resumen del modelo:

Código

summary(xy_mod2)


Call:
lm(formula = y ~ x1, data = xy_datos2)

Residuals:
   Min     1Q Median     3Q    Max 
 -6.01  -1.87  -0.37   2.06   6.32 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -3.881      0.402   -9.66  7.9e-13 ***
x1             2.890      0.442    6.53  3.9e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.8 on 48 degrees of freedom
Multiple R-squared:  0.471, Adjusted R-squared:  0.459 
F-statistic: 42.7 on 1 and 48 DF,  p-value: 3.85e-08

.. el gráfico con los tamaños de efecto:

Código

ci_df <- data.frame(param = names(xy_mod2$coefficients), 
                    est = xy_mod2$coefficients, confint(xy_mod2))

ggplot(ci_df, aes(x=param, y=est)) + 
  geom_hline(yintercept = 0, color="red", lty = 2) +
  geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + 
  labs(x = "Parámetro", y = "Tamaño de efecto") + 
  coord_flip()

… y los gráficos diagnósticos del modelo:

Código

plot_model(xy_mod2, type = "diag")

[[1]]

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


[[2]]


[[3]]

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

Interpretación del modelo

Cuadro con los coeficientes (estimados):

            Estimate Std. Error t value   Pr(>|t|)
(Intercept)  -3.8807    0.40188 -9.6562 7.8616e-13
x1            2.8900    0.44249  6.5311 3.8537e-08

El modelo nos dice que $beta_1$ (el tamaño de efecto de ‘x1’) es 2.88996 y que es significativamente diferente de 0 (p = 3.85365^{-8})
Los valores simulados de los parámetros de regresión pueden compararse con el resumen del modelo lm() para tener una idea de la precisión del modelo:
- $\beta_1$ (the effect size of ‘x1’) was set to 3 and was estimated as 2.89 by the model

Caso de estudio

Keenan EL, Odom KJ, Araya-Salas M, Horton KG, Strimas-Mackey M, Meatte MA, Mann NI, Slater PJ, Price JJ, and Templeton CN. 2020. Breeding season length predicts duet coordination and consistency in Neotropical wrens (Troglodytidae). Proceeding of the Royal Society B. 20202482

“… la coordinación y la consistencia de los duetos son mayores en las especies con temporadas de apareamiento especialmente largas.”

Ejercicio 4

Aumente el tamaño de la muestra (n) a 1000 o más
¿Cómo cambiaron los estimados del tamaño de efecto ($\beta$)?
¿Cómo cambió el error estándar del tamaño de efecto?
Ahora cambie n a 15 y compruebe de nuevo las estimaciones del modelo (esta vez compruebe también el valor p)

6.4 Varios predictores: regresión múltiple

La regresión lineal múltiple es una extensión del modelo de regresión lineal simple que puede tomar varios predictores:

$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1} + \cdots + \beta_{n} * x_{n}$

La fórmula parece compleja, pero sólo quiere decir que cualquier parámetro adicional tendrá su propia estimación ($\beta$). La fórmula para una regresión lineal de dos predictores se ve así:

$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1} + \beta_{2} * x_{2}$

… y se puede simular así:

Código

# semila
set.seed(123)

# numero de observaciones
n <- 50
b0 <- -4
b1 <- 3
b2 <- -2
error <- rnorm(n = n, mean = 0, sd = 3)

# variables aleatorias
x1 <-  rnorm(n = n, mean = 0, sd = 1)
x2 <-  rnorm(n = n, mean = 0, sd = 1)
y <- b0 + b1 * x1 + b2 * x2 + error

# crear un data frame
xy_datos_multp <- data.frame(x1, x2, y)

# construir el modelo
xy_mod_multp <- lm(formula = y ~ x1 + x2, data = xy_datos_multp)

summary(xy_mod_multp)


Call:
lm(formula = y ~ x1 + x2, data = xy_datos_multp)

Residuals:
   Min     1Q Median     3Q    Max 
-5.986 -1.893 -0.363  2.002  6.413 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -3.865      0.417   -9.27  3.5e-12 ***
x1             2.901      0.453    6.41  6.4e-08 ***
x2            -1.932      0.414   -4.67  2.6e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.83 on 47 degrees of freedom
Multiple R-squared:  0.612, Adjusted R-squared:  0.595 
F-statistic:   37 on 2 and 47 DF,  p-value: 2.23e-10

… graficar los tamaños de efecto:

Código

ci_df <- data.frame(param = names(xy_mod_multp$coefficients), 
                    est = xy_mod_multp$coefficients, confint(xy_mod_multp))

ggplot(ci_df, aes(x=param, y=est)) + 
  geom_hline(yintercept = 0, color="red", lty = 2) +
  geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + 
  labs(x = "Parámetro", y = "Tamaño de efecto") + 
  coord_flip()

… y los gráficos diagnósticos:

Código

plot_model(xy_mod_multp, type = "diag")

[[1]]

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


[[2]]


[[3]]

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

Interpretación del modelo

Cuadro con los coeficientes (estimados):

            Estimate Std. Error t value   Pr(>|t|)
(Intercept)  -3.8652    0.41695 -9.2702 3.4856e-12
x1            2.9015    0.45260  6.4108 6.4137e-08
x2           -1.9324    0.41422 -4.6650 2.5851e-05

El modelo encontró que $\beta_1$ (el tamaño de efecto de ‘x1’) es 2.9015 y que es significativamente diferente de 0 (p = 6.41368^{-8}). Este es el efecto de ‘x1’ sobre ‘y’ una vez removida la variación explicada por ‘x2’.
También se encontró que el $\beta_2$ (el tamaño de efecto de ‘x2’) es -1.93235 y que también es significativamente diferente de 0 (p = 2.58513^{-5}). Este es el efecto de ‘x2’ sobre ‘y’ una vez removida la variación explicada por ‘x1’.
Los valores simulados de los parámetros de regresión pueden compararse con el resumen del modelo lm() para tener una idea de la precisión del modelo:
- $\beta_1$ se fijó en 3 y se estimó como 2.901.
- $\beta_2$ (el tamaño de efecto de ‘x2’) se fijó en -2 y se estimó como -1.932.

Caso de estudio

Araya-Salas M, P González-Gómez, K Wojczulanis-Jakubas, V López III & T Wright. 2018. Spatial memory is as important as weapon and body size for territorial ownership in a lekking hummingbird. Scientific Reports. 13, e0189969

“La memoria espacial, el tamaño corporal y el largo de la punta del pico … predijeron positivamente la probabilidad de adquirir y defender un territorio.”

Hay un punto importante que es necesario enfatizar aquí: la regresión múltiple estima el efecto de un predictor después de tener en cuenta el efecto de los demás predictores del modelo. En otras palabras, los nuevos predictores del modelo tratarán de explicar la variación de los datos que no fue explicada por los otros predictores. Así que el resultado de la regresión múltiple no es equivalente a los resultados de las regresiones lineales simples sobre los mismos predictores. Esto puede demostrarse fácilmente corriendo esas regresiones:

Código

# construir modelos
x1y_mod <- lm(formula = y ~ x1, data = xy_datos)
x2y_mod <- lm(formula = y ~ x2, data = xy_datos)

# atajo a los coeficientes
coef(xy_mod)

(Intercept)          x1 
   0.039774   -0.036679

Código

coef(x1y_mod)

(Intercept)          x1 
   0.039774   -0.036679

Código

coef(x2y_mod)

(Intercept)          x2 
   0.041320    0.027239

Las estimaciones de las mismas variables varían considerablemente entre la regresión múltiple y las regresiones de un solo predictor.

Este punto se demuestra además por el hecho de que, si uno de los predictores no tiene ninguna influencia en la respuesta, el efecto del predictor adicional convergerá a su efecto en una regresión lineal simple. Para simular este escenario, fijamos b2 en 0:

Código

# definir semilla
set.seed(123)

# numero de observaciones
n <- 50
b0 <- -4
b1 <- 3
b2 <- 0
error <- rnorm(n = n, mean = 0, sd = 1)

#  variables aleatorias
x1 <-  rnorm(n = n, mean = 0, sd = 1)
x2 <-  rnorm(n = n, mean = 0, sd = 1)
y <- b0 + b1 * x1 + b2 * x2 + error

# crear data frame
xy_datos <- data.frame(x1, x2, y)

# construir modelos
xy_mod <- lm(formula = y ~ x1 + x2, data = xy_datos)
x1y_mod <- lm(formula = y ~ x1, data = xy_datos)

# shortcut to coefficients
coef(xy_mod)

(Intercept)          x1          x2 
  -3.955064    2.967166    0.022549

Código

coef(x1y_mod)

(Intercept)          x1 
    -3.9602      2.9633

El estimado de $\beta_1$ fue casi el mismo en la regresión múltiple (2.96717) y en la regresión de un solo predictor (2.96332)

Por comodidad, hemos utilizado coef() para extraer sólo los estimado de la regresión, pero los valores son los mismos que obtenemos con summary(model).

Ejercicio 5

La siguiente simulación genera un par de predictores continuos altamente colineales:

Código

# definir semilla
set.seed(123)

# numero de observaciones
n <- 30
b0 <- -4
b1 <- 3
b2 <- -2
error <- rnorm(n = n, mean = 0, sd = 1)

#  variables aleatorias colineales
x1 <-  rnorm(n = n, mean = 0, sd = 1)
x2 <-  x1 + rnorm(n = n, mean = 0, sd = 0.2) # hacer x2 colineal con x1
y <- b0 + b1 * x1 + b2 * x2 + error

# crear data frame
xy_datos <- data.frame(x1, x2, y)

cor(x1, x2)

[1] 0.97812

Código

# graficar
ggplot(data = xy_datos, aes(x = x1, y = x2)) +
    geom_point(size = 3) +

Construya un modelo lineal multiple (con lm()) con ‘y’ como respuesta y ‘x1’ y ‘x2’ como predictores
Compare los betas estimados por el modelo con los usados para generar los datos. ¿Es una buena estimación?. También observe los valores de p. ¿Esperaría que fueran significativos?

6.5 Añadir un predictor categórico

Para los predictores categóricos podemos crear primero una variable binaria (0, 1) y luego añadir etiquetas a cada valor:

Código

# definir semilla
set.seed(13)

# numero de observaciones
n <- 50
b0 <- -3
b1 <- 2
error <- rnorm(n = n, mean = 0, sd = 3)

#  variables aleatorias
x1_num <- sample(0:1, size = n, replace = TRUE)
y <- b0 + b1 * x1_num + error

x1 <- factor(x1_num, labels = c("a", "b"))

# crear data frame
xy_datos_cat <- data.frame(x1, x1_num, y)

head(xy_datos_cat)

x1	x1_num	y
b	1	0.66298
a	0	-3.84082
a	0	2.32549
b	1	-0.43804
a	0	0.42758
a	0	-1.75342

Y así es como se escribe formalmente:

$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1}$

Lo mismo que con los predictores continuos.

Podemos explorar el patrón de los datos utilizando un diagrama de cajas (boxplot):

Código

# graficar
ggplot(xy_datos_cat, aes(x = x1, y = y)) + 
  geom_boxplot()

… y obtener los estimados del modelo:

Código

# construir modelos
xy_mod_cat <- lm(formula = y ~ x1, data = xy_datos_cat)

summary(xy_mod_cat)


Call:
lm(formula = y ~ x1, data = xy_datos_cat)

Residuals:
   Min     1Q Median     3Q    Max 
-5.898 -1.909 -0.094  1.809  5.506 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -2.997      0.558   -5.37  2.3e-06 ***
x1b            1.814      0.842    2.16    0.036 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.95 on 48 degrees of freedom
Multiple R-squared:  0.0882,    Adjusted R-squared:  0.0693 
F-statistic: 4.65 on 1 and 48 DF,  p-value: 0.0362

… graficar los tamaños de efecto:

Código

ci_df <- data.frame(param = names(xy_mod_cat$coefficients), 
                    est = xy_mod_cat$coefficients, confint(xy_mod_cat))

ggplot(ci_df, aes(x=param, y=est)) + 
  geom_hline(yintercept = 0, color="red", lty = 2) +
  geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + 
  labs(x = "Parámetro", y = "Tamaño de efecto") + 
  coord_flip()

… y los gráficos diagnósticos del modelo:

Código

plot_model(xy_mod_cat, type = "diag")[[2]]

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

Interpretación del modelo

Cuadro con los coeficientes (estimados):

            Estimate Std. Error t value   Pr(>|t|)
(Intercept)  -2.9974    0.55825 -5.3693 2.2677e-06
x1b           1.8140    0.84160  2.1554 3.6172e-02

El modelo encontró que $\beta_1$ (el tamaño de efecto de ‘x1’) es 1.814 y que es significativamente diferente de 0 (p = 0.03617)
Los valores simulados de los parámetros de regresión pueden compararse con el resumen del modelo lm() para tener una idea de la precisión del modelo:
- $\beta_1$ se fijó en 2 y se estimó como 1.814.
Tenga en cuenta que en este caso el intercepto se refiere a la estimación del nivel ‘a’ en el predictor categórico, que se utilizó como base:

Código

# graficar
ggplot(xy_datos_cat, aes(x = x1, y = y)) + 
  geom_boxplot() +
geom_hline(yintercept = xy_mod_cat$coefficients[1], col = "blue")

Por lo tanto, el intercepto es el mismo que la media de y para la categoría ‘a’:

Código

mean(xy_datos_cat$y[xy_datos_cat$x1 == "a"])

[1] -2.9974

Observe también que la etiqueta del estimado es ‘x1b’, no ‘x1’ como en los predictores continuos. Esto se debe a que en este caso el estimado se refiere a la diferencia entre los dos niveles de la variable categórica (‘a’ y ‘b’). Más concretamente, nos dice que en promedio las observaciones de la categoría ‘b’ son 1.814 más altas que las observaciones de la categoría ‘a’.

Caso de estudio

Rico-Guevara A, & M Araya-Salas. 2015. Bills as daggers? A test for sexually dimorphic weapons in a lekking hummingbird species. Behavioral Ecology. 26 (1): 21-29.

“Los machos con puntas de pico más grandes y puntiagudas tuvieron más éxito en ganar control de territorios en el lek.”

Ejercicio 6

Los datos desbalanceados cuando hay categorías (es decir, algunas categorías tienen muchas más observaciones que otras) pueden ser problemáticos para la inferencia estadística. Modifique el código arriba para simular un juego de datos muy desbalanceado y compruebe la precisión del modelo.

Variables indicadoras (Dummy variables)

En un modelo de regresión, los predictores categóricos también se representan como vectores numéricos. Más concretamente, los predictores categóricos se codifican como 0s y 1s, en los que 1 significa “pertenece a la misma categoría” y 0 “pertenece a una categoría diferente”. Mantuvimos el vector numérico original (‘x1_num’) al simular el juego de datos con el predictor categórico:

Código

head(xy_datos_cat)

x1	x1_num	y
b	1	0.66298
a	0	-3.84082
a	0	2.32549
b	1	-0.43804
a	0	0.42758
a	0	-1.75342

Observe que las “b” de la columna “x1” se convierten en 1 en la columna “x1_num” y las “a” se convierten en 0. Esto se denomina variable indicadora y el proceso se conoce como ‘codificación indicadora’ (dummy coding).

En realidad, podemos utilizar el vector numérico en el modelo de regresión y obtener exactamente los mismos resultados:

Código

# summary model with categorical variable
summary(xy_mod_cat)$coefficients

            Estimate Std. Error t value   Pr(>|t|)
(Intercept)  -2.9974    0.55825 -5.3693 2.2677e-06
x1b           1.8140    0.84160  2.1554 3.6172e-02

Código

# construir modelos with dummy variable
xy_mod_num <- lm(formula = y ~ x1_num, data = xy_datos_cat)

# summary with dummy coding
summary(xy_mod_num)$coefficients

            Estimate Std. Error t value   Pr(>|t|)
(Intercept)  -2.9974    0.55825 -5.3693 2.2677e-06
x1_num        1.8140    0.84160  2.1554 3.6172e-02

Las cosas se complican un poco más cuando se codifica un predictor categórico con más de dos niveles. Pero la lógica es la misma.

6.6 Interacciones

Una interacción estadística se refiere a un efecto de una variable predictora que está mediado por una segunda variable.

$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1} + \beta_{2} * x_{2} + \beta_{3} * x_{1} * x_{2}$

Esto es más fácil de entender si se observa la interacción de una variable continua y una binaria:

Código

# definir semilla
set.seed(123)

# numero de observaciones
n <- 50
b0 <- -4
b1 <- 3
b2 <- 1.7
b3 <- -3
error <- rnorm(n = n, mean = 0, sd = 3)

#  variables aleatorias
x1 <- rbinom(n = n, size = 1, prob = 0.5)
x2 <-  rnorm(n = n, mean = 0, sd = 1)

# interaccion se añade como el producto dex1 y x2
y <- b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2 + error

x1 <- factor(x1, labels = c("a", "b"))

# crear data frame
xy_datos_intr <- data.frame(x1, x2, y)

head(xy_datos_intr)

x1	x2	y
b	1.02557	-4.0147
a	-0.28477	-5.1746
a	-1.22072	-1.3991
b	0.18130	-1.0242
a	-0.13889	-3.8483
b	0.00576	4.1377

Código

# construir modelos
xy_mod_intr <- lm(formula = y ~ x1 + x2 + x1 * x2, data = xy_datos_intr)

# guardar resumen para graficar lineas de mejor ajuste
xy_summ_intr <- summary(xy_mod_intr)

xy_summ_intr


Call:
lm(formula = y ~ x1 + x2 + x1 * x2, data = xy_datos_intr)

Residuals:
   Min     1Q Median     3Q    Max 
-6.222 -1.664 -0.158  1.650  6.383 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -4.193      0.570   -7.35  2.7e-09 ***
x1b            3.597      0.806    4.46  5.2e-05 ***
x2             1.327      0.682    1.95   0.0576 .  
x1b:x2        -2.972      0.962   -3.09   0.0034 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.83 on 46 degrees of freedom
Multiple R-squared:  0.396, Adjusted R-squared:  0.357 
F-statistic:   10 on 3 and 46 DF,  p-value: 3.29e-05

También ayuda graficar los datos:

Código

# graficar
ggplot(data = xy_datos_intr, aes(x = x2, y = y, color = x1)) +
    geom_point(size = 3) +
    geom_smooth(method = "lm", se = FALSE)

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

… y los tamaños de efecto:

Código

ci_df <- data.frame(param = names(xy_mod_intr$coefficients), 
                    est = xy_mod_intr$coefficients, confint(xy_mod_intr))

ggplot(ci_df, aes(x=param, y=est)) + 
  geom_hline(yintercept = 0, color="red", lty = 2) +
  geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + 
  labs(x = "Parámetro", y = "Tamaño de efecto") + 
  coord_flip()

También deberíamos revisar los gráficos de diagnóstico:

Código

plot_model(xy_mod_intr, type = "diag")

[[1]]


[[2]]

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


[[3]]


[[4]]

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

Interpretación del modelo

Cuadro con los coeficientes:

            Estimate Std. Error t value   Pr(>|t|)
(Intercept)  -4.1930    0.57023 -7.3532 2.6978e-09
x1b           3.5974    0.80607  4.4629 5.1950e-05
x2            1.3274    0.68161  1.9474 5.7610e-02
x1b:x2       -2.9720    0.96227 -3.0885 3.4061e-03

El modelo encontró que $\beta_1$ (el tamaño de efecto de ‘x1-b’ a ‘x1-a’) es 3.59741 y que es significativamente diferente de 0 (valor p = 5.19496^{-5})
El modelo encontró que $\beta_2$ (el tamaño de efecto de ‘x2’) es 1.32735 y que es significativamente diferente de 0 (p = 0.05761). Esto es en realidad la pendiente de la relación entre x2 e y cuando x1 = `a’
El modelo encontró que $\beta_3$ (el tamaño de efecto del término de interacción ‘x1 * x2’) es -2.97196 y que es significativamente diferente de 0 (p = 0.00341). Esta es la diferencia entre las pendientes de x2 vs y cuando x1 = ‘a’ y x2 vs y cuando x1 = ‘b’.
Los valores simulados para los parámetros de regresión pueden compararse con el resumen del modelo lm() para tener una idea de la precisión del modelo:
- $\beta_1$ se simuló con un valor de 3 y se estimó en 3.597
- $\beta_2$ se simuló con un valor de 1.7 y se estimó en 1.327
- $\beta_3$ se simuló con un valor de -3 y se estimó en -2.972

Caso de estudio

Chirino F, B Wilink, Araya-Salas M. In prep. Climatic factors affecting vocal activity in lemur leaf frogs.

“El aumento de la luz de la luna disminuye la actividad vocal de Agalychnis lemur aunque esta relación es mediada por la temperatura.”

Ejercicio 7

Modifique el código usado para simular un único predictor asociado aumentando gradualmente el error. Esto se hace aumentando el argumento ‘sd’ en error <- rnorm(n = n, sd = 2).
Mire cómo los errores más grandes afectan la inferencia (debe volver a correr los modelos)

Ahora reemplace el término de error con error <- rexp(n = n, rate = 0.2). Esto crea un error con una distribución exponencial (por lo tanto no normal). Se supone que esto es problemático para el poder inferencial de estos modelos. Compare las estimaciones que obtuvo con los valores de la simulación (‘b0’ y ‘b1’). Explore la distribución de los residuos (plot(nombre_del_modelo)) para los modelos de error ‘normal’ y ‘exponencial’.

Se supone que la colinealidad (la presencia de predictores correlacionados) afecta a la estabilidad de la regresión múltiple. El siguiente código crea dos predictores altamente colineales (‘x1’ y ‘x2’). La última línea de código muestra la correlación entre ellos.

Código

# definir semilla
set.seed(123)

# numero de observaciones
n <- 50
b0 <- -4
b1 <- 3
b2 <- -2
error <- rnorm(n = n, mean = 0, sd = 1)

#  variables aleatorias
x1 <-  rnorm(n = n, mean = 0, sd = 1)

# hacer x2 muy parecido a x1
x2 <- x1 + rnorm(n = n, mean = 0, sd = 0.3)

cor(x1, x2)

[1] 0.94642

Construya un modelo de regresión múltiple para estos datos (y ~ x1 + x2). Puede utilizar el mismo código que en la sección Añadir más de un predictor: regresión múltiple.
¿Cómo afecta a la inferencia la presencia de predictores colineales? Realice los gráficos de diagnóstico para este modelo (plot(nombre_del_modelo)).
Simule un juego de datos con tres predictores en los que sólo dos de ellos son altamente colineales. Ajuste un modelo de regresión múltiple (y ~ x1 + x2 + x3) para esos datos y observe cómo la colinealidad afecta a la estimación del predictor no colineal.

7 Extender los modelos lineales a estructuras de datos más complejas

7.1 Modelos lineales generalizados (GLM)

Los GLM nos permiten modelar la asociación a variables respuesta que no siguen una distribución normal. Además, permiten modelar distribuciones que se asemejan más al proceso que generó los datos. El siguiente código crea un juego de datos con una respuesta que representa cuentas (por lo tanto, no normal):

Código

# definir semilla
set.seed(1234)

# numero de muestas
n <- 50

# coeficientes
b0 <- 1.2
b1 <- 1.3
b2 <- 0


#generar variables
y <- rpois(n = n, lambda = 6.5) # lambda = tasa promedio de exitos 
x2 <- seq(-0.5, 0.5, , length(y))
x1 <- (log(y) - b0 - b2 * x2) / b1

# crear data frame
xy_datos_pois <- data.frame(x1, x2, y)

head(xy_datos_pois)

x1	x2	y
0.14330	-0.50000	4
0.57378	-0.47959	7
0.57378	-0.45918	7
0.57378	-0.43878	7
0.76710	-0.41837	9
0.57378	-0.39796	7

También grafiquemos ‘x1’ vs ‘y’:

Código

# graficar
ggplot(xy_datos_pois, aes(x = x1, y = y)) + 
  geom_point()

La relación no parece muy lineal ni la varianza parece ser constante a través de ‘x1’.

Podemos relajar el requisito de la distribución normal con GLMs. glm() es una función básica de R que nos ayuda a hacer el truco. Para este ejemplo la distribución más apropiada es Poisson. Esto se puede establecer en el argumento `family’ así:

Código

glm_pois  <- glm(formula = y ~ x1 + x2, data = xy_datos_pois, family = poisson())

Que es equivalente a esto:

$\hat{Y} \sim_{poisson} \beta_{o} + \beta_{1} * x_{1} + \beta_{2} * x_{2}$

Como puede ver el único argumento extra comparado con lm() es family. El resto es simplemente la “fórmula” y los “datos” con los que ya estamos familiarizados. Así que, de nuevo, podemos aprovechar nuestros conocimientos sobre los modelos lineales para extenderlos a estructuras de datos más complejas.

También necesitamos ejecutar summary() para obtener el resultado del modelo:

Código

summary(glm_pois)


Call:
glm(formula = y ~ x1 + x2, family = poisson(), data = xy_datos_pois)

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
-1.49e-08  -1.49e-08   0.00e+00   0.00e+00   2.98e-08  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 1.20e+00   1.29e-01    9.31  < 2e-16 ***
x1          1.30e+00   2.20e-01    5.92  3.2e-09 ***
x2          1.39e-16   1.90e-01    0.00        1    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3.9471e+01  on 49  degrees of freedom
Residual deviance: 2.2204e-16  on 47  degrees of freedom
AIC: 186.9

Number of Fisher Scoring iterations: 3

Interpretación del modelo

Cuadro de los coeficientes (estimados):

              Estimate Std. Error    z value   Pr(>|z|)
(Intercept) 1.2000e+00    0.12889 9.3106e+00 1.2709e-20
x1          1.3000e+00    0.21960 5.9200e+00 3.2198e-09
x2          1.3882e-16    0.19016 7.3003e-16 1.0000e+00

El modelo nos dice que $\beta_1$ (el tamaño de efecto de ‘x1’) es 1.3 y que es significativamente diferente de 0 (p = 3.21981^{-9}). Esto se interpreta en realidad como que un aumento de 1 unidad de ‘x1’ resulta en ‘y’ (tasa) en un factor de exp(1.3) = 3.6693.
El modelo también nos dice que $\beta_2$ (el tamaño de efecto de ‘x2’) es 1.38824^{-16} y que es significativamente diferente de 0 (p = 1). Esto significa que un aumento en 1 unidad de ‘x2’ resulta en ‘y’ (tasa) en un factor de exp(1.38824^{-16}) = 1.

Caso de estudio

Topp, E. N., Tscharntke, T., & Loos, J. (2022). Fire and landscape context shape plant and butterfly diversity in a South African shrubland. Diversity and Distributions, 28(3), 357-371.

“La riqueza de especies de mariposas es de tres a cuatro veces mayor cuando aumenta el hábitat natural en el paisaje circundante (en un radio de 2 km), mientras que la abundancia de mariposas se asociaba negativamente con el aumento del tiempo transcurrido desde el último incendio.”

Ejercicio 8

Intente correr un modelo lm() (con una distribución gaussiana), compare los resultados y compruebe los residuos (plot_model(model_name, type = "diag"))

Existen muchas otras funciones de distribución y enlace:

7.2 Modelos de efectos variables (modelos mixtos)

A veces nuestros conjuntos de datos incluyen estructuras con varios niveles de organización. Por ejemplo, cuando tomamos muestras de varios individuos de diferentes poblaciones. En esos casos, la variación en el nivel estructural superior (poblaciones) podría impedir la detección de patrones en el nivel inferior (individuos).

Vamos a simular algunos datos que se asemejan a ese escenario. Tenemos dos predictores continuos (x1) y una respuesta continua (y). Cada muestra procede de 1 de 8 poblaciones diferentes (poblacion):

Código

# x<- 1
# definir semilla
set.seed(28)

# numero de observaciones
n <- 300
b0 <- 1
b1 <- 1.3
poblacion <- sample(0:8, size = n, replace = TRUE)
error <- rnorm(n = n, mean = 0, sd = 2)

#  variables aleatorias
x1 <-  rnorm(n = n, mean = 0, sd = 1)
y <- b0 + poblacion * 2 + b1 * x1 + error

# add letters
poblacion <- letters[poblacion + 1]

# create data set
xy_datos_mixto <- data.frame(x1, y, poblacion)

head(xy_datos_mixto, 10)

x1	y	poblacion
1.54607	3.3378	a
0.56834	3.1399	a
-0.63855	15.2824	i
0.98424	2.6610	a
1.05355	3.7312	b
-0.82294	3.0189	a
-1.32373	3.2186	c
2.17854	19.0542	i
2.26484	15.4549	i
-0.77843	11.6670	h

Podemos explorar la relación entre ‘y’ y ‘x1’ con un gráfico:

Código

ggplot(data = xy_datos_mixto, aes(x = x1, y = y)) + 
  geom_point()

¿Puede ver claramente el patrón de asociación entre las dos variables que hemos utilizado para simular los datos? Podemos seguir explorando los datos con un modelo de regresión lineal simple:

Código

summary(lm(y ~ x1, data = xy_datos_mixto))


Call:
lm(formula = y ~ x1, data = xy_datos_mixto)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.755  -5.119   0.186   4.602  12.409 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    8.852      0.325   27.21   <2e-16 ***
x1             0.633      0.329    1.92    0.055 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.63 on 298 degrees of freedom
Multiple R-squared:  0.0123,    Adjusted R-squared:  0.00895 
F-statistic:  3.7 on 1 and 298 DF,  p-value: 0.0554

A pesar de haber simulado un $\beta_1$ distinto de cero, no tenemos ninguna asociación significativa según este modelo y la estimación de $\beta_1$ está muy lejos de la simulada. Esta pobre inferencia se debe a que estamos ignorando una característica importante de nuestros datos, la agrupación de las muestras en “poblaciones”.

Los modelos de efectos mixtos (también conocidos como modelos multinivel o modelos de efectos variables) pueden ayudarnos a tener en cuenta estas características adicionales, mejorando significativamente nuestro poder de inferencia. Coloreemos cada una de las poblaciones para ver cómo covarían las variables en cada subgrupo de datos:

Código

ggplot(data = xy_datos_mixto, aes(x = x1, y = y, color = poblacion)) + 
  geom_point()

Parece haber un patrón claro de asociación positiva entre x1 e y. El patrón se hace un poco más evidente si mostramos cada población en su propio panel:

Código

ggplot(data = xy_datos_mixto, aes(x = x1, y = y, color = poblacion)) +
  geom_point() +
  facet_wrap( ~ poblacion) +
  geom_smooth(method = "lm", se = FALSE)

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

Construyamos un modelo de efectos mixtos utilizando la población como intercepto variable:

Código

mix_eff_mod <- lmer(formula = y ~ x1 + (1 | poblacion), data = xy_datos_mixto)

Que es equivalente a esto:

$\hat{Y} \sim_{gausiana} (\beta_{o} + \beta_{grupo}) + \beta_{1} * x_{1}$

Podemos ver el resultado del modelo igual que con un modelo lm():

Código

summary(mix_eff_mod)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: y ~ x1 + (1 | poblacion)
   Data: xy_datos_mixto

REML criterion at convergence: 1296.4

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.524 -0.657 -0.032  0.612  3.245 

Random effects:
 Groups    Name        Variance Std.Dev.
 poblacion (Intercept) 30.07    5.48    
 Residual               3.76    1.94    
Number of obs: 300, groups:  poblacion, 9

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)    8.788      1.831   8.003     4.8   0.0014 ** 
x1             1.317      0.115 290.068    11.4   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
   (Intr)
x1 -0.003

El modelo detectó correctamente el patrón simulado y la estimación de $\beta_1$ (1.317) es muy cercana al valor simulado.

Caso de estudio

Araya-Salas M, P González-Gómez, K Wojczulanis-Jakubas, V López III & T Wright. 2018. Spatial memory is as important as weapon and body size for territorial ownership in a lekking hummingbird. Scientific Reports. 13, e0189969

“… el tamaño corporal mostró una correlación negativa con la frecuencia baja del canto.”

Interpretación del modelo

Cuadro con los coeficientes (estimados):

            Estimate Std. Error       df t value   Pr(>|t|)
(Intercept)   8.7881     1.8315   8.0031  4.7983 1.3569e-03
x1            1.3166     0.1151 290.0679 11.4395 2.9326e-25

El modelo encontró que $\beta_1$ (el tamaño de efecto de ‘x1’) es 1.31665 y que es significativamente diferente de 0 (p = 2.93263^{-25})
Los valores simulados de los parámetros de regresión pueden compararse con el resumen del modelo lmer() para tener una idea de la precisión del modelo:
- $\beta_1$ se fijó en 1.3 y se estimó como 1.317.
La varianza del intercepto para cada población fue de 29.95624 y la desviación estándar de 5.47323
Los interceptos estimados para cada población son:

Código

ranef(mix_eff_mod)$poblacion

	(Intercept)
a	-8.22835
b	-6.17134
c	-4.23048
d	-1.30838
e	0.24573
f	2.22552
g	3.91556
h	5.75609
i	7.79566

Podemos compararlos con los promedios para cada población en los datos:

Código

# sumar interceptos por poblacion a el intercepto total
interceptos <- ranef(mix_eff_mod)$poblacion + fixef(mix_eff_mod)[1]

# calcular los promedios por poblacion en los datos 
prom_pobs <- aggregate(y ~ poblacion, xy_datos_mixto, mean)

# unir
cbind(interceptos, prom_pobs= prom_pobs[,2])

	(Intercept)	prom_pobs
a	0.55975	0.86151
b	2.61677	2.76860
c	4.55762	4.90112
d	7.47972	7.73174
e	9.03384	8.87649
f	11.01363	10.77650
g	12.70366	12.62651
h	14.54420	14.62110
i	16.58376	16.42520

La varianza entre los niveles del factor aleatorio (poblacion) es de 29.95624 y la desviación estándar de 5.47323

7.2.1 La paradoja de Simpson

La paradoja de Simpson es un fenómeno en estadística en el que una tendencia aparece en varios grupos de datos pero desaparece o se invierte cuando se combinan los grupos. Se ha utilizado para ilustrar el tipo de resultados engañosos que puede generar el ignorar la estructura en los datos y las relaciones causales entre variables.

Tomado de wikipedia

La estructura de los datos debido a observaciones pertenecientes al mismo grupo puede ser tomada en cuenta usando el grupo como un factor aleatorio, tal y como se hizo en el ejemplo anterior.

Modelos lineales como abordaje estadístico estándar

$\hat{Y} \sim_{función\ de\ enlace} (\beta_{o} + \beta_{grupo}) + \beta_{1} * x_{1} + \cdots + \beta_{n} * x_{n}$

Este modelo básico incluye una función de enlace que cuando es gausiana (i.e. normal) equivale a un modelo lineal y cuando no, a un modelo generalizado
Cuando el $\beta_{grupo}$ es diferente entre grupos (i.e. cuando el intercepto difiere entre grupos) es un modelo mixto con intercepto variable (i.e. aleatorio)

Este manual solo pretende sugerir un abordaje estadístico centrado en los modelos de regresión desde el cual se pueden llevar a cabo adecuadamente la mayoría de los análisis estadísticos que usamos en biología. El manual no tiene como objetivo explicar en detalle las particularidades de los modelos generalizados y mixtos. Note que estos modelos, principalmente los mixtos, pueden adaptarse a otras estructuras de datos mas complejas no cubiertas aquí, como lo son las pendientes aleatorias, pendientes e interceptos aleatorios, estructuras de correlación (i.e. autocorrelación espacial o temporal, pedigrís o árboles filogenéticos), medidas repetidas y el uso de pseudoreplicas, entre otras.

7.3 Referencias

Información de la sesión

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=es_CR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=es_CR.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=es_CR.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] car_3.1-2         carData_3.0-5     sjPlot_2.8.15     lmerTest_3.1-3   
 [5] lme4_1.1-35.1     Matrix_1.4-0      viridis_0.6.2     viridisLite_0.4.2
 [9] ggplot2_3.4.4     knitr_1.45       

loaded via a namespace (and not attached):
 [1] tidyr_1.2.0         jsonlite_1.8.7      splines_4.1.2      
 [4] modelr_0.1.11       assertthat_0.2.1    yaml_2.3.7         
 [7] remotes_2.4.2.1     bayestestR_0.13.1   numDeriv_2016.8-1.1
[10] pillar_1.9.0        backports_1.4.1     lattice_0.20-45    
[13] glue_1.6.2          digest_0.6.33       minqa_1.2.4        
[16] colorspace_2.1-0    sandwich_3.0-1      htmltools_0.5.7    
[19] pkgconfig_2.0.3     broom_0.7.12        purrr_0.3.4        
[22] xtable_1.8-4        mvtnorm_1.1-3       scales_1.2.1       
[25] emmeans_1.8.9       tibble_3.2.1        mgcv_1.8-39        
[28] farver_2.1.1        generics_0.1.2      sjlabelled_1.2.0   
[31] sketchy_1.0.2       TH.data_1.1-0       withr_2.5.2        
[34] cli_3.6.1           survival_3.2-13     magrittr_2.0.3     
[37] crayon_1.5.2        estimability_1.4.1  evaluate_0.23      
[40] fansi_1.0.5         nlme_3.1-155        MASS_7.3-55        
[43] tools_4.1.2         lifecycle_1.0.4     multcomp_1.4-18    
[46] munsell_0.5.0       ggeffects_1.3.2     packrat_0.9.2      
[49] compiler_4.1.2      rlang_1.1.2         grid_4.1.2         
[52] nloptr_2.0.0        rstudioapi_0.13     htmlwidgets_1.5.4  
[55] labeling_0.4.3      rmarkdown_2.25      xaringanExtra_0.7.0
[58] boot_1.3-28         gtable_0.3.4        codetools_0.2-18   
[61] abind_1.4-5         sjstats_0.18.2      DBI_1.1.3          
[64] sjmisc_2.8.9        R6_2.5.1            gridExtra_2.3      
[67] zoo_1.8-9           dplyr_1.0.8         performance_0.8.0  
[70] fastmap_1.1.1       utf8_1.2.4          insight_0.19.6     
[73] stringi_1.7.12      Rcpp_1.0.11         vctrs_0.6.4        
[76] tidyselect_1.1.2    xfun_0.41           coda_0.19-4

--- title: Modelos de regresión --- ```{=html} <style> body { counter-reset: source-line 0; } pre.numberSource code { counter-reset: none; } </style> ``` ```{r,echo=FALSE,message=FALSE} options("digits"=5) options("digits.secs"=3) ``` ::: {.alert .alert-info} # Objetivos del manual {.unnumbered .unlisted} - Familiarizarse con el formato para reportes dinámicos Rmarkdown/quarto - Ser capaz de documentar el manejo y análisis de datos en R usando Rmarkdown/quarto ::: ```{r, echo = FALSE, message=FALSE} library(knitr) library(ggplot2) library(viridis) library(lmerTest) library(sjPlot) library(car) # ggplot settings geom_histogram <- function(...) ggplot2::geom_histogram(..., fill = viridis(10, alpha = 0.5)[8], show.legend = FALSE, bins = 20, color = "black") geom_smooth <- function(...) ggplot2::geom_smooth(..., color = viridis(10, alpha = 0.5)[8]) geom_boxplot <- function(...) ggplot2::geom_boxplot(..., fill = viridis(10, alpha = 0.5)[7]) geom_pointrange <- function(...) ggplot2::geom_pointrange(..., show.legend = FALSE, color = viridis(10, alpha = 0.5)[7], size = 2) plot_model <- function(...) sjPlot::plot_model(xy_mod, type = "diag", colors = viridis(10, alpha = 0.5)[7]) theme_set(theme_classic(base_size = 20)) ``` ```{=html} <style> body { counter-reset: source-line 0; } pre.numberSource code { counter-reset: none; } </style> ``` # Modelos lineales como marco unificador Tradicionalmente, los modelos estadísticos se han enseñado como herramientas desconectadas sin una relación clara entre ellas. Sólo hace falta darle un vistazo a un libro clásico de estadística para biólogos (Sokal & Rolhf): ```{r, out.width = "100%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/sokal_rohlf.png") ``` Sin embargo, la mayoría de esos modelos estadísticos comunes no son más que casos especiales de modelos lineales. Por lo tanto, aprenderlos como tales puede simplificar mucho las cosas. Este enfoque tiene varias ventajas: - En primer lugar, **todo se reduce a *y=a⋅x+b***, lo que simplifica en gran medida el aprendizaje - Esto también significa que no hay **necesidad de aprender sobre los parámetros, las hipótesis y la interpretación de los resultados para cada caso especial**. - Los modelos lineales se han ampliado para **tener en cuenta distribuciones y estructuras de datos complejas** (por ejemplo, modelos mixtos, modelos lineales generalizados, modelos cero-inflados, etc.) proporcionando una **plataforma más flexible**. - Los modelos lineales se **aplican en todos los paradigmas estadísticos** (por ejemplo, frecuentista, bayesiano) ::: {.alert .alert-info} # Objetivo del manual - Comprender la inferencia estadística a través de una única herramienta de modelaje (modelos lineales en sentido amplio) - Familiarizarse con la construcción de modelos lineales - Extender los modelos lineales a diferentes estructuras de datos ::: ------------------------------------------------------------------------ Paquetes a utilizar en este manual: ```{r} # instalar/cargar paquetes sketchy::load_packages( c("ggplot2", "viridis", "lmerTest", "sjPlot") ) ``` ------------------------------------------------------------------------ # Cómo simular datos ## Generación de números aleatorios en R La estadística nos permite inferir patrones en los datos. Solemos utilizar conjuntos de datos reales para enseñar estadística. Sin embargo, puede ser circular entender el funcionamiento interno de una herramienta estadística probando su capacidad para inferir un patrón que no estamos seguros de encontrar en los datos (y no tenemos idea del mecanismo que produjo ese patrón). **Las simulaciones nos permiten crear escenarios controlados en los que conocemos con seguridad los patrones** presentes en los datos y los procesos subyacentes que los han generado. R ofrece algunas funciones básicas para la simulación de datos. Las más utilizadas son las funciones generadoras de números aleatorios. Los nombres de estas funciones comienzan con *r* (`r____()`). Por ejemplo, `runif()`: ```{r, echo=FALSE} set.seed(7) ``` ```{r} # simular variable uniforme unif_var <- runif(n = 100, min = 0, max = 10) ``` El resultado es un vector numérico de longitud 100 (`n = 100`): ```{r} # imprimir variable unif_var ``` Podemos explorar el resultado graficando un histograma: ```{r} # crear histograma ggplot(data = data.frame(unif_var), mapping = aes(x = unif_var)) + geom_histogram() ``` Muestra una distribución uniforme que va de 0 a 10. También podemos simular números aleatorios procedentes de una distribución normal utilizando `rnorm()`: ```{r} # crear una variable normal norm_var <- rnorm(n = 1000, mean = 2, sd = 1) # graficar histograma ggplot(data = data.frame(norm_var), mapping = aes(x = norm_var)) + geom_histogram() ``` Tenga en cuenta que todas las funciones generadoras de números aleatorios tienen el argumento 'n', que determina la longitud del vector generado (es decir, el número de números aleatorios), además de algunos argumentos adicionales relacionados con parámetros específicos de la distribución. Las variables continuas (es decir, los vectores numéricos) pueden convertirse en variables discretas (es decir, números enteros) simplemente redondeándolas: ```{r} v1 <- rnorm(n = 5, mean = 10, sd = 3) v1 round(x = v1, digits = 0) ``` ::: {.alert .alert-info} Ejercicio 1 - ¿Qué hacen las funciones `rbinom()` y `rexp()`? - Ejecutela y haga histogramas de sus resultados - ¿Qué hacen los argumentos 'mean' y 'sd' en `rnorm()`? Juege con diferentes valores y comprueba el histograma para hacerse una idea de su efecto en la simulación ::: ## Generación de variables categóricas La forma más sencilla de generar variables categóricas es utilizar el vector de ejemplo `letters' (o`LETTERS') para asignar niveles de categoría. Podemos hacerlo utilizando la función `rep()`. Por ejemplo, el siguiente código crea un vector categórico (caracteres) con dos niveles, cada uno con 4 observaciones: ```{r} rep(letters[1:2], each = 4) ``` También podemos replicar este patrón utilizando el argumento 'times'. Este código replica el vector anterior 2 veces: ```{r} rep(letters[1:2], each = 4, times = 2) ``` Otra opción es simular una variable a partir de una distribución binomial y luego convertirla en un factor: ```{r} # correr rbinom binom_var <- rbinom(n = 50, size = 1, prob = 0.5) binom_var ``` ```{r} # convertir a factor categ_var <- factor(binom_var, labels = c("a", "b")) categ_var ``` ## Muestreo aleatorio La otra herramienta importante de R para jugar con datos simulados es `sample()`. Esta función permite tomar muestras de tamaños específicos de vectores. Por ejemplo, tomemos el ejemplo del vector 'letters': ```{r} letters ``` Podemos tomar una muestra de este vector como es: ```{r} # tomar muestra sample(x = letters, size = 10) ``` El argumento 'size' nos permite determinar el tamaño de la muestra. Tenga en cuenta que obtendremos un error si el tamaño es mayor que el propio vector: ```{r, error=TRUE} sample(x = letters, size = 30) ``` Esto sólo puede hacerse cuando el muestreo es con reemplazo (replacement). El muestreo con reemplazo puede aplicarse estableciendo el argumento `replace = TRUE`: ```{r} sample(x = letters, size = 30, replace = TRUE) ``` ## Iterar un proceso A menudo, las simulaciones deben repetirse varias veces para descartar resultados espurios debidos al azar o simplemente para probar diferentes parámetros. Las funciones de simulación de datos mencionadas anteriormente pueden ejecutarse varias veces (por ejemplo, iteradas) utilizando la función `replicate()`: ```{r} # replicar repl_rnorm <- replicate(n = 3, expr = rnorm(2), simplify = FALSE) # ver clase class(repl_rnorm) # imprimir repl_rnorm ``` ## Hacer que las simulaciones sean reproducibles El último truco que necesitamos para ejecutar simulaciones en R es la capacidad de reproducir una simulación (es decir, obtener exactamente los mismos datos y resultados simulados). Esto puede ser útil para que otros investigadores puedan ejecutar nuestros análisis exactamente de la misma manera. Esto puede hacerse fácilmente con la función `set.seed()`. Pruebe a ejecutar el siguiente código. Debería obtener la misma salida: ```{r} # definir semilla set.seed(10) # crear variable uniforme runif(n = 2) ``` ------------------------------------------------------------------------ # Crear juegos de datos ## Juegos de datos con variables numéricas y categóricas Ahora que sabemos cómo simular variables continuas y categóricas. Podemos juntarlas para crear conjuntos de datos simulados. Esto se puede hacer utilizando la función `data.frame()`: ```{r} # crear variable categorica grupo <- rep(letters[1:2], each = 3) # crear variable continuaa size <- rnorm(n = 6, mean = 5, sd = 1) # poner juntas en un data frame df <- data.frame(grupo, size) # imprimir df ``` Por supuesto, podríamos añadir más variables a este cuadro de datos: ```{r} # crear variable categorica grupo <- rep(letters[1:2], each = 3) individuo <- rep(LETTERS[1:6]) # crear variables continuas size <- rnorm(n = 6, mean = 5, sd = 1) weight <- rnorm(n = 6, mean = 100, sd = 10) # poner todo en un data frame df <- data.frame(grupo, individuo, size, weight) # imprimir df ``` Y eso es un juego de datos simulados en su forma más básica. Se parece mucho al tipo de datos con los que trabajamos en biología. ------------------------------------------------------------------------ # Cómo utilizar datos simulados para entender el comportamiento de las herramientas estadísticas ## Prueba de concepto: *el Teorema Central del Límite* El [Teorema del Límite Central](https://en.wikipedia.org/wiki/Central_limit_theorem) afirma que, si tomamos muestras aleatorias de una población, los promedios de esas muestras seguirán una distribución normal, aunque la población no esté distribuida normalmente. Además, la distribución normal resultante debe tener un promedio cercano al promedio de la población. El teorema es un concepto clave para la estadística inferencial, ya que implica que los métodos estadísticos que funcionan para las distribuciones normales pueden ser aplicables a muchos problemas que implican otros tipos de distribuciones. No obstante, el objetivo aquí es sólo mostrar cómo se pueden utilizar las simulaciones para entender el comportamiento de los métodos estadísticos. Para comprobar si esas afirmaciones básicas sobre el Teorema del Límite Central son ciertas, podemos utilizar datos simulados en R. Vamos a simular una población de 1000 observaciones con una distribución uniforme: ```{r, eval=FALSE} # simular popublacion uniforme unif_pop <- runif(1000, min = 0, max = 10) # ver histograma ggplot(data = data.frame(unif_pop), mapping = aes(x = unif_pop)) + geom_histogram() ``` ```{r, echo=FALSE} set.seed(10) # simulate uniform population unif_pop <- runif(1000, min = 0, max = 10) # ver distribucion ggplot(data = data.frame(unif_pop), mapping = aes(x = unif_pop)) + geom_histogram() ``` Podemos tomar muestras aleatorias usando `sample()` así: ```{r} sample(x = unif_pop, size = 30) ``` Este proceso puede ser replicado varias veces con `replicate()`: ```{r} # replicar samples <- replicate(n = 100, expr = mean(sample(x = unif_pop, size = 30))) ``` El código anterior toma 100 muestras con 30 valores cada una. Ahora podemos comprobar la distribución de las muestras: ```{r, eval=FALSE} # ver distribucion/ histograma ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram() ``` ```{r, echo=FALSE} # ver distribucion ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram() ``` ... asi como el promedio: ```{r} mean(samples) ``` Como era de esperar, las muestras siguen una distribución normal con una media cercana a la media de la población, que es: ```{r} mean(unif_pop) ``` Probemos con una distribución más compleja. Por ejemplo, una distribución bimodal: ```{r, eval=FALSE} # usar semilla set.seed(123) # simular variables norm1 <- rnorm(n = 1000, mean = 10, sd = 3) norm2 <- rnorm(n = 1000, mean = 20, sd = 3) # juntar en una sola variable bimod_pop <- c(norm1, norm2) # ver histograma ggplot(data = data.frame(bimod_pop), mapping = aes(x = bimod_pop)) + geom_histogram() ``` ```{r, echo=FALSE} # definir semilla set.seed(123) norm1 <- rnorm(n = 1000, mean = 10, sd = 3) norm2 <- rnorm(n = 1000, mean = 20, sd = 3) bimod_pop <- c(norm1, norm2) # ver distribucion ggplot(data = data.frame(bimod_pop), mapping = aes(x = bimod_pop)) + geom_histogram() ``` ```{r, eval=FALSE} # replicar muestreo samples <- replicate(200, mean(sample(bimod_pop, 10))) # ver histograma ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram() ``` ```{r, echo=FALSE} samples <- replicate(200, mean(sample(bimod_pop, 10))) # ver distribucion ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram() ``` ```{r} # ver promedios mean(samples) mean(bimod_pop) ``` ::: {.alert .alert-info} Ejercicio 2 - Intenta explorar el Teorema del Límite Central como en el caso anterior, pero esta vez utilizando: 1. Una distribución exponencial (`rexp()`) 2. Una distribución binomial (`rbinom()`) - Para cada distribución: grafique un histograma y compare los promedios de la población y de las muestras ::: ------------------------------------------------------------------------ ::: {.alert .alert-warning} <center>Regresión lineal simple</center> Las regresiones lineales se basan en la ecuación lineal *a = mx + b* que aprendimos en el colegio. La representación formal tiene este aspecto: <center>$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1}$</center> - $\hat{Y}$: variable respuesta - $\beta_{o}$: intercepto - $\beta_{1}$: estimado de la magnitud del efecto de $x_{1}$ en $\hat{Y}$ (también conocido como tamaño de efecto, coeficiente o simplemente "estimado") - $x_{1}$: variable predictora El objetivo más común de una regresión lineal es la estimación de los valores de $\beta_{*}$. Esto se consigue encontrando la línea recta que mejor se ajusta y que representa la asociación entre un predictor y la respuesta: ```{r, out.width = "100%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/regression.jpg") ``` Estos $\beta_{*}$ son el tamaño de efecto estimado del predictor correspondiente (por ejemplo, $\beta_{1}$ es el tamaño de efecto $x_{1}$). Su valor representa el cambio medio en $\hat{Y}$ (en unidades $\hat{Y}$) para una unidad de cambio en $\beta_{*}$. Por lo tanto, la hipótesis nula es que esos $\beta_{*}$ no son diferentes de 0: <center>$\qquad \mathcal{H}_0: \hat{Y} = \beta_0 + 0 * x_{1}$</center> lo que equivale a esto: <center>$\qquad \mathcal{H}_0: \hat{Y} = \beta_0$</center> ```{r, echo = FALSE, eval = FALSE} For instance, a regression model with this output: set.seed(123) # numero de observaciones n <- 50 b0 <- 1 b1 <- 3 error <- rnorm(n = n, sd = 2) # variables aleatorias B1 <- rnorm(n = n, mean = 0, sd = 1) y <- b0 + b1 * B1 + error # crear data frame xy_datos <- data.frame(B1, y) # construir modelos xy_mod <- lm(formula = y ~ B1, data = xy_datos) summary(xy_mod)$coefficients Results in the following predictive model: <center>$\hat{Y} \sim 1.079 + 2.926 * x_{1}$</center> ``` Por ejemplo, un modelo de regresión con este resultado: ```{r, out.width = "90%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/model_betas_esp.png") ``` ::: ------------------------------------------------------------------------ # Regresión lineal en R Para sacar el máximo provecho de los modelos lineales necesitamos sentirnos cómodos con ellos. Empezaremos explorando la función de regresión lineal de R `lm()`. En R la mayoría de los modelos lineales y sus extensiones comparten formatos comunes de entrada de datos y resultados, lo que facilita su aplicación una vez que entendemos sus fundamentos. Utilizaremos el juego de datos 'trees' que viene por defecto con R. 'trees' proporciona medidas del **diámetro** (etiquetado como 'Girth'), **altura** y **volumen** de 31 cerezos talados: ```{r} head(trees) ``` La función básica de R para construir un modelo lineal es `lm()`. Veamos los componentes básicos de un modelo de regresión utilizando `lm()`: ```{r, out.width = "100%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/lm.png") ``` Podemos correr este modelo para ver el resultado: ```{r} reg_mod <- lm(formula = Height ~ Girth, data = trees) summary(reg_mod) ``` Esto es lo que significan los elementos de la salida: - **Llamada** (call): la función y los parámetros que se utilizaron para crear el modelo - **Residuales** (residuals): distribución de los residuales. Los residuales son la diferencia entre lo que predijo el modelo y el valor real de *y*. Esta es una representación gráfica de los residuales: ```{r, echo = FALSE} # añadir valores predichos trees$predicted <- reg_mod$fitted.values # graficar ggplot(trees, aes(x = Girth, y = Height)) + geom_smooth(method = "lm", se = FALSE) + labs(x = "Diámetro (Girth)", y = "Alto (height)") + geom_segment(aes(xend = Girth, yend = predicted), alpha = .2) + geom_point() ``` - **Coeficientes** (coefficients): contiene los tamaños de efecto ('Estimates'), una medida de su incertidumbre ('Standar error'), la estadística asociada ('t value') y el valor p ('Pr(\>\|t\|)'). Los tamaño de efecto o estimados se dan como el cambio promedio en *y* por cada aumento de 1 unidad en *x*. Para este ejemplo el estimado es de 1,0544 $in / in$. A pesar de que las unidades se anulan (\`1,0544 $in / in$ = 1,0544) tenerlas en cuenta sigue siendo biológicamente significativo. Significan que, en promedio, y un aumento de 1 pulgada en la circunferencia que se espera un aumento de 1,0544 pulgadas en la altura. - **Error estándar de los residuales** (residual standard error): se explica por sí mismo. El error estándar de los residuales - **R-cuadrado múltiple** (multiple R-squared): el coeficiente de determinación, que pretende ser una medida de lo bien que su modelo se ajusta a los datos - **R-cuadrado ajustado** (adjusted R-squared): similar al 'R-cuadrado múltiple' pero penalizado por el número de parámetros - **Estadístico F** (F-statistic): estadístico para una prueba global que comprueba si al menos uno de sus coeficientes es distinto de cero - **Valor de p**: probabilidad de una prueba global que comprueba si al menos uno de sus coeficientes es distinto de cero Utilizaremos `lm()` para mostrar la flexibilidad de los modelos de regresión. Los componentes de regresión se añadirán gradualmente para que podamos tomarnos el tiempo de entender cada uno de ellos así como los correspondientes cambios en la salida de la regresión. ------------------------------------------------------------------------ ## Modelo solo con intercepto Primero vamos a crear una variable numérica de respuesta: ```{r} # definir semilla set.seed(123) # numero de observaciones n <- 50 # variables aleatorias y <- rnorm(n = n, mean = 0, sd = 1) # put it in a data frame y_data <- data.frame(y) ``` Esta única variable puede introducirse en un **modelo de regresión sólo con intercepto**. Para ello, debemos suministrar la fórmula del modelo y los datos a `lm()`: ```{r} # run model y_mod <- lm(formula = y ~ 1, data = y_data) ``` Lo que equivale a: <center>$\hat{Y} \sim \beta_{o}$</center> Podemos obtener el resumen por defecto de los resultados del modelo ejecutando `summary()` en el objeto de salida 'y_mod': ```{r} summary(y_mod) ``` Puede ser bastante informativo graficar los tamaños de efecto (aunque en este caso sólo tenemos uno): ```{r} ci_df <- data.frame(param = names(y_mod$coefficients), est = y_mod$coefficients, confint(y_mod)) ggplot(ci_df, aes(x=param, y=est)) + geom_hline(yintercept = 0, color="red", lty = 2) + geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + labs(x = "Parámetro", y = "Tamaño de efecto") + coord_flip() ``` ::: {.alert .alert-success} <center>Interpretación del modelo</center> Para evaluar la importancia de la asociación nos centramos en la tabla de coeficientes: ```{r, echo= FALSE} summ_y_mod <- summary(y_mod) summ_y_mod$coefficients ``` - En este ejemplo no hay predictores en el modelo, por lo que sólo tenemos una estimación para el intercepto ($\beta_0$) - El modelo nos dice que el intercepto se estima en `r summ_y_mod$coefficients[1, 1]` y que este valor no es significativamente diferente de 0 (valor p = `r summ_y_mod$coefficients[1, 4]`) - En este caso el intercepto es simplemente la media de la variable de respuesta ```{r} mean(y_data$y) ``` - Rara vez tenemos predicciones sobre el intercepto, por lo que tendemos a ignorar esta estimación. ::: ::: {.alert .alert-warning} <center>Caso de estudio</center> - Ruhs, E. C., Martin, L. B., & Downs, C. J. (2020). [**The impacts of body mass on immune cell concentrations in birds**](https://royalsocietypublishing.org/doi/pdf/10.1098/rspb.2020.0655?download=true). *Proceedings of the Royal Society B*, 287(1934), 20200655. "*Encontramos que un modelo sólo con intercepto (intercept-only model) explicaba mejor las concentraciones de linfocitos y eosinófilos en las aves, indicando que las concentraciones de estos tipos de células eran independientes de la masa corporal.*" ```{r, out.width = "85%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/leukocytes.png") ``` ::: ::: {.alert .alert-info} Ejercicio 3 - Cambie el argumento 'mean' en la llamada de la función "rnorm()\` ([línea 130](#cb53-8)) por un valor distinto de 0 y observe cómo cambian los valores en la tabla de coeficientes - Cambia el argumento `sd` en la llamada a la función `rnorm()` ([línea 130](#cb53-8)) por un valor más alto y observe cómo cambian los valores en la tabla de coeficientes ::: ------------------------------------------------------------------------ ## Añadir un predictor no asociado Podemos crear 2 variables numéricas no relacionadas así: ```{r} # definir semilla set.seed(123) # numero de observaciones n <- 50 # variables aleatorias y <- rnorm(n = n, mean = 0, sd = 1) x1 <- rnorm(n = n, mean = 0, sd = 1) # crear data frame xy_datos <- data.frame(x1, y) ``` Estas dos variables pueden introducirse en un modelo de regresión para evaluar la asociación entre ellas: ```{r} # construir model xy_mod <- lm(formula = y ~ x1, data = xy_datos) # graficar ggplot(xy_datos, aes(x = x1, y = y)) + geom_smooth(method = "lm", se = FALSE) + geom_point() # graficar points ``` Que es equivalente a esto: <center>$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1}$</center> Imprimamos el resumen de este modelo: ```{r} summary(xy_mod) ``` ... y graficar los tamaños de efecto: ```{r} ci_df <- data.frame(param = names(xy_mod$coefficients), est = xy_mod$coefficients, confint(xy_mod)) ggplot(ci_df, aes(x=param, y=est)) + geom_hline(yintercept = 0, color="red", lty = 2) + geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + labs(x = "Parámetro", y = "Tamaño de efecto") + coord_flip() ``` Deberíamos "diagnosticar" la idoneidad del modelo inspeccionando más de cerca la distribución de los residuales. La función `plot_model()` del paquete `sjPlot` hace un buen trabajo para crear gráficos de diagnóstico para modelos lineales: ```{r, out.width= "80%", fig.height=4} plot_model(xy_mod, type = "diag") ``` ::: {.alert .alert-success} <center>Interpretación del modelo</center> Cuadro con coeficientes: ```{r, echo= FALSE} summ_xy_mod <- summary(xy_mod) summ_xy_mod$coefficients ``` - En este ejemplo hemos añadido un predictor al modelo, por lo que hemos obtenido una estimado adicional (y una fila extra, 'x1') - El modelo nos dice que la estimación de 'x1' es `r summ_xy_mod$coefficients[2, 1]` y que no es significativamente diferente de 0 (p = `r summ_xy_mod$coefficients[2, 4]`) ::: ------------------------------------------------------------------------ ## Simular un predictor asociado {#simular-un-predictor-asociado} Podemos utilizar la fórmula del modelo lineal anterior para simular dos variables continuas asociadas así: ```{r} # definir semilla set.seed(123) # numero de observaciones n <- 50 b0 <- -4 b1 <- 3 error <- rnorm(n = n, sd = 3) # variables aleatorias x1 <- rnorm(n = n, mean = 0, sd = 1) y <- b0 + b1 * x1 + error # crear data frame xy_datos2 <- data.frame(x1, y) ``` Note que **también hemos añadido un término de error**, por lo que la asociación no es perfecta. Vamos a correr el modelo y graficar la asociación entre las dos variables: ```{r} # construir model xy_mod2 <- lm(formula = y ~ x1, data = xy_datos2) # graficar ggplot(xy_datos2, aes(x = x1, y = y)) + geom_smooth(method = "lm", se = FALSE) + geom_point() # graficar points ``` La fórmula es la misma que la del modelo anterior: <center>$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1}$</center> Este es el resumen del modelo: ```{r} summary(xy_mod2) ``` .. el gráfico con los tamaños de efecto: ```{r} ci_df <- data.frame(param = names(xy_mod2$coefficients), est = xy_mod2$coefficients, confint(xy_mod2)) ggplot(ci_df, aes(x=param, y=est)) + geom_hline(yintercept = 0, color="red", lty = 2) + geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + labs(x = "Parámetro", y = "Tamaño de efecto") + coord_flip() ``` ... y los gráficos diagnósticos del modelo: ```{r, out.width= "80%", fig.height=4} plot_model(xy_mod2, type = "diag") ``` ::: {.alert .alert-success} <center>Interpretación del modelo</center> Cuadro con los coeficientes (estimados): ```{r, echo= FALSE} # resumen summ_xy_mod2 <- summary(xy_mod2) summ_xy_mod2$coefficients ``` - El modelo nos dice que $beta_1$ (el tamaño de efecto de 'x1') es `r summ_xy_mod2$coefficients[2, 1]` y que es significativamente diferente de 0 (p = `r summ_xy_mod2$coefficients[2, 4]`) - Los valores simulados de los parámetros de regresión pueden compararse con el resumen del modelo `lm()` para tener una idea de la precisión del modelo: - $\beta_1$ (the effect size of 'x1') was set to `r b1` and was estimated as `r round(summ_xy_mod2$coefficients[2, 1], 3)` by the model ::: ::: {.alert .alert-warning} <center>Caso de estudio</center> - Keenan EL, Odom KJ, Araya-Salas M, Horton KG, Strimas-Mackey M, Meatte MA, Mann NI, Slater PJ, Price JJ, and Templeton CN. 2020. [**Breeding season length predicts duet coordination and consistency in Neotropical wrens (Troglodytidae)**](https://royalsocietypublishing.org/doi/pdf/10.1098/rspb.2020.2482). *Proceeding of the Royal Society B*. 20202482 ```{r, out.width = "100%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/duets.png") ``` "*... la coordinación y la consistencia de los duetos son mayores en las especies con temporadas de apareamiento especialmente largas.*" ```{r, out.width = "70%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/duet_reg.png") ``` ::: ::: {.alert .alert-info} Ejercicio 4 - Aumente el tamaño de la muestra (`n`) a 1000 o más - ¿Cómo cambiaron los estimados del tamaño de efecto ($\beta$)? - ¿Cómo cambió el error estándar del tamaño de efecto? - Ahora cambie `n` a 15 y compruebe de nuevo las estimaciones del modelo (esta vez compruebe también el valor p) ::: ------------------------------------------------------------------------ ## Varios predictores: regresión múltiple La regresión lineal múltiple es una extensión del modelo de regresión lineal simple que puede tomar varios predictores: <center>$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1} + \cdots + \beta_{n} * x_{n}$</center> La fórmula parece compleja, pero sólo quiere decir que cualquier parámetro adicional tendrá su propia estimación ($\beta$). La fórmula para una regresión lineal de dos predictores se ve así: <center>$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1} + \beta_{2} * x_{2}$</center> ... y se puede simular así: ```{r} # semila set.seed(123) # numero de observaciones n <- 50 b0 <- -4 b1 <- 3 b2 <- -2 error <- rnorm(n = n, mean = 0, sd = 3) # variables aleatorias x1 <- rnorm(n = n, mean = 0, sd = 1) x2 <- rnorm(n = n, mean = 0, sd = 1) y <- b0 + b1 * x1 + b2 * x2 + error # crear un data frame xy_datos_multp <- data.frame(x1, x2, y) # construir el modelo xy_mod_multp <- lm(formula = y ~ x1 + x2, data = xy_datos_multp) summary(xy_mod_multp) ``` ... graficar los tamaños de efecto: ```{r} ci_df <- data.frame(param = names(xy_mod_multp$coefficients), est = xy_mod_multp$coefficients, confint(xy_mod_multp)) ggplot(ci_df, aes(x=param, y=est)) + geom_hline(yintercept = 0, color="red", lty = 2) + geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + labs(x = "Parámetro", y = "Tamaño de efecto") + coord_flip() ``` ... y los gráficos diagnósticos: ```{r, out.width= "80%", fig.height=4} plot_model(xy_mod_multp, type = "diag") ``` ::: {.alert .alert-success} <center>Interpretación del modelo</center> Cuadro con los coeficientes (estimados): ```{r, echo= FALSE} summ_xy_mod_multp <- summary(xy_mod_multp) summ_xy_mod_multp$coefficients ``` - El modelo encontró que $\beta_1$ (el tamaño de efecto de 'x1') es `r summ_xy_mod_multp$coefficients[2, 1]` y que es significativamente diferente de 0 (p = `r summ_xy_mod_multp$coefficients[2, 4]`). Este es el efecto de 'x1' sobre 'y' una vez removida la variación explicada por 'x2'. - También se encontró que el $\beta_2$ (el tamaño de efecto de 'x2') es `r summ_xy_mod_multp$coefficients[3, 1]` y que también es significativamente diferente de 0 (p = `r summ_xy_mod_multp$coefficients[3, 4]`). Este es el efecto de 'x2' sobre 'y' una vez removida la variación explicada por 'x1'. - Los valores simulados de los parámetros de regresión pueden compararse con el resumen del modelo `lm()` para tener una idea de la precisión del modelo: - $\beta_1$ se fijó en `r b1` y se estimó como `r round(summ_xy_mod_multp$coefficients[2, 1], 3)`. - $\beta_2$ (el tamaño de efecto de 'x2') se fijó en `r b2` y se estimó como `r round(summ_xy_mod_multp$coefficients[3, 1], 3)`. ::: ::: {.alert .alert-warning} <center>Caso de estudio</center> - Araya-Salas M, P González-Gómez, K Wojczulanis-Jakubas, V López III & T Wright. 2018. [**Spatial memory is as important as weapon and body size for territorial ownership in a lekking hummingbird**](https://www.nature.com/articles/s41598-018-20441-x). *Scientific Reports*. 13, e0189969 "*La memoria espacial, el tamaño corporal y el largo de la punta del pico ... predijeron positivamente la probabilidad de adquirir y defender un territorio.*" ```{r, out.width = "100%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/lbh_cognition.png") ``` ::: Hay un punto importante que es necesario enfatizar aquí: **la regresión múltiple estima el efecto de un predictor después de tener en cuenta el efecto de los demás predictores del modelo**. En otras palabras, los nuevos predictores del modelo tratarán de explicar la variación de los datos que no fue explicada por los otros predictores. Así que **el resultado de la regresión múltiple no es equivalente a los resultados de las regresiones lineales simples** sobre los mismos predictores. Esto puede demostrarse fácilmente corriendo esas regresiones: ```{r} # construir modelos x1y_mod <- lm(formula = y ~ x1, data = xy_datos) x2y_mod <- lm(formula = y ~ x2, data = xy_datos) # atajo a los coeficientes coef(xy_mod) coef(x1y_mod) coef(x2y_mod) ``` Las estimaciones de las mismas variables varían considerablemente entre la regresión múltiple y las regresiones de un solo predictor. Este punto se demuestra además por el hecho de que, si uno de los predictores no tiene ninguna influencia en la respuesta, el efecto del predictor adicional convergerá a su efecto en una regresión lineal simple. Para simular este escenario, fijamos b2 en 0: ```{r} # definir semilla set.seed(123) # numero de observaciones n <- 50 b0 <- -4 b1 <- 3 b2 <- 0 error <- rnorm(n = n, mean = 0, sd = 1) # variables aleatorias x1 <- rnorm(n = n, mean = 0, sd = 1) x2 <- rnorm(n = n, mean = 0, sd = 1) y <- b0 + b1 * x1 + b2 * x2 + error # crear data frame xy_datos <- data.frame(x1, x2, y) # construir modelos xy_mod <- lm(formula = y ~ x1 + x2, data = xy_datos) x1y_mod <- lm(formula = y ~ x1, data = xy_datos) # shortcut to coefficients coef(xy_mod) coef(x1y_mod) ``` El estimado de $\beta_1$ fue casi el mismo en la regresión múltiple (`r coef(xy_mod)[2]`) y en la regresión de un solo predictor (`r coef(x1y_mod)[2]`) Por comodidad, hemos utilizado `coef()` para extraer sólo los estimado de la regresión, pero los valores son los mismos que obtenemos con `summary(model)`. ::: {.alert .alert-info} Ejercicio 5 La siguiente simulación genera un par de predictores continuos altamente colineales: ```{r} # definir semilla set.seed(123) # numero de observaciones n <- 30 b0 <- -4 b1 <- 3 b2 <- -2 error <- rnorm(n = n, mean = 0, sd = 1) # variables aleatorias colineales x1 <- rnorm(n = n, mean = 0, sd = 1) x2 <- x1 + rnorm(n = n, mean = 0, sd = 0.2) # hacer x2 colineal con x1 y <- b0 + b1 * x1 + b2 * x2 + error # crear data frame xy_datos <- data.frame(x1, x2, y) cor(x1, x2) ``` ```{r, eval = FALSE} # graficar ggplot(data = xy_datos, aes(x = x1, y = x2)) + geom_point(size = 3) + ``` ```{r, echo = FALSE} # graficar ggplot(data = xy_datos, aes(x = x1, y = x2)) + geom_point(size = 3) + theme(panel.background = element_rect(fill = "#d9edf7"), plot.background = element_rect(fill = "#d9edf7", colour=NA)) ``` - Construya un modelo lineal multiple (con `lm()`) con 'y' como respuesta y 'x1' y 'x2' como predictores - Compare los betas estimados por el modelo con los usados para generar los datos. ¿Es una buena estimación?. También observe los valores de p. ¿Esperaría que fueran significativos? ::: ------------------------------------------------------------------------ ## Añadir un predictor categórico Para los predictores categóricos podemos crear primero una variable binaria (0, 1) y luego añadir etiquetas a cada valor: ```{r} # definir semilla set.seed(13) # numero de observaciones n <- 50 b0 <- -3 b1 <- 2 error <- rnorm(n = n, mean = 0, sd = 3) # variables aleatorias x1_num <- sample(0:1, size = n, replace = TRUE) y <- b0 + b1 * x1_num + error x1 <- factor(x1_num, labels = c("a", "b")) # crear data frame xy_datos_cat <- data.frame(x1, x1_num, y) head(xy_datos_cat) ``` Y así es como se escribe formalmente: <center>$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1}$</center> Lo mismo que con los predictores continuos. Podemos explorar el patrón de los datos utilizando un diagrama de cajas (boxplot): ```{r} # graficar ggplot(xy_datos_cat, aes(x = x1, y = y)) + geom_boxplot() ``` ... y obtener los estimados del modelo: ```{r} # construir modelos xy_mod_cat <- lm(formula = y ~ x1, data = xy_datos_cat) summary(xy_mod_cat) ``` ... graficar los tamaños de efecto: ```{r} ci_df <- data.frame(param = names(xy_mod_cat$coefficients), est = xy_mod_cat$coefficients, confint(xy_mod_cat)) ggplot(ci_df, aes(x=param, y=est)) + geom_hline(yintercept = 0, color="red", lty = 2) + geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + labs(x = "Parámetro", y = "Tamaño de efecto") + coord_flip() ``` ... y los gráficos diagnósticos del modelo: ```{r, out.width= "80%", fig.height=4} plot_model(xy_mod_cat, type = "diag")[[2]] ``` ::: {.alert .alert-success} <center>Interpretación del modelo</center> Cuadro con los coeficientes (estimados): ```{r, echo= FALSE} summ_xy_mod_cat <- summary(xy_mod_cat) summ_xy_mod_cat$coefficients ``` - El modelo encontró que $\beta_1$ (el tamaño de efecto de 'x1') es `r summ_xy_mod_cat$coefficients[2, 1]` y que es significativamente diferente de 0 (p = `r summ_xy_mod_cat$coefficients[2, 4]`) - Los valores simulados de los parámetros de regresión pueden compararse con el resumen del modelo `lm()` para tener una idea de la precisión del modelo: - $\beta_1$ se fijó en `r b1` y se estimó como `r round(summ_xy_mod_cat$coefficients[2, 1], 3)`. - Tenga en cuenta que en este caso el intercepto se refiere a la estimación del nivel 'a' en el predictor categórico, que se utilizó como base: ```{r, eval=FALSE} # graficar ggplot(xy_datos_cat, aes(x = x1, y = y)) + geom_boxplot() + geom_hline(yintercept = xy_mod_cat$coefficients[1], col = "blue") ``` ```{r, echo=FALSE} ggplot(xy_datos_cat, aes(x = x1, y = y)) + geom_boxplot() + geom_hline(yintercept = xy_mod_cat$coefficients[1], col = "blue") + theme(panel.background = element_rect(fill = "#dff0d8"), plot.background = element_rect(fill = "#dff0d8", colour=NA)) ``` - Por lo tanto, el intercepto es el mismo que la media de *y* para la categoría 'a': ```{r} mean(xy_datos_cat$y[xy_datos_cat$x1 == "a"]) ``` - Observe también que la etiqueta del estimado es 'x1b', no 'x1' como en los predictores continuos. Esto se debe a que en este caso el estimado se refiere a la diferencia entre los dos niveles de la variable categórica ('a' y 'b'). Más concretamente, nos dice que en promedio las observaciones de la categoría 'b' son `r round(summ_xy_mod_cat$coefficients[2, 1], 3)` más altas que las observaciones de la categoría 'a'. ::: ::: {.alert .alert-warning} <center>Caso de estudio</center> - Rico-Guevara A, & M Araya-Salas. 2015. [**Bills as daggers? A test for sexually dimorphic weapons in a lekking hummingbird species**](https://academic.oup.com/beheco/article/26/1/21/2262689). *Behavioral Ecology*. 26 (1): 21-29. "*Los machos con puntas de pico más grandes y puntiagudas tuvieron más éxito en ganar control de territorios en el lek.*" ```{r, out.width = "75%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/daggers_lbh.png") ``` ::: ::: {.alert .alert-info} Ejercicio 6 - Los datos desbalanceados cuando hay categorías (es decir, algunas categorías tienen muchas más observaciones que otras) pueden ser problemáticos para la inferencia estadística. Modifique el código [arriba](#cb99-1) para simular un juego de datos muy desbalanceado y compruebe la precisión del modelo. ::: ------------------------------------------------------------------------ ::: {.alert .alert-warning} <center>Variables indicadoras (Dummy variables)</center> En un modelo de regresión, los predictores categóricos también se representan como vectores numéricos. Más concretamente, los predictores categóricos se codifican como 0s y 1s, en los que 1 significa "pertenece a la misma categoría" y 0 "pertenece a una categoría diferente". Mantuvimos el vector numérico original ('x1_num') al simular el juego de datos con el predictor categórico: ```{r} head(xy_datos_cat) ``` Observe que las "b" de la columna "x1" se convierten en 1 en la columna "x1_num" y las "a" se convierten en 0. Esto se denomina variable indicadora y el proceso se conoce como 'codificación indicadora' (dummy coding). En realidad, podemos utilizar el vector numérico en el modelo de regresión y obtener exactamente los mismos resultados: ```{r} # summary model with categorical variable summary(xy_mod_cat)$coefficients # construir modelos with dummy variable xy_mod_num <- lm(formula = y ~ x1_num, data = xy_datos_cat) # summary with dummy coding summary(xy_mod_num)$coefficients ``` Las cosas se complican un poco más cuando se codifica un predictor categórico con más de dos niveles. Pero la lógica es la misma. ::: ```{r, eval = FALSE, echo = FALSE} This approach can be extended to simulate categorical variables with more than 2 levels: # definir semilla set.seed(123) # numero de observaciones n <- 50 b0 <- -4 b1 <- 3 error <- rnorm(n = n, mean = 0, sd = 1) # variables aleatorias x1 <- rbinom(n = n, size = 2, prob = c(0.33, 0.33)) y <- b0 + b1 * x1 + error x1 <- factor(x1, labels = c("a", "b", "c")) # crear data frame xy_datos <- data.frame(x1, y) # boxplot(formula = y ~ x1, data = xy_datos) ggplot(xy_datos, aes(x = x1, y = y)) + geom_boxplot() + geom_hline(yintercept = xy_mod$coefficients[1], col = "blue")   And these are the results: # construir modelos xy_mod <- lm(formula = y ~ x1, data = xy_datos) summary(xy_mod) ``` ------------------------------------------------------------------------ ## Interacciones Una interacción estadística se refiere a un efecto de una variable predictora que está mediado por una segunda variable. <center>$\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1} + \beta_{2} * x_{2} + \beta_{3} * x_{1} * x_{2}$</center> Esto es más fácil de entender si se observa la interacción de una variable continua y una binaria: ```{r} # definir semilla set.seed(123) # numero de observaciones n <- 50 b0 <- -4 b1 <- 3 b2 <- 1.7 b3 <- -3 error <- rnorm(n = n, mean = 0, sd = 3) # variables aleatorias x1 <- rbinom(n = n, size = 1, prob = 0.5) x2 <- rnorm(n = n, mean = 0, sd = 1) # interaccion se añade como el producto dex1 y x2 y <- b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2 + error x1 <- factor(x1, labels = c("a", "b")) # crear data frame xy_datos_intr <- data.frame(x1, x2, y) head(xy_datos_intr) # construir modelos xy_mod_intr <- lm(formula = y ~ x1 + x2 + x1 * x2, data = xy_datos_intr) # guardar resumen para graficar lineas de mejor ajuste xy_summ_intr <- summary(xy_mod_intr) xy_summ_intr ``` También ayuda graficar los datos: ```{r, eval = FALSE} # graficar ggplot(data = xy_datos_intr, aes(x = x2, y = y, color = x1)) + geom_point(size = 3) + geom_smooth(method = "lm", se = FALSE) ``` ```{r, echo = FALSE, out.width = "80%"} # graficar ggplot(data = xy_datos_intr, aes(x = x2, y = y, color = x1)) + geom_point(size = 3) + scale_color_viridis_d(end = 0.9, begin = 0.2, alpha = 0.6) + ggplot2::geom_smooth(method = "lm", se = FALSE) ``` ... y los tamaños de efecto: ```{r} ci_df <- data.frame(param = names(xy_mod_intr$coefficients), est = xy_mod_intr$coefficients, confint(xy_mod_intr)) ggplot(ci_df, aes(x=param, y=est)) + geom_hline(yintercept = 0, color="red", lty = 2) + geom_pointrange(aes(ymin = X2.5.., ymax = X97.5..)) + labs(x = "Parámetro", y = "Tamaño de efecto") + coord_flip() ``` También deberíamos revisar los gráficos de diagnóstico: ```{r, out.width= "80%", fig.height=4} plot_model(xy_mod_intr, type = "diag") ``` ::: {.alert .alert-success} <center>Interpretación del modelo</center> Cuadro con los coeficientes: ```{r, echo= FALSE} xy_summ_intr$coefficients ``` - El modelo encontró que $\beta_1$ (el tamaño de efecto de 'x1-b' a 'x1-a') es `r xy_summ_intr$coefficients[2, 1]` y que es significativamente diferente de 0 (valor p = `r xy_summ_intr$coefficients[2, 4]`) - El modelo encontró que $\beta_2$ (el tamaño de efecto de 'x2') es `r xy_summ_intr$coefficients[3, 1]` y que es significativamente diferente de 0 (p = `r xy_summ_intr$coefficients[3, 4]`). Esto es en realidad la pendiente de la relación entre x2 e y cuando x1 = \`a' - El modelo encontró que $\beta_3$ (el tamaño de efecto del término de interacción 'x1 \* x2') es `r xy_summ_intr$coefficients[4, 1]` y que es significativamente diferente de 0 (p = `r xy_summ_intr$coefficients[4, 4]`). Esta es la diferencia entre las pendientes de x2 *vs* y cuando x1 = 'a' y x2 *vs* y cuando x1 = 'b'. - Los valores simulados para los parámetros de regresión pueden compararse con el resumen del modelo `lm()` para tener una idea de la precisión del modelo: - $\beta_1$ se simuló con un valor de `r b1` y se estimó en `r round(xy_summ_intr$coefficients[2, 1], 3)` - $\beta_2$ se simuló con un valor de `r b2` y se estimó en `r round(xy_summ_intr$coefficients[3, 1], 3)` - $\beta_3$ se simuló con un valor de `r b3` y se estimó en `r round(xy_summ_intr$coefficients[4, 1], 3)` ::: ------------------------------------------------------------------------ ::: {.alert .alert-warning} <center>Caso de estudio</center> - Chirino F, B Wilink, Araya-Salas M. In prep. **Climatic factors affecting vocal activity in lemur leaf frogs**. "*El aumento de la luz de la luna disminuye la actividad vocal de* Agalychnis lemur *aunque esta relación es mediada por la temperatura.*" ```{r, out.width = "70%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/lemur.png") ``` ::: ::: {.alert .alert-info} Ejercicio 7 - Modifique el [código usado para simular un único predictor asociado](#simular-un-predictor-asociado) aumentando gradualmente el error. Esto se hace aumentando el argumento 'sd' en `error <- rnorm(n = n, sd = 2)`. - Mire cómo los errores más grandes afectan la inferencia (debe volver a correr los modelos) ```{r, eval = FALSE, echo = FALSE} # definir semilla set.seed(123) # numero de observaciones n <- 50 b0 <- -4 b1 <- 3 error <- rnorm(n = n, sd = 10) # variables aleatorias x1 <- rnorm(n = n, mean = 0, sd = 1) y <- b0 + b1 * x1 + error # crear data frame xy_datos <- data.frame(x1, y) # construir modelos summary(lm(formula = y ~ x1, data = xy_datos)) ``` - Ahora reemplace el término de error con `error <- rexp(n = n, rate = 0.2)`. Esto crea un error con una distribución exponencial (por lo tanto no normal). Se supone que esto es problemático para el poder inferencial de estos modelos. Compare las estimaciones que obtuvo con los valores de la simulación ('b0' y 'b1'). Explore la distribución de los residuos (`plot(nombre_del_modelo)`) para los modelos de error 'normal' y 'exponencial'. ```{r, eval = FALSE, echo = FALSE} # definir semilla set.seed(123) # numero de observaciones n <- 50 b0 <- -4 b1 <- 3 error <- rnorm(n = n, sd = 5) sd(error) error <- rexp(n = n, rate = 0.2) sd(error) # variables aleatorias x1 <- rnorm(n = n, mean = 0, sd = 1) y <- b0 + b1 * x1 + error # crear data frame xy_datos <- data.frame(x1, y) # construir modelos mod <- lm(formula = y ~ x1, data = xy_datos) coef(mod) resd_mod <- resid(mod) # ver distribucion ggplot(data = data.frame(resd_mod), mapping = aes(x = resd_mod)) + geom_histogram() ``` - Se supone que la colinealidad (la presencia de predictores correlacionados) afecta a la estabilidad de la regresión múltiple. El siguiente código crea dos predictores altamente colineales ('x1' y 'x2'). La última línea de código muestra la correlación entre ellos. ```{r} # definir semilla set.seed(123) # numero de observaciones n <- 50 b0 <- -4 b1 <- 3 b2 <- -2 error <- rnorm(n = n, mean = 0, sd = 1) # variables aleatorias x1 <- rnorm(n = n, mean = 0, sd = 1) # hacer x2 muy parecido a x1 x2 <- x1 + rnorm(n = n, mean = 0, sd = 0.3) cor(x1, x2) ``` ```{r, eval=FALSE, echo=FALSE} xy_datos <- data.frame(x1, x2, y) # construir modelos xy_mod <- lm(formula = y ~ x1 + x2, data = xy_datos) summary(xy_mod) ``` - Construya un modelo de regresión múltiple para estos datos (y \~ x1 + x2). Puede utilizar el mismo código que en la sección [Añadir más de un predictor: regresión múltiple](#añadir-más-de-un-predictor:-regresión-múltiple). - ¿Cómo afecta a la inferencia la presencia de predictores colineales? Realice los gráficos de diagnóstico para este modelo (`plot(nombre_del_modelo)`). - Simule un juego de datos con tres predictores en los que sólo dos de ellos son altamente colineales. Ajuste un modelo de regresión múltiple (y \~ x1 + x2 + x3) para esos datos y observe cómo la colinealidad afecta a la estimación del predictor no colineal. ::: ------------------------------------------------------------------------ # Extender los modelos lineales a estructuras de datos más complejas ## Modelos lineales generalizados (GLM) Los GLM nos permiten modelar la asociación a variables respuesta que no siguen una distribución normal. Además, permiten modelar distribuciones que se asemejan más al proceso que generó los datos. El siguiente código crea un juego de datos con una respuesta que representa cuentas (por lo tanto, no normal): ```{r} # definir semilla set.seed(1234) # numero de muestas n <- 50 # coeficientes b0 <- 1.2 b1 <- 1.3 b2 <- 0 #generar variables y <- rpois(n = n, lambda = 6.5) # lambda = tasa promedio de exitos x2 <- seq(-0.5, 0.5, , length(y)) x1 <- (log(y) - b0 - b2 * x2) / b1 # crear data frame xy_datos_pois <- data.frame(x1, x2, y) head(xy_datos_pois) ``` También grafiquemos 'x1' *vs* 'y': ```{r} # graficar ggplot(xy_datos_pois, aes(x = x1, y = y)) + geom_point() ``` La relación no parece muy lineal ni la varianza parece ser constante a través de 'x1'. Podemos relajar el requisito de la distribución normal con GLMs. `glm()` es una función básica de R que nos ayuda a hacer el truco. Para este ejemplo la distribución más apropiada es *Poisson*. Esto se puede establecer en el argumento \`family' así: ```{r} glm_pois <- glm(formula = y ~ x1 + x2, data = xy_datos_pois, family = poisson()) ``` Que es equivalente a esto: <center>$\hat{Y} \sim_{poisson} \beta_{o} + \beta_{1} * x_{1} + \beta_{2} * x_{2}$</center> Como puede ver el único argumento extra comparado con `lm()` es `family`. El resto es simplemente la "fórmula" y los "datos" con los que ya estamos familiarizados. Así que, de nuevo, podemos aprovechar nuestros conocimientos sobre los modelos lineales para extenderlos a estructuras de datos más complejas. También necesitamos ejecutar `summary()` para obtener el resultado del modelo: ```{r} summary(glm_pois) ``` ::: {.alert .alert-success} <center>Interpretación del modelo</center> Cuadro de los coeficientes (estimados): ```{r, echo= FALSE} summ_glm_pois <- summary(glm_pois) summ_glm_pois$coefficients ``` - El modelo nos dice que $\beta_1$ (el tamaño de efecto de 'x1') es `r summ_glm_pois$coefficients[2, 1]` y que es significativamente diferente de 0 (p = `r summ_glm_pois$coefficients[2, 4]`). Esto se interpreta en realidad como que un aumento de 1 unidad de 'x1' resulta en 'y' (tasa) en un factor de exp(`r summ_glm_pois$coefficients[2, 1]`) = `r exp(summ_glm_pois$coefficients[2, 1])`. - El modelo también nos dice que $\beta_2$ (el tamaño de efecto de 'x2') es `r summ_glm_pois$coefficients[3, 1]` y que es significativamente diferente de 0 (p = `r summ_glm_pois$coefficients[3, 4]`). Esto significa que un aumento en 1 unidad de 'x2' resulta en 'y' (tasa) en un factor de exp(`r summ_glm_pois$coefficients[3, 1]`) = `r exp(summ_glm_pois$coefficients[3, 1])`. ::: ::: {.alert .alert-warning} <center>Caso de estudio</center> - Topp, E. N., Tscharntke, T., & Loos, J. (2022). [**Fire and landscape context shape plant and butterfly diversity in a South African shrubland**](https://onlinelibrary.wiley.com/doi/pdf/10.1111/ddi.13257). *Diversity and Distributions*, 28(3), 357-371. "*La riqueza de especies de mariposas es de tres a cuatro veces mayor cuando aumenta el hábitat natural en el paisaje circundante (en un radio de 2 km), mientras que la abundancia de mariposas se asociaba negativamente con el aumento del tiempo transcurrido desde el último incendio.*" ```{r, out.width = "100%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/butterflies.png") ``` ::: ::: {.alert .alert-info} Ejercicio 8 - Intente correr un modelo `lm()` (con una distribución gaussiana), compare los resultados y compruebe los residuos (`plot_model(model_name, type = "diag")`) ::: Existen muchas otras funciones de distribución y enlace: ```{r, out.width = "80%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/link-functions.jpg") ``` ------------------------------------------------------------------------ ## Modelos de efectos variables (modelos mixtos) A veces nuestros conjuntos de datos incluyen estructuras con varios niveles de organización. Por ejemplo, cuando tomamos muestras de varios individuos de diferentes poblaciones. En esos casos, la variación en el nivel estructural superior (poblaciones) podría impedir la detección de patrones en el nivel inferior (individuos). Vamos a simular algunos datos que se asemejan a ese escenario. Tenemos dos predictores continuos (x1) y una respuesta continua (y). Cada muestra procede de 1 de 8 poblaciones diferentes (poblacion): ```{r} # x<- 1 # definir semilla set.seed(28) # numero de observaciones n <- 300 b0 <- 1 b1 <- 1.3 poblacion <- sample(0:8, size = n, replace = TRUE) error <- rnorm(n = n, mean = 0, sd = 2) # variables aleatorias x1 <- rnorm(n = n, mean = 0, sd = 1) y <- b0 + poblacion * 2 + b1 * x1 + error # add letters poblacion <- letters[poblacion + 1] # create data set xy_datos_mixto <- data.frame(x1, y, poblacion) head(xy_datos_mixto, 10) ``` Podemos explorar la relación entre 'y' y 'x1' con un gráfico: ```{r} ggplot(data = xy_datos_mixto, aes(x = x1, y = y)) + geom_point() ``` ¿Puede ver claramente el patrón de asociación entre las dos variables que hemos utilizado para simular los datos? Podemos seguir explorando los datos con un modelo de regresión lineal simple: ```{r} summary(lm(y ~ x1, data = xy_datos_mixto)) ``` A pesar de haber simulado un $\beta_1$ distinto de cero, no tenemos ninguna asociación significativa según este modelo y la estimación de $\beta_1$ está muy lejos de la simulada. Esta pobre inferencia se debe a que estamos ignorando una característica importante de nuestros datos, la agrupación de las muestras en "poblaciones". Los modelos de efectos mixtos (también conocidos como modelos multinivel o modelos de efectos variables) pueden ayudarnos a tener en cuenta estas características adicionales, mejorando significativamente nuestro poder de inferencia. Coloreemos cada una de las poblaciones para ver cómo covarían las variables en cada subgrupo de datos: ```{r} ggplot(data = xy_datos_mixto, aes(x = x1, y = y, color = poblacion)) + geom_point() ``` Parece haber un patrón claro de asociación positiva entre x1 e y. El patrón se hace un poco más evidente si mostramos cada población en su propio panel: ```{r, out.width= "100%"} ggplot(data = xy_datos_mixto, aes(x = x1, y = y, color = poblacion)) + geom_point() + facet_wrap( ~ poblacion) + geom_smooth(method = "lm", se = FALSE) ``` Construyamos un modelo de efectos mixtos utilizando la población como intercepto variable: ```{r} mix_eff_mod <- lmer(formula = y ~ x1 + (1 | poblacion), data = xy_datos_mixto) ``` Que es equivalente a esto: <center>$\hat{Y} \sim_{gausiana} (\beta_{o} + \beta_{grupo}) + \beta_{1} * x_{1}$</center> Podemos ver el resultado del modelo igual que con un modelo `lm()`: ```{r} summary(mix_eff_mod) ``` El modelo detectó correctamente el patrón simulado y la estimación de $\beta_1$ (`r round(fixef(mix_eff_mod)[2], 3)`) es muy cercana al valor simulado. ::: {.alert .alert-warning} <center>Caso de estudio</center> - Araya-Salas M, P González-Gómez, K Wojczulanis-Jakubas, V López III & T Wright. 2018. [**Spatial memory is as important as weapon and body size for territorial ownership in a lekking hummingbird**](https://www.nature.com/articles/s41598-018-20441-x). *Scientific Reports*. 13, e0189969 "*... el tamaño corporal mostró una correlación negativa con la frecuencia baja del canto.*" ```{r, out.width = "80%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/lbh_song_size.png") ``` ::: ::: {.alert .alert-success} <center>Interpretación del modelo</center> Cuadro con los coeficientes (estimados): ```{r, echo= FALSE} summ_mix_eff_mod <- summary(mix_eff_mod) summ_mix_eff_mod$coefficients ``` - El modelo encontró que $\beta_1$ (el tamaño de efecto de 'x1') es `r summ_mix_eff_mod$coefficients[2, 1]` y que es significativamente diferente de 0 (p = `r summ_mix_eff_mod$coefficients[2, 5]`) - Los valores simulados de los parámetros de regresión pueden compararse con el resumen del modelo `lmer()` para tener una idea de la precisión del modelo: - $\beta_1$ se fijó en `r b1` y se estimó como `r round(summ_mix_eff_mod$coefficients[2, 1], 3)`. - La varianza del intercepto para cada población fue de `r var(ranef(mix_eff_mod)$poblacion[,1])` y la desviación estándar de `r sd(ranef(mix_eff_mod)$poblacion[,1])` - Los interceptos estimados para cada población son: ```{r} ranef(mix_eff_mod)$poblacion ``` - Podemos compararlos con los promedios para cada población en los datos: ```{r} # sumar interceptos por poblacion a el intercepto total interceptos <- ranef(mix_eff_mod)$poblacion + fixef(mix_eff_mod)[1] # calcular los promedios por poblacion en los datos prom_pobs <- aggregate(y ~ poblacion, xy_datos_mixto, mean) # unir cbind(interceptos, prom_pobs= prom_pobs[,2]) ``` - La varianza entre los niveles del factor aleatorio (`poblacion`) es de `r var(ranef(mix_eff_mod)$poblacion[,1])` y la desviación estándar de `r sd(ranef(mix_eff_mod)$poblacion[,1])` ::: ### La paradoja de Simpson La paradoja de Simpson es un fenómeno en estadística en el que una tendencia aparece en varios grupos de datos pero desaparece o se invierte cuando se combinan los grupos. Se ha utilizado para ilustrar el tipo de resultados engañosos que puede generar el ignorar la estructura en los datos y las relaciones causales entre variables. ```{r, out.width = "80%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/simpsons_paradox.gif") ``` Tomado de wikipedia La estructura de los datos debido a observaciones pertenecientes al mismo grupo puede ser tomada en cuenta usando el grupo como un factor aleatorio, tal y como se hizo en el ejemplo anterior. ::: {.alert .alert-info} <center>Modelos lineales como abordaje estadístico estándar </center> ```{r, out.width = "40%", echo = FALSE, fig.align= "center"} knitr::include_graphics("./images/swiss_knife.png") ``` <center>$\hat{Y} \sim_{función\ de\ enlace} (\beta_{o} + \beta_{grupo}) + \beta_{1} * x_{1} + \cdots + \beta_{n} * x_{n}$</center> - Este modelo básico incluye una función de enlace que cuando es gausiana (i.e. normal) equivale a un modelo lineal y cuando no, a un modelo generalizado - Cuando el $\beta_{grupo}$ es diferente entre grupos (i.e. cuando el intercepto difiere entre grupos) es un modelo mixto con intercepto variable (i.e. aleatorio) ::: Este manual solo pretende sugerir un abordaje estadístico centrado en los modelos de regresión desde el cual se pueden llevar a cabo adecuadamente la mayoría de los análisis estadísticos que usamos en biología. **El manual no tiene como objetivo explicar en detalle las particularidades de los modelos generalizados y mixtos**. Note que estos modelos, principalmente los mixtos, pueden adaptarse a otras estructuras de datos mas complejas no cubiertas aquí, como lo son las pendientes aleatorias, pendientes e interceptos aleatorios, estructuras de correlación (i.e. autocorrelación espacial o temporal, pedigrís o árboles filogenéticos), medidas repetidas y el uso de pseudoreplicas, entre otras. ------------------------------------------------------------------------ ## Referencias - [Richard McElreath's Statistical Rethinking book](https://github.com/rmcelreath/stat_rethinking_2022) - [R's rbinom -- Simulate Binomial or Bernoulli trials](https://www.programmingr.com/examples/neat-tricks/sample-r-function/r-rbinom/) - [R's rnorm -- selecting values from a normal distribution](https://www.programmingr.com/examples/neat-tricks/sample-r-function/r-rnorm/) - [R's exp -- Simulating Exponential Distributions](https://www.programmingr.com/examples/neat-tricks/sample-r-function/rexp/) - [Simulating data in R](https://aosmith.rbind.io/2018/08/29/getting-started-simulating-data/) ------------------------------------------------------------------------ Información de la sesión ```{r session info, echo=F} sessionInfo() ```