Capítulo 5 Modelos en R

5.1 Paquetes necesarios para este capítulo

Para este capítulo necesitas tener instalado el paquete tidyverse, broom y MuMIn.

En este capítulo se explicará como generar modelos en R, el como obtener información y tablas a partir de los modelos con el paquete Broom (Robinson and Hayes 2018) y una leve introducción a la selección de modelos con el paquete MuMIn (Barton 2018)

Dado que este libro es un apoyo para el curso BIO4022, esta clase puede también ser seguida en este link. El video de la clase se encontrará disponible en este link.

5.2 Modelos estadísticos

Un modelo estadístico intenta explicar las causas de un suceso basado en un muestreo de la población total. El supuesto es que si la muestra que obtenemos de la población es representativa de esta, podremos inferir las causas de la variación de la población midiendo variables explicativas. En general tenemos una variable respuesta (fenómeno que queremos explicar), y una o varias variables explicativas que generarían deterministamente parte de la variabilidad en la variable respuesta.

5.2.1 Ejemplo

Tomemos el ejemplo de la base de datos CO2 presente en R (Potvin, Lechowicz, and Tardif 1990). Supongamos que nos interesa saber que factores afectan la captación de \(CO_2\) en las plantas.

Tabla 5.1: Primeras 20 variables de la base de datos CO2.
Plant Type Treatment conc uptake
Qn1 Quebec nonchilled 95 16.0
Qn1 Quebec nonchilled 175 30.4
Qn1 Quebec nonchilled 250 34.8
Qn1 Quebec nonchilled 350 37.2
Qn1 Quebec nonchilled 500 35.3
Qn1 Quebec nonchilled 675 39.2
Qn1 Quebec nonchilled 1000 39.7
Qn2 Quebec nonchilled 95 13.6
Qn2 Quebec nonchilled 175 27.3
Qn2 Quebec nonchilled 250 37.1
Qn2 Quebec nonchilled 350 41.8
Qn2 Quebec nonchilled 500 40.6
Qn2 Quebec nonchilled 675 41.4
Qn2 Quebec nonchilled 1000 44.3
Qn3 Quebec nonchilled 95 16.2
Qn3 Quebec nonchilled 175 32.4
Qn3 Quebec nonchilled 250 40.3
Qn3 Quebec nonchilled 350 42.1
Qn3 Quebec nonchilled 500 42.9
Qn3 Quebec nonchilled 675 43.9

En la tabla 5.1 vemos las primeras 20 observaciones de esta base de datos. Vemos que dentro de los factores que tenemos para explicar la captación de \(CO_2\) estan:

  • Type: Subespecie de la planta (Missisipi o Quebec)
  • Treatment: Tratamiento de la plnata, enfriado (chilled) o no enfriado (nonchilled)
  • conc: Concentración ambiental de \(CO_2\), en mL/L.

Una posible explicación que nos permitiría intentar explicar este fenómeno, es que las plantas de distintas subespecies, tendrán distinta captación de \(CO_2\), lo cual exlploramos en el gráfico 5.1:

Captación de CO2 por plantas dependiente de su subespecie

Figura 5.1: Captación de CO2 por plantas dependiente de su subespecie

Vemos que se observa una tendencia a que las plantas con origen en Quebec capten más \(CO_2\) que las que estan en el Mississippi, pero ¿Podemos decir efectivamente que ambas poblaciónes tienen medias distintas medias? Es ahí donde entran los modelos.

5.2.2 Representando un modelo en R

En R la mayoría de los modelos se representan con el siguiente codigo:

En este modelo, tenemos la variable respuesta Y, la cual puede estar explcada por una o multiples variables explicativas X, es por esto que el simbolo ~ se lee explicado por, donde lo que esta a su izquerada es la variable respuesta y a la derecha la variable explicativa. Los datos se encuentran en un data frame y finalmente usaremos alguna función, que identificará algún modelo. Algunas de estas funciones las encontramos en la tabla 5.2

Tabla 5.2: Algunos modelos que podemos generar en R
Modelos Funcion
Prueba de t t.test()
ANOVA aov()
Modelo lineal simple lm()
modelo lineal generalizado glm()
Modelo aditivo gam()
Modelo no lineal nls()
modelos lineales mixtos lmer()
Boosted regression trees gbm()

5.2.3 Volvamos al ejemplo de las plantas

Para este ejemplo usaremos un modelo lineal simple, para esto siguiendo la tabla 5.2 usaremos la función lm:

5.2.3.1 usando broom para sacarle más a tu modelo

El paquete broom (Robinson and Hayes 2018) es un paquete adyacente al tidyverse (por lo que debes cargarlo aparte del tidyverse), el cual nos permite tomar información de modelos generados en formato tidy. Hoy veremos 3 funciones de broom, estas son glance, tidy y augment.

5.2.3.1.1 glance

la función glance, nos entregará información general del modelo, como el valor de p, el \(R^2\), log-likelihood, grados de libertad, y/o otros parametros dependiendo del modelo a utilizar. Esta información nos es entregada en un formato de data frame, como vemos en el código siguiente y en la tabla 5.3

Tabla 5.3: Información del modelo fi1 entregada por la función glance
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.346713 0.3387461 8.794012 43.5191 0 1 -300.8007 607.6014 614.8939 6341.441 82 84
5.2.3.1.2 tidy

la función tidy, nos entregará información sobre los parametros del modelo, esto es el intercepto, la pendiente y/o interacciones, como vemos en el código siguiente y en la tabla 5.4

Tabla 5.4: Información del modelo fi1 entregada por la función glance
term estimate std.error statistic p.value
(Intercept) 33.54286 1.356945 24.719384 0
TypeMississippi -12.65952 1.919011 -6.596901 0
5.2.3.1.3 augment

la función augment, nos entregará para cada observación de nuestro modelo, varios parametros importantes como el valor predicho, los residuales, el distancia de cook entre otros, esto nos sirve principalmente para estudiar los supuestos de nuestro modelo. A continuación vemos el uso de la función augment y 20 de sus observaciones en la tabla 5.5

Tabla 5.5: Información del modelo fi1 entregada por la función augment
uptake Type .fitted .resid .std.resid .hat .sigma .cooksd
32.4 Mississippi 20.88333 11.5166667 1.3254778 0.0238095 8.752828 0.0214255
12.0 Mississippi 20.88333 -8.8833333 -1.0224018 0.0238095 8.791552 0.0127476
37.5 Quebec 33.54286 3.9571429 0.4554360 0.0238095 8.836932 0.0025295
30.6 Mississippi 20.88333 9.7166667 1.1183119 0.0238095 8.780397 0.0152515
31.1 Mississippi 20.88333 10.2166667 1.1758580 0.0238095 8.773216 0.0168615
30.3 Quebec 33.54286 -3.2428571 -0.3732274 0.0238095 8.840611 0.0016988
21.9 Mississippi 20.88333 1.0166667 0.1170103 0.0238095 8.847391 0.0001670
10.5 Mississippi 20.88333 -10.3833333 -1.1950400 0.0238095 8.770741 0.0174161
7.7 Mississippi 20.88333 -13.1833333 -1.5172980 0.0238095 8.723037 0.0280755
31.8 Mississippi 20.88333 10.9166667 1.2564225 0.0238095 8.762547 0.0192512
30.9 Mississippi 20.88333 10.0166667 1.1528396 0.0238095 8.776132 0.0162078
18.0 Mississippi 20.88333 -2.8833333 -0.3318490 0.0238095 8.842186 0.0013430
30.0 Mississippi 20.88333 9.1166667 1.0492567 0.0238095 8.788531 0.0134261
10.6 Mississippi 20.88333 -10.2833333 -1.1835308 0.0238095 8.772231 0.0170823
22.0 Mississippi 20.88333 1.1166667 0.1285196 0.0238095 8.847238 0.0002014
41.4 Quebec 33.54286 7.8571429 0.9042954 0.0238095 8.803900 0.0099726
13.7 Mississippi 20.88333 -7.1833333 -0.8267452 0.0238095 8.811176 0.0083355
39.6 Quebec 33.54286 6.0571429 0.6971295 0.0238095 8.821870 0.0059267
16.2 Quebec 33.54286 -17.3428571 -1.9960265 0.0238095 8.630502 0.0485869
19.9 Mississippi 20.88333 -0.9833333 -0.1131739 0.0238095 8.847439 0.0001562

5.2.3.2 Selección de modelos usando broom y el AIC

El AIC, o Criterio de informacion de Akaike (Aho, Derryberry, and Peterson 2014), es una medida de cuanta información nos entrega un modelo dada su complejidad. Esta última medida a partir del número de parámetros que tiene. Cuanto más bajo sea el AIC, mejor comparativamente es un modelo, y en general, un modelo que sea dos unidades de AIC menor que otro modelo, será considerado un modelo que es significativamente mejor que otro.

La formula del criterio de selección de Akaike es la que vemos en la ecuación (5.1).

\[\begin{equation} AIC = 2 K - 2 \ln{(\hat{L})} \tag{5.1} \end{equation}\]

Donde \(K\) es el número de parametros, lo cual podemos ver con tidy, si vemos la tabla 5.4, vemos que el modelo Fit1 tiene 2 parametros, esto es \(K\) es igual a 2.

El log-likelihood del modelo (\(\ln{(\hat{L})}\)) es el ajuste que este tiene a los datos. Cuanto más positivo es este valor mejor se ajusta el modelo a los datos, y cuanto mas negativo es, menos se ajusta a los datos, en nuestro modelo, usando glance, podemos ver que el valor del log-likelyhood del modelo es de -300.8 (ver tabla 5.4).

Por lo tanto remplazando la ecuación (5.1), obtenemos 605.6, que es un valor muy cercano a los 608, que aparecen en el glance del modelo (tabla 5.4).

5.2.3.2.1 Modelos candidatos

Veamos la figura 5.2. para pensar cuales podrían ser modelos interesantes a explorar.

Gráfico exploratorio para generar modelos de la base de datos CO2

Figura 5.2: Gráfico exploratorio para generar modelos de la base de datos CO2

Referencias

Aho, Ken, DeWayne Derryberry, and Teri Peterson. 2014. “Model Selection for Ecologists: The Worldviews of Aic and Bic.” Ecology 95 (3). Wiley Online Library: 631–36.

Barton, Kamil. 2018. MuMIn: Multi-Model Inference. https://CRAN.R-project.org/package=MuMIn.

Potvin, Catherine, Martin J Lechowicz, and Serge Tardif. 1990. “The Statistical Analysis of Ecophysiological Response Curves Obtained from Experiments Involving Repeated Measures.” Ecology 71 (4). Wiley Online Library: 1389–1400.

Robinson, David, and Alex Hayes. 2018. Broom: Convert Statistical Analysis Objects into Tidy Tibbles. https://CRAN.R-project.org/package=broom.