Introducción

El propósito de este tutorial es introducir la implementación de modelos de ANOVA y Regresión Lineal en R conjuntamente, para evidenciar que básicamente son lo mismo (con algunas reservas). ¿Cuáles reservas?. Cuando nos referimos vagamente a modelos de Regresión Lineal usualmente deberíamos utilizar el término modelo Lineal General (glm, ojo, no confundir con GLM - o Modelos Lineales Generalizados). El modelo lineal general (glm) incluye una familia de modelos estadísticos diferentes: ANOVA, ANCOVA, MANOVA, MANCOVA, regresión lineal ordinaria, prueba \(t\) y prueba \(\mathcal{F}\). Estrictamente, el glm es una generalización de la regresión lineal múltiple al caso de más de una variable dependiente. Pero, volviendo al foco del tutorial, ¿en qué se diferencian el ANOVA, la regresión lineal ordinaria y las demás técnicas de la familia?. La regresión es un glm que incluye covariables continuas. Hablamos de análisis de covarianza (ANCOVA) cuando las covariables continuas se combinan con variables ficticias (“dummy”, i.e. variables utilizadas habitualmente para representar variables categóricas que tienen más de dos niveles - llamadas “factores” en un ANOVA). Si sólo se utilizan variables ficticias, nos referimos a análisis de varianza (ANOVA). MANOVA y MANCOVA constituyen las versiones multivariadas del ANOVA y ANCOVA, respectivamente, que consideran más de una variable respuesta (o dependiente) simultáneamente (pero las dejaremos de lado, por ahora, para enfocarnos en este tutorial en la comparación entre ANOVA y Regresión Lineal). Como veremos, el ANOVA se distingue además porque en el procedimiento para probar coeficientes significativos utiliza la descomposición de la varianza en componentes del término del modelo (Cuadrado Medio del Modelo) y el componente del término de error (Cuadrado Medio del Error). Un ANOVA pregunta a través de una prueba de hipótesis de \(\mathcal{F}\) sobre un parámetro si la relación entre la varianza explicada y la varianza no explicada es suficientemente alta. En tanto, en una regresión lineal una prueba \(t\) sobre el parámetro de un modelo de regresión pregunta si el efecto del parámetro es suficientemente distinto de cero, obteniéndose además la dirección del cambio. Aunque casi todo los modelos de ANOVA pueden implementarse como regresión lineal, el ANOVA es útil para entender y estructurar diseños complejos (como los modelos multinivel).

Modelos lineales

Los datos usualmente son analizados para crear modelos que expliquen cómo están organizados los elementos de los datos y cómo se relacionan entre sí y con las propiedades de las entidades del mundo real. Generalmente asumimos que las observaciones, que conforman los datos, tienen un componente previsible y otro aleatorio (no previsible). La parte previsible puede expresarse matemáticamente como una función de parámetros desconocidos, mientras que para la parte aleatoria asumimos una distribución probabilística determinada. Entonces buscamos estimar los parámetros desconocidos basándonos en las muestras. El Análisis de Varianza (ANOVA) y la Regresión Lineal son dos modelos aditivos con la forma:

Observación = Previsible + Aleatorio

O, más formalmente:

\[Y_{i} = \theta + \varepsilon_{i}\]

donde, por ejemplo:
\(Y_i\) es la producción de masas de huevos (puestas) de la i-ésima lapita pulmonada (\(Siphonaria\) sp.) adulta del intermareal rocoso;
\(\theta\) es el efecto fijo, común a todos las lapita (por ejemplo, la densidad de lapitas adultas); y
\(\varepsilon_{ijk}\): es el efecto residual no controlado (error aleatorio), que también contiene la información de variables no consideradas en el modelo (como, por ejemplo, la estación del año, etc).
Para simplificar se imponen como condiciones al modelo que los residuos tomen valores positivos o negativos y que su media, calculada sobre todos los individuos de la población (no solamente los de la muestra), sea 0. Asumimos entonces que:

\[E(\varepsilon) = 0\], y

\[Var(\varepsilon) = \sigma_{\varepsilon}^{2}\]

Entonces, para entender el comportamiento del número de puestas necesitaremos estimar los parámetros del modelo \(\varepsilon\) y \(\sigma_{\varepsilon}^{2}\) mediante el método de cuadrados mínimos, que minimiza la suma de los cuadrados residuales:

\[\varepsilon =Y_{i} - \theta, \ \forall \ i = 1, 2,\ …n\]

Esto significa que el mejor valor de \(\theta\) será aquel que produzca residuos pequeños para las \(n\) observaciones de la muestra. Consecuentemente se define la Suma de Cuadrados de los residuos, que se minimiza (o sea, que es mínima) para el mejor valor de \(\theta\). Entonces, el objetivo del método es encontrar la solución \(\hat\theta\) que minimiza la \(SC(\theta)\). En otras palabras, la media muestral es el valor que minimiza la suma de los cuadrados residuales. Concluyendo, el modelo:

\[Y_{i} = \theta + \varepsilon_{i},\]

que describe el número de puestas por lapita pulmonada, sugiere el valor \(\hat\theta = \bar y\) como el mejor estimador puntual (i.e. la media) y \(s_{\varepsilon}^{2}\) como la mejor medida de variabilidad (la varianza). En tanto, \(SC(\hat\theta)\) representa el cuadrado de la cantidad de información no explicada por el modelo.

Sin embargo, para hacer inferencias sobre la población (a partir de la muestra) necesitamos asumir los siguientes supuestos adicionales:

El error muestral tiene una distribución normal con media 0 y varianza \(\sigma_{\varepsilon}^{2}\): \[\varepsilon \sim N (0, \sigma_{\varepsilon}^{2})\]

Los errores asociados con cualquier par de observaciones son independientes: \[E(\varepsilon_i \varepsilon_j = 0, \ \forall \ i \neq j)\]

Estos dos supuestos sobre la distribución de los errores explica por qué las pruebas gráficas y de hipótesis se hacen sobre los residuales. En los diagnósticos gráficos esperamos ver que los residuos se distribuyan en un plano cartesiano como una “nube de puntos”, sin tendencias o patrones que sugieran una distribución no aleatoria o falta de independencia (explicada por una autocorrelación espacial o temporal de las obseervaciones, por ejemplo).

Hasta acá definimos a nuestro modelo para UNA población:

\[Y_{i} = \theta + \varepsilon_{i},\]

Noten que \(\theta\) no varía, ya que corresponde a un efecto fijo común a todos las lapitas (que en nuestro ejemplo es una estación del año - digamos, “verano” - o una densidad, 30 lapitas por 225 \(cm^{2}\)). Por eso hablamos de que nuestro modelo está definido para sólo una población. Ahora, si permitimos que \(\theta\) pueda tomar dos “valores” diferentes, por ejemplo, “verano” y “primavera”, nuestro modelo es:

\[Y_{i,j} = \theta_{i} + \varepsilon_{i,j}, \ i = 1,2, \ j=1, 2, …, n_i\]

donde:
\(Y_{i,j}\) es la variable respuesta del j-ésimo individuo en el \(i\)-ésimo grupo. O sea, el número de masas de huevos de la \(j\)-ésima lapita pulmonada adulta del intermareal rocoso, en verano (o primavera);
\(\theta\) es el efecto fijo, común a todos las lapitas del grupo \(i\) (con \(i=1, 2, …, n_i\). En el ejemplo, supongamos \(i = 1 = primavera, i = 2 = verano\) una estación del año); y
\(\varepsilon_{ij}\): es el efecto residual no controlado (error aleatorio) del \(j\)-ésimo individuo en el \(i\)-ésimo grupo.

Análogamente al caso de una población,necesitaremos estimar los parámetros \(\theta_{1}\) y \(\theta_{2}\). Pero además, ahora querremos verificar si ambos parámetros difieren. En el ejemplo, si conocer la estación reproductiva de las lapitas nos ayuda a predecir el número de puestas de una lapita.

Similarmente al caso de una población, para hacer inferencias sobre las poblaciiones (a partir de la muestra) necesitamos asumir una distribución probabilística para la parte aleatoria (los errores):

El error muestral tiene una distribución normal con media 0 y varianza \(\sigma_{\varepsilon}^{2}\), tal que:

\[\varepsilon_{i,j} \sim N (0, \sigma_{\varepsilon}^{2}), \ \forall i = 1,2\ \ \text{y} \ \ j = 1, 2, …, n_{i}\]

Existe independencia de las observaciones dentro de cada población. Esto es, los errores asociados con cualquier par de observaciones dentro de un grupo son independientes: \[E(\varepsilon_{i,j} \varepsilon_{j,k}) = 0, \ \forall\ j \neq k \ \ \text{e} \ \ i = 1, 2\]

Existe independencia de las observaciones entre las poblaciones. Esto es, los errores asociados con cualquier par de observaciones pertenecientes a grupos distintos son independientes: \[E(\varepsilon_{1,j} \varepsilon_{j,2}) = 0, \ \forall \ j\ \ \text{y}\ \ k\]

El modelo para dos poblaciones puede generalizarse a más poblaciones, como:

\[Y_{i,j} = \theta_{i} + \varepsilon_{i,j}, \ \text{con}\ \ i=1,2,…,k \ \ \text{y}\ \ j=1,2,…,n_i\]

donde \(k\) es el número de poblaciones, y \(n\) es el número de observaciones individuales en cada población \(i\).

ANOVA y RL: Teoría

A partir de acá haremos explícitas las diferencias de notación entre modelos de ANOVA y de Regresión Lineal (RL). Recordamos que en ambos tipos de modelos las variables dependientes son cuantitativas. Se diferencian, sin embargo, en el tipo de variables indepedientes consideradas. Los modelos de ANOVA se utilizan cuando las variables independientes son categóricas. Se denominan factores, y cada factor se representa por diferentes categorías o niveles, como “Alta”, “Intermedia” y “Baja” para el factor “Densidad”, o niveles “Control”, “Ambiente” y “Enriquecido” de “Concentraciones de nitrógeno inorgánico disuelto”). En cambio, los modelos de RL consideran variables cuantitativas (no cualitativas, o categóricas), de “razón” (para las que el valor 0 no es arbitrario, i.e. indica ausencia de la variable) como, por ejemplo, la densidad o la concentración de nitrógeno inorgánico disuelto, o de “intervalo” (para las que el valor 0 es arbitrario, como la escala de medición de la temperatura en grados Celsius; no así en la escala Kelvin).

La notación más familiar de un modelo de ANOVA de 1 vía (o, equivalentemente, de un factor) es:

\[y_{i,j,k} = \mu + \alpha_{i} + \varepsilon_{i,j,k}, \ \text{con}\ \ i=1,2,…,l \ \ \text{y}\ \ j=1,2,…,m\]

donde:
\(y_{i,j,k}\) es el valor de la variable respuesta de la k-ésima observación en el tratamiento \(i\) (notar que, en el contexto de un factor, un tratamiento es equivalente a un nivel del factor considerado en el modelo);
\(\mu\) es el parrámetro que representa la media general de todos los tratamientos. \(\alpha_{i}\) es el parámetro que representa el efecto fijo del nivel \(i\) del factor \(\alpha\) (con \(i=1, 2, …, l\); y
\(\varepsilon_{i,k}\): es una variable aleatoria que refleja los efectos no controlados que actúan sobre las unidades experimentales. Residual o error aleatorio de la \(k\)-ésima observación en el \(i\)-ésimo tratamiento.

Noten que los diferentes \(k\) \(\alpha_{i}\)s representan las \(k\) poblaciones, que son los “niveles” del “factor” \(\alpha\).

Comparativamente, en su notación presentada más frecuentemente, el modelo análogo de Regresión Lineal (RL) con un predictor (o variable independiente) es:

\[Y_{i} = \beta_{0} + \beta_{1} * X_{1i} + \varepsilon_{i}, \ \text{con}\ \ i=1,2,…,k\]

donde:
\(Y_{i}\) es la variable respuesta del \(i\)-ésimo individuo;
\(\beta_{0}\) representa al parámetro denominadoo “intercepto”, que es el valor medio de \(Y\) cuando todos los \(X_i = 0\). \(\beta_{1}\) el coeficiente de \(X_1i\) es un parámetro que representa la pendiente de la recta de regresión para la variable independiente \(X_1i\), y que puede interpretarse como la tasa de cambio de \(Y_{i}\) resultante de un incremento unitario en \(X_{1i}\) cuando el resto de variables independientes se mantienen constantes (en este caso tenemos sólo la variable independiente \(X_1i\));
\(X_{1i}\) es el valor de la variable independiente (explicativa o predictora) del individuo \(i\); y
\(\varepsilon_{ij}\): es el error aleatorio (residual) del \(i\)-ésimo individuo.

Ambos modelos (ANOVA y RL) pueden extenderse para estudiar el efecto de dos (o más) factores o predictores, respectivamente. En el caso de un ANOVA de 2-vías, o bifactorial, como:

\[y_{i,j,k} = \mu + \alpha_{i} + \beta_{j} + \gamma_{i,j} + \varepsilon_{i,j,k}, \ \text{con}\ \ i=1,2,…,l \ \ \text{y}\ \ j=1,2,…,m\]

donde:
\(y_{i,j,k}\) es el valor de la variable respuesta de la k-ésima observación en el tratamiento \(i,j\). Notar que ahora se denomina “tratamiento” a la combinación cruzada de los niveles de los factores involucrados en el experimento aplicada a una unidad experimental u observación;
\(\mu\) es el parrámetro que representa la media general de todos los tratamientos. \(\alpha_{i}\) es el parámetro que representa el efecto fijo del nivel \(i\) del factor \(\alpha\) (con \(i=1, 2, …, l\);
\(\beta_{j}\) es el parámetro que representa el efecto fijo del nivel \(j\) del factor \(\beta\) (con \(j=1, 2, …, m\);
\(\gamma_{i,j}\) es el parámetro que representa el efecto que produce la combinación de los niveles \(i\) del factor \(\alpha\) con los niveles \(j\) del factor \(\beta\). Es un efecto adicional al que produce cada uno de esos niveles separadamente; y
\(\varepsilon_{i,j,k}\): es una variable aleatoria que refleja los efectos no controlados que actúan sobre las unidades experimentales. Residual o error aleatorio de la \(k\)-ésima observación en el \(i,j\)-ésimo tratamiento.

Análogamente, el modelo de RL con dos predictores es:

\[Y_{i} = \beta_{0} + \beta_{1} * X_{1i} + \beta_{2} * X_{2i} + \beta_{3} * X_{1i}X_{2i} + \varepsilon_{i}, \ \text{con}\ \ i=1,2,…,k\]

donde lo novedoso son:
\(\beta_{2}\) que es el parámetro representado por el coeficiente de la variable independiente \(X{2}\); y
\(\beta_{3}\) que es el parámetro que representa el coeficiente de regresión de la interacción entre las variables \(X_1\) y \(X_2\). Es una interacción doble, porque es la interacción entre dos variables independientes.

ANOVA y RL en R

Veremos cómo implementar ambos tipos de modelos (ANOVA y RL) en R mediante ejemplos prácticos.

Ejemplo de ANOVA: lapitas del intermareal rocoso

El Profesor Gerry examinó los efectos de la estación del año (dos niveles, primavera y verano) y la densidad de lapitas adultas (cuatro niveles, 8, 15, 30 y 45 animales - lapitas - por 225 \(cm^{2}\)) sobre la producción individual de masas de huevos (puestas) de la lapa pulmonada (\(Siphonaria\) sp.) del intermareal rocoso. Las lapitas (de unos 10 mm de longitud de caparazón) fueron encerradas en cercados de malla de acero inoxidable de 225 \(cm^{2}\) fijados a la plataforma rocosa. Se aplicaron ocho combinaciones de tratamiento (cuatro densidades en cada una de las dos estaciones) y tres réplicas de cercados por combinación de tratamiento. Las cuatro densidades se utilizaron en ambas estaciones, por lo que se trata de un diseño factorial o cruzado. Una de las cuestiones importantes que se plantearon con este experimento fue si el efecto de la densidad sobre el número de masas de huevos por lapa depende de la estación. Gerry sospechó que el efecto de la densidad sería mayor en verano (cuando su alimento, algas, es escaso), que en primavera (cuando las algas abundan).

Presentación del modelo de ANOVA

Para definir el modelo que usaremos para evaluar las sospechas de Gerry, antes que nada, identificaremos (y definiremos) a la población (biológica y estadística), el tipo la variable respuesta (o dependiente) y las variables independientes, la unidad experimental y el tamaño de la muestra. Además enunciaremos las hipótesis que tendremos que poner a prueba.

Población: todas las cuadrículas (celdas) de 225 \(cm^{2}\) en el intermareal rocoso de la zona de estudio. Notar que la población estadística es el número de puestas en todas las cuadrículas del intermareal rocoso en el área de estudio. Tendré tantas poblaciones diferentes como número de tratamientos aplique, o sea, 2 estaciones * 4 densidades = 8 tratamientos o poblaciones diferentes cuyos parámetros querré estimar (y comparar) a partir de las muestras.

Unidad experimental: es la entidad más pequeña a la que se aplica un tratamiento, i.e. un corralito (o sea, el área encerrada por el corralito). No confundir con las unidades de observación, que son las unidades de medición (una puesta). El tamaño muestral es el número de réplicas de la unidad experimental por tratamiento (n = 3).

Variable dependiente (respuesta): suelen diferenciarse cuatro diferentes tipos de variables por su escala de medición: dos tipos de variables cualitativas - o categóricas - (“nominales” u “ordinales”), y dos cuantitativas (“razón” o de “intervalo”). Las variables cuantitativas suelen también clasificarse en “discretas” o “contínuas”, según si provienen de un conteo o de una medición, respectivamente. En nuestro caso, la variable dependiente es cuantitativa discreta en una escala de razón (donde el cero tiene un sentido de ausencia absoluta), \(Y\): “Producción individual, en número de puestas por corral de 225 \(cm^{2}\), de la lapa pulmonada (\(Siphonaria\) sp.) del intermareal rocoso”.

Variables independientes: dos variables cualitativas o categóricas (factores). \(X_1\): “Estación del año”, con 2 niveles (“Primavera”, “Verano”); \(X_2\): “Densidad de lapitas adultas”, con 4 niveles (“8”, “15”, “30” y “45” lapas adultas por 225 \(cm^{2}\)). La “densidad de lapitas adultas” es, claramente, una variable cuualitativa ordinal (tiene sentido ordenar los niveles de menor a mayor densidad). Para la variable “estación del año” ese orden no es tan claro, pero podría ser razonable considerarla como ordinal argumentando que ambos niveles siguen una secuencia temporal de eventos predecibles (la primavera ocurre antes que el verano, con consecuencias predecibles para los organismos).

Por simplicidad inicialmente consideraremos ambos factores “fijos”. Si consideramos a los efectos como fijos significa que estamos interesados en comparar específicamente los niveles seleccionados. En bioestadística usualmente se enseña que un factor es aleatorio si los niveles usados en el experimento se eligen al azar de un conjunto grande de niveles posibles. Por tanto, si se repitiera el experimento seguramente escogeríamos un grupo diferente de niveles. La elección de un modelo mixto o un modelo fijo es una opción de modelización que debe tener en cuenta el diseño experimental (ver entrada de Ben Bolker para más detalles). Resumiendo, los efectos fijos son constantes en todos los individuos, y los efectos aleatorios varían. Los efectos son fijos si son interesantes en sí mismos, o aleatorios si hay interés en la población subyacente. Por otra parte, si no existen suficientes niveles del factor en los que basar una estimación de la varianza de la población de efectos significa que probablemente debería tratarse a la variable como efectos fijos. En el mejor de los casos, tratar factores con un número reducido de niveles como aleatorios conducirá a estimaciones muy pequeñas y/o imprecisas de los efectos aleatorios y, en el peor de los casos, conducirá a diversas dificultades numéricas como la falta de convergencia, estimaciones de varianza cero, etc. Basándonos en esta síntesis de la revisión de Bolker sería recomendable tratar a “densidad” como un efecto aleatorio y a “estación” como fijo, pero consideraremos a ambas como fijas por simplicidad.

En un experimento factorial, como el del ejemplo, es recomendable comenzar analizando primero la interacción entre los niveles de ambos factores (Estación \(*\) Densidad). Si graficamos “box-plots” de los datos vemos que el efecto de la estación sobre cada nivel del factor Densidad es aproximadamente constante, sugieriendo que no existe una interacción fuerte (figura 1). El número de puestas en primavera es mayor al del verano con aproxiadamente la misma diferencia para los cuatro niveles de densidad analizados (figura 1). El segundo y el tercer box-plot ilustran las diferencias entre los niveles de los factores por separado (“Efectos Principales”) (figura 1). En el segundo box-plot se observa que el número de puestas de primavera es mayor que en verano. El tercer box-plot muestra una reducción aproximadamente constante en el número de puestas con el incremento de la densidad de adultos en los corralitos (figura 1).

rm(list=ls())         # Limpia memoria.
# Lectura de archivos .csv.
lapitas <-read.csv("./0Data/EcolMarina/lapa01.csv",header=T) # sep = ",",
head(lapitas) # muestra las 6 primeras filas (alternativa rápida y
# sintética de View).
str(lapitas) # debo asegurarme que Densidad y Estación sean factores.  
#> 'data.frame':    24 obs. of  3 variables:
#>  $ Densidad: int  8 8 8 8 8 8 15 15 15 15 ...
#>  $ Estacion: chr  "P" "P" "P" "V" ...
#>  $ Huevos  : num  2.88 2.62 1.75 2.12 1.5 ...
# No lo son, debo convertir esas columnas a factores. 
# Doy formato de factor a variables independientes del modelo.
lapitas$Densidad = as.factor(lapitas$Densidad)
lapitas$Estacion = as.factor(lapitas$Estacion)
# Confirmo la correcta estructura de la matriz (data.frame).
str(lapitas)
#> 'data.frame':    24 obs. of  3 variables:
#>  $ Densidad: Factor w/ 4 levels "8","15","30",..: 1 1 1 1 1 1 2 2 2 2 ...
#>  $ Estacion: Factor w/ 2 levels "P","V": 1 1 1 2 2 2 1 1 1 2 ...
#>  $ Huevos  : num  2.88 2.62 1.75 2.12 1.5 ...
# View(lapitas)
# Visualizo los datos.

par(mfrow=c(2,2)) # Para visualizar 2 gráficos en una fila (2 columnas).

# Visualizo diferencias entre los tratamientos.
# Primero ambos tratamientos juntos (Interaccin):
boxplot(Huevos~Estacion*Densidad, data=lapitas,
        col=c("red","blue"),
        main="Puestas de lapitas", xlab="Estación:Densidad",
        ylab="Número de Masas de Huevos") # frame.plot=TRUE,axes=FALSE

# Ahora por separado (Efectos Principales).
# Estación:
boxplot(Huevos~Estacion, data=lapitas,
        col=c("red","blue"),
        main="Puestas de lapitas", xlab="Estación",
        ylab="Número de Masas de Huevos") # frame.plot=TRUE,axes=FALSE

# Densidad:
boxplot(Huevos~Densidad, data=lapitas,
        col=c("red","blue", "green", "magenta"),
        main="Puestas de lapitas", xlab="Densidad",
        ylab="Número de Masas de Huevos") # frame.plot=TRUE,axes=FALSE
\label{fig:figs} Figura 1. Producción individual de masas de huevos (puestas) de la lapa pulmonada ($Siphonaria$ sp.) del intermareal rocoso en dos estaciones del año (niveles: primavera y verano) y cuatro densidades de lapitas adultas (niveles: 8, 15, 30 y 45 animales - lapitas - por $225 cm^2$)

Figura 1. Producción individual de masas de huevos (puestas) de la lapa pulmonada (\(Siphonaria\) sp.) del intermareal rocoso en dos estaciones del año (niveles: primavera y verano) y cuatro densidades de lapitas adultas (niveles: 8, 15, 30 y 45 animales - lapitas - por \(225 cm^2\))

Hipótesis estadísticas: Como mencioné antes, es conveniente comenzar con el análisis de la interacción. Una interacción significativa en el experimento nos diría que los niveles del factor Estación tienen efectos diferentes según el nivel del factor Densidad que se aplique. De hecho, Gerry espera que el número de puestas sea afectado más con la estación cuando se aumenta la densidad de adultos. Podríamos representar al efecto que el factor Estación (P y V) está teniendo cuando estamos aplicando la Densidad \(D_1\) (=8) como: \(P*D_1-V*D_1\). Análogamente para \(D_2\) (=15): \(P*D_2-V*D_2\). Y así sucesivamente hasta \(D_4\) (=45): \(P*D_4-V*D_4\). Si no hay interacción estos efectos deberían ser equivalentes entre sí, y sus diferencias iguales a cero. Puedo enunciar la hipótesis nula de la interacción simplemente como:

\(H0_{E*D}\): No hay interacción entre los factores Estación y Densidad.

Similarmente, para los efectos principales:

\(H0_{E}\): No hay diferencia en las medias del factor Estación.
\(H0_{D}\): No hay diferencia en las medias del factor Densidad.

Ahora especifico el modelo de ANOVA para contrastar esas hipótesis nulas. Formalmente:

\[y_{ijk} = \mu + \alpha_i + \beta_j + \gamma_{ij} + \varepsilon_{ijk}\]

donde:
\(y_{i,j,k}\) es el número de puestas de lapas en el \(k-ésimo\) corralito bajo el tratamiento \(i, j\) de las combinaciones de los niveles de los factores Estación y Densidad, respectivamente.
\(\mu\) es un parámetro que representa la media general de todos los tratamientos.
\(\alpha_i\) parámetro que representa el efecto del nivel \(i\) del Factor \(I\) (Estación), siendo \(\displaystyle\sum_{i}\alpha_{i} = 0\).
\(\beta_j\) parámetro que representael efecto del nivel \(j\) del Factor \(J\) (Densidad), siendo \(\displaystyle\sum_{j}\beta_{j} = 0\).
\(\gamma_{ij}\) parámetro que representa el efecto que produce la combinación del nivel \(i\) de Estación y el nivel \(j\) del Factor Densidad, adicional al que produce cada uno de estos niveles separadamente \(\displaystyle\sum_{i}\gamma_{ij=0}, \displaystyle\sum_{j}\gamma_{ij=0}\).
\(\varepsilon_{ijk}\): efectos no controlados que actúan sobre las unidades experimentales, o sea, los corralitos (es el error aleatorio).

Entonces, las hipótesis nulas se pueden formular más formalmente como:

\(H0_{I*J}: \gamma_{ij} = 0\), \(\forall\ i,j\)
\(H0_{I}: \alpha_{i} = 0\), \(\forall\ i\)
\(H0_{J}: \beta_{j} = 0\), \(\forall\ j\)

Más concretamente, ¿qué son los “efectos” de los niveles de los factores principales (\(\alpha \ \text{y} \ \beta\)) y de la interacción (\(\gamma\))?. Son la diferencias entre la media del tratamiento y la media global:

\(\alpha_{i} = \mu_{i\bullet} - \mu\), y sus estimadores \(\hat{\alpha}_{i} = \bar{y}_{i\bullet} - \bar{y}\). Análogamente:
\(\beta_{j} = \mu_{\bullet j} - \mu\), y sus estimadores \(\hat{\beta}_{j} = \bar{y}_{\bullet j} - \bar{y}\).
\(\gamma_{ij} = \mu_{ij} - \mu_{i\bullet} - \mu_{\bullet j + \mu}\), y sus estimadores \(\hat{\gamma}_{ij} = \bar{y}_{ij} - \bar{y}_{i\bullet} - \bar{y}_{\bullet j} + \bar{y}\).

Donde,
\(\mu_{ij}\) es la media del tratamiento con una combinación de niveles \(ij\) de los factores \(I\) (Estación) y \(J\) (Densidad).
\(\mu_{i\bullet}\) es la media de todas las observaciones del nivel \(i\) del factor Estación. El punto significa que promedia sobre todos los niveles de \(j\) (Densidad) para un dado nivel de \(i\) (Estación).
\(\mu_{\bullet j}\) es la media de todas las observaciones del nivel \(j\) del factor Densidad. El punto significa que promedia sobre todos los niveles de \(i\) (Estación) para un dado nivel de \(j\) (Densidad).
\(\mu\) es la media general de todos los tratamientos.

Equivalentemente, las hipótesis nulas basadas en las diferencias entre medias poblacionales pueden formularse como:

\(H0_{Estación}: \mu_{Est.1} = \mu_{Est.2} = \mu_{Estación}\),
\(H0_{Densidad}: \mu_{Den.1} = \mu_{Den.2} = \mu_{Den.3} = \mu_{Den.4} = \mu_{Densidad}\),

La tabla general para un ANOVA bifactorial (con efectos fijos) es:

F. de V. \(g.l.\) S.C. C.M. \(\mathcal{F}\)
Densidad \(I-1\) \(\displaystyle\sum_{i=1}^{I}nJ(\bar{y}_{i \bullet}-\bar{y})^2\) \(\displaystyle \frac{SC_{Dens.}}{I-1}\) \(\displaystyle \frac{CM_{Dens.}}{CM_{Dentro}}\)
Estación \(J-1\) \(\displaystyle\sum_{J=1}^{J}nI(\bar{y}_{\bullet j}-\bar{y})^2\) \(\displaystyle \frac{SC_{Estac.}}{J-1}\) \(\displaystyle \frac{CM_{Estac.}}{CM_{Dentro}}\)
D * E \((I-1)\cdot(J-1)\) \(\displaystyle\sum_{i=1}^{I}\displaystyle\sum_{J=1}^{J}n(\bar{y}_{ij}-\bar{y}_{i \bullet}-\bar{y}_{\bullet j}-\bar{y})^2\) \(\displaystyle \frac{SC_{Dens.\text{x}Estac.}}{(I-1)\cdot(J-1)}\) \(\displaystyle \frac{CM_{Dens.\text{x}Estac.}}{CM_{Dentro}}\)
Dentro \(I \cdot J(n-1)\) \(\displaystyle\sum_{i=1}^{I}\displaystyle\sum_{J=1}^{J}\displaystyle\sum_{k=1}^{n}(y_{ijk}-\bar{y}_{ij})^2\) \(\displaystyle \frac{SC_{Dentro}}{I \cdot J(n-1)}\)
———— ——————- —————- ————– ————-
TOTAL \(N-1\) \(\displaystyle\sum_{i=1}^{I}\displaystyle\sum_{J=1}^{J}\displaystyle\sum_{k=1}^{n}(y_{ijk}-\bar{y})^2\)

Habiendo presentado el modelo de ANOVA apropiado y formulado las hipótesis estadísticas, ahora pasamos a contrastar las hipótesis con la ayuda de R.

Implementación del ANOVA en R

Para ajustar un ANOVA a los datos del ejercicio usaremos la clásica función aov() del paquete “stats”. La función aov() es realmente una “envoltura” de la función lm() del mismo paquete, usada para ajustar modelos lineales de diseños experimentales balanceados o no. La principal diferencia entre aov() y lm() radica en la forma en que print, summary - y demás funciones usadas para resumir los resultados de los análisis - describen el ajuste. Todas estas funciones aplicadas sobre un objeto “aov” expresan los resultados en el lenguaje tradicional del análisis de varianza (en lugar de en el de los modelos de regresión lineales). Esta diferencia se entenderán mejor cuando comparemos el summary() de la aplicación de aov() y lm() sobre los mismos datos y esecificando el mismo modelo.

Cuadro 1: Antes de presentar la implementación del código en R, es útil recordar aspectos de la notación seguida por R para expresar los modelos fijos o mixtos, cruzados (= factoriales) o anidados (= jerárquicos). Esto es, la forma de transcribir las expresiones matemáticas o “fórmulas” en código de R. La barra inclinada “/” denota anidamiento tal que \(A/B\) significa que \(B\) está anidado dentro de \(A\). Cuando se utiliza “/”, R lo expande automáticamente al efecto principal de la variable principal más la interacción entre la principal y la anidada. Por ejemplo, \(A/B\) se expande automáticamente a \(A + A:B\). Notar que \(A:B\) corresponde al término de la interacción entre los niveles del factor \(A\) y los del factor \(B\). Esto es similar a cómo \(A*B\) se expande automáticamente a \(A + B + A:B\) pero, con la anidación, la variable en el nido nunca aparece fuera de su nido (es decir, no puede haber un efecto principal de \(B\) por sí mismo). El diseño del ejemplo de las lapitas considera a “Estación” y “Densidad” como factores con sus niveles cruzados (diseño factorial o cruzado). Además, en nuestro ejemplo consideramos ambas variables (factores) con efectos fijos. Suponiendo que hubiésemos considerado a “Densidad” como aleatoria, entonces nuestro modelo sería un ANOVA bifactorial mixto, y “Densidad”, así como la interacción entre “Estación:Densidad” serían términos aleatorios del modelo mixto. Para resolverlo la función lmer del paquete lme4 es la más conveniente. La notación \((1 | \text{Densidad})\) denotada por el término entre paréntesis que contiene una barra condicionante se trata de una distribución condicional de posibles interceptos de los niveles de Estación para cada nivel de “Densidad” (la “granularidad” del efecto aleatorio se especifica después de la barra vertical “|”). Todas las observaciones que compartan el mismo nivel de Densidad (“8”, “15”, “30” y “45” lapas adultas por 225 \(cm^{2}\)) tendrán el mismo efecto aleatorio \(\text{Densidad}_i\), y en lmer escribiríamos: \(Y \sim \text{Estación}+(1|\text{Densidad})+(1|\text{Estación:Densidad})\). El “1” significa que querremos tener un intercepto aleatorio por nivel de Densidad y para la interacción entre el factor fijo y el aleatorio. En particular, el término \((1|\text{Estación:Densidad})\) significa que existe un intercepto aleatorio para cada combinación entre Estación y Densidad. Si quisiéramos especificar el intercepto y la pendiente aleatorios para cada nvel de Densidad y la interacción entre los niveles de ambos factores: \(Y \sim \text{Estación}+(\text{Estación}|\text{Densidad})+(\text{Estación}|\text{Estación:Densidad})\). Análogamente, el término \((\text{Estación}|\text{Estación:Densidad})\) significa que existe una pendiente (y un intercepto) aleatorios para el efecto de la Estación para cada combinación de Estación por Densidad. Esto quedará más claro cuando analicemos un modelo de ANOVA mixto. Ahora sí paso a presentar los resultados.

Pasamos a especificar las variables independientes y la variable dependiente con una fórmula con el formato: \(y \sim x_1 + x_2\), donde \(y\) es la variable dependiente, y \(x_1, x_2\) son las variables independientes, o factores, del ANOVA.

A partir de la salida del modelo corrido en R debajo (notar la notación \(\text{Densidad} * \text{Estacion}\), explicada en el Cuadro 1) confirmamos lo que había anticipado con los box-plots, que ambos factores tienen efectos significativos sobre la producción de huevos, pero la interacción entre los niveles de ambos factores no lo es.

En R:

# ANOVA de 2-factores cruzados:
aov.lapitas <- aov(Huevos ~ Densidad * Estacion, data = lapitas)
summary(aov.lapitas)
#>                   Df Sum Sq Mean Sq F value   Pr(>F)    
#> Densidad           3  5.284   1.761   9.669 0.000704 ***
#> Estacion           1  3.250   3.250  17.842 0.000645 ***
#> Densidad:Estacion  3  0.165   0.055   0.301 0.823955    
#> Residuals         16  2.915   0.182                     
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La interacción entre los niveles de ambos factores es despreciable, no hay evidencia de interacción entre los niveles de Estación y los de Densidad (\(p = 0.823955\)). Esto se representa en la figura 2 por las curvas de interacción paralelas al visualizar cómo el número de puestas se diferencia entre los niveles de la Estación (Primavera y Verano) para cada nivel del factor Densidad. En otras palabras, la diferencia entre las medias de los niveles de Estación se mantienen casi iguales para todos los niveles del factor Densidad, en oposición a lo que esperaba Gerry (figura 2). Los efectos de la densidad de lapas sobre el número de puestas no depende de la estación del año considerada.


# Visualizo Interacción no significativa:
interaction.plot(x.factor = lapitas$Densidad,
                 trace.factor = lapitas$Estacion,
                 response = lapitas$Huevos,
                 type = "b", # "l"
                 pch = c(19,21,23),
                 main ="Interacción entre factores",
                 xlab = expression("Densidad de lapitas"),
                 ylab = expression("Número de masas de huevos"),
                 trace.label = "Estación",  # label for legend
                 xpd = FALSE) #,  # 'clip' legend at border)
\label{fig:figs} Figura 2. Interacción entre los factores Estación del año y Densidad de Adultos en la producción individual de puestas de la lapa pulmonada del intermareal rocoso

Figura 2. Interacción entre los factores Estación del año y Densidad de Adultos en la producción individual de puestas de la lapa pulmonada del intermareal rocoso

En contraste, por un lado puede afirmarse que existen diferencias entre las medias de los niveles de Densidad (\(p = 0.000704\)) y, por otra parte, también existen diferencias entre las medias de los niveles de Estación (\(p = 0.000645\)). Los efectos principales, tanto Densidad como Estación, son altamente significativos (\(p < 0.001\)) y rechazamos las hipótesis nulas (\(H_0\)) planteadas para las comparaciones entre las medias de los niveles de ambos factores por separado. La figura 1 da indicios de la magnitud de las diferencias entre las medias de los niveles de Densidad y Estación. A partir de la prueba de hipótesis con el ANOVA de 2 vías ejecutado en R podemos concluir que el número medio de puestas de lapitas en Primavera es significativamente mayor al número medio de Verano. Podemos concluir esto porque el factor tiene sólo dos niveles. En cambio, la prueba de hipótesis para el factor Densidad no nos permite identificar el origen de las diferencias entre las medias de sus niveles porque ese factor tiene más de dos niveles (4 niveles).

CONCLUSIÓN: Hay evidencias sólidas de efectos de la densidad y de la estación del año en la producción de masas de huevos, pero el efecto de la densidad es independiente del efecto de la estación del año. Los efectos combinados de los niveles de ambos factores se suman (o se restan) y por eso decimos que sus efectos son “aditivos”, no interactúan entre sí.

Ejemplo de Regresión lineal: relación de la abundancia de aves con la fragmentación boscosa

Para investigar los efectos de la fragmentación de bosques, una investigadora relacionó la abundancia de aves (número de aves por parcela) de un manglar (bosque de mangles) con una serie de variables del paisaje descriptoras de los cambios en el bosque o que presumiblemente afectan a la abundancia de aves. Estas variables son:

\(X_1\): el área de la parcela boscosa (en hectáreas) - esto es, el área del parche de bosque inmerso en un área sin árboles;
\(X_2\): el tiempo de aislamiento de la parcela boscosa (en años) - esto es, el tiempo transcurrido desde que se desforestó el área que rodea a la parcela;
\(X_3\): la distancia a la parcela boscosa más cercana (en kilometros) - esto es, distancia a otra parcela boscosa que se salvó de la desforestación;
\(X_4\): distancia a la parcela boscosa más grande (en km);
\(X_5\): la intensidad del pastoreo en cada parcela boscosa (toma valores de 1 a 5). Es el pastoreo por ganado criado en las inmediaciones; y
\(X_6\): la altitud (en metros) - esto es, altura de la parcela respecto del nivel del mar.

Las variables fueron medidas en un total de 56 parcelas de manglar en Belice.

Análisis prliminar en R

Primero exploramos la correlación entre pares de variables con la función “pairs.panels” del paquete “psych” (figura 3):

# Preparación de datos e Introducción del Problema:
rm(list=ls())
graphics.off()
aves <- read.table("./0Data/EcolMarina/aves01.csv", header = T, sep = ",")
str(aves)
#> 'data.frame':    56 obs. of  7 variables:
#>  $ abund   : num  5.3 2 1.5 17.1 13.8 14.1 3.8 2.2 3.3 3 ...
#>  $ area    : num  0.1 0.5 0.5 1 1 1 1 1 1 1 ...
#>  $ tiempo  : int  17 65 85 19 67 20 30 65 20 85 ...
#>  $ dist.1  : int  39 234 104 66 246 234 467 284 156 311 ...
#>  $ dist.2  : int  39 234 311 66 246 285 467 1829 156 571 ...
#>  $ pastoreo: int  2 5 5 3 5 3 5 5 4 5 ...
#>  $ altura  : int  160 60 140 160 140 130 90 60 130 130 ...
# View(aves)

library(psych)
pairs.panels(aves[,c("abund", "area", "tiempo", "dist.1", "dist.2", "pastoreo", "altura")], 
             pch = 21,           # pch symbol
             cex = 1.5,          # tamaño de la fuente
             jiggle = FALSE,     # If TRUE, data points are jittered
             method = "pearson", # correlation method
             cor = TRUE,         # If TRUE, reports correlations
             stars = TRUE,       # If TRUE, adds significance level with stars
             hist.col = 4,       # Histograms color
             lm = FALSE,         # If TRUE, plots linear fit rather than the LOESS
             # (smoothed) fit
             ci = FALSE,          # If TRUE, adds confidence intervals
             density = TRUE,     # show density plots
             ellipses = TRUE,    # show correlation ellipses
             smooth = TRUE,      # If TRUE, draws loess smooths
             scale = FALSE,      # If TRUE, scales the correlation text font
             )
\label{fig:figs} Figura 3. Correlograma para evaluar la relación entre agentes de disturbio en un manglar de Belice afectado por la fragmentación boscosa

Figura 3. Correlograma para evaluar la relación entre agentes de disturbio en un manglar de Belice afectado por la fragmentación boscosa

Del análisis preliminar de la correlación entre las variables (incluída la variable dependiente - o respuesta - “Abundancia de aves”) vemos que, particularmente la variable dependiente, tiene una fuerte correlación con dos variables independientes - o predictoras - que son: “Tiempo” y “Pastoreo” (figura 3). Comprobamos además que la correlación entre ellas es alta y significativa, por lo que probablemente alguna de ellas podría ser excluída de un modelo de Regresión lineal en una selección de variables por pasos por ser potencialmente colineales.

La función “rcorr()” [paquete Hmisc] puede ser usada para computar el nivel de significancia del los coeficientes de correlación de Pearson y Spearman, confirmando los resultados y conclusiones preliminares:

library(Hmisc)
#> 
#> Attaching package: 'Hmisc'
#> The following object is masked from 'package:psych':
#> 
#>     describe
#> The following objects are masked from 'package:base':
#> 
#>     format.pval, units
cor.aves <- rcorr(as.matrix(aves), type = "pearson") # Pearson's seems to be the default
# Extract the correlation coefficients
cor.aves$r
#>                abund         area       tiempo     dist.1      dist.2
#> abund     1.00000000  0.255970206 -0.503357741  0.2361125  0.08715258
#> area      0.25597021  1.000000000  0.001494192  0.1083429  0.03458035
#> tiempo   -0.50335774  0.001494192  1.000000000 -0.1132175  0.08331686
#> dist.1    0.23611248  0.108342870 -0.113217524  1.0000000  0.31717234
#> dist.2    0.08715258  0.034580346  0.083316857  0.3171723  1.00000000
#> pastoreo -0.68251138 -0.310402417  0.635567104 -0.2558418 -0.02800944
#> altura    0.38583617  0.387753885 -0.232715406 -0.1101125 -0.30602220
#>             pastoreo     altura
#> abund    -0.68251138  0.3858362
#> area     -0.31040242  0.3877539
#> tiempo    0.63556710 -0.2327154
#> dist.1   -0.25584182 -0.1101125
#> dist.2   -0.02800944 -0.3060222
#> pastoreo  1.00000000 -0.4071671
#> altura   -0.40716705  1.0000000
# Extract p-values
cor.aves$P
#>                 abund        area       tiempo     dist.1     dist.2
#> abund              NA 0.056889437 7.677519e-05 0.07978817 0.52301973
#> area     5.688944e-02          NA 9.912798e-01 0.42671752 0.80025683
#> tiempo   7.677519e-05 0.991279840           NA 0.40608130 0.54153936
#> dist.1   7.978817e-02 0.426717521 4.060813e-01         NA 0.01722718
#> dist.2   5.230197e-01 0.800256832 5.415394e-01 0.01722718         NA
#> pastoreo 6.896940e-09 0.019897121 1.420679e-07 0.05701822 0.83763768
#> altura   3.315592e-03 0.003150012 8.434937e-02 0.41915824 0.02180537
#>              pastoreo      altura
#> abund    6.896940e-09 0.003315592
#> area     1.989712e-02 0.003150012
#> tiempo   1.420679e-07 0.084349374
#> dist.1   5.701822e-02 0.419158243
#> dist.2   8.376377e-01 0.021805368
#> pastoreo           NA 0.001843259
#> altura   1.843259e-03          NA

Las dos variables regresoras más fuertemente correlacionadas con la abundancia de aves son tiempo \(= -68\), pastoreo \(= -50\). Arbitrariamente y con fines puramente exploratorios ajustamos un modelo de RLM (Regresión Lineal Múltiple) que relacione la abndancia de aves del bosque (\(Y\)) con las dos variables regresoras con mayor correlación lineal con \(Y\), y lo llamamos “mod1”:

\[\hat{y} = b_0 + b_1 * x_1 + b_2 * x_2\]

Ajuste de RLM en R e interpretación de los resultados

mod1 <- lm(abund ~ tiempo + pastoreo, data=aves)
anova(mod1)
summary(mod1)
#> 
#> Call:
#> lm(formula = abund ~ tiempo + pastoreo, data = aves)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -19.4688  -4.7551   0.3169   4.4733  19.3304 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 34.48115    2.41674  14.268  < 2e-16 ***
#> tiempo      -0.04898    0.05415  -0.905     0.37    
#> pastoreo    -4.43984    0.94182  -4.714  1.8e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 7.931 on 53 degrees of freedom
#> Multiple R-squared:  0.4739, Adjusted R-squared:  0.4541 
#> F-statistic: 23.87 on 2 and 53 DF,  p-value: 4.049e-08

El ajuste de lm() y aov() es idéntico, pero el informe derivado de anova() y summary() de un modelo ajustado con lm() es diferente. Es importante aclarar que esa diferencia se debe a que las pruebas \(\mathcal{F}_{global}\) reportadas por anova() son secuenciales. La \(\mathcal{F}_{global}\) para “tiempo” evalúa la importancia del “tiempo” en presencia de únicamente el intercepto. En otras palabras, se cuantifica la explicación de la variabilidad de la abundancia de aves para la primera variable predictora analizada (i.e. Tiempo) sin considerar la proporción de la variabilidad explicada por la segunada variable en la tabla de ANOVA (i.e. Pastoreo). La \(\mathcal{F}_{global}\) del pastoreo evalúa su efecto en presencia de el intercepto y el “tiempo”. Vemos así que, en función de lo presentado en la tabla de ANOVA (salida de anova(mod1)), ambas variables independientes son significativas. En contraste, la función summary reporta los resultados no secuencialmente. Las pruebas \(t\) del summary son el impacto marginal de las variables en cuestión, dada la presencia de todas las demás variables. Entonces, es importante valorar la importancia de las variables independientes en el modelo de RLM a partir de la tabla del “Summary”, summary(mod1) (no de la producida por anova(mod1)).

Si comparásemos los resultados de las pruebas de hipótesis para los coeficientes de las variables independientes del summary() con los cálculos a mano de los coeficientes parciales del modelo confirmaríamos que son casi iguales (tanto los valores para los estadísticos \(t\), como los valores \(p\)). El análisis de los coeficientes de regresión obtenidos con summary() sugiere que podríamos llegar a conclusiones erradas si incluímos ambas variables simultáneamente. El valor de los coeficientes en summary() indica cuánto cambia en promedio la variable dependiente (Abundancia de Aves) ante un cambio de una unidad en una variable independiente (por ejemplo, Pastoreo), manteniendo constantes las demás variables del modelo (en el ejemplo, Tiempo). Entonces, esperamos en promedio una disminución de 4.4 aves por un aumento de una unidad en la intensidad del pastoreo, cuando Tiempo se mantiene constante. Análogamente, esperamos una disminución promedio de 0.05 aves por un incremento de un año en el aislamiento de la parcela (Tiempo), cuando la intensidad del Pastoreo se mantiene constante. Evidentemente el efecto del Tiempo en la abundancia de aves es pequeño (y no significativo) cuando la intensidad de Pastoreo no cambia. Entonces, los resultados apoyan la idea de una colinealidad entre ambas variables como consecuencia de su elevada correlación lineal. Será necesaria una selección por pasos (“stepwise”) de las variables predictoras de la Abundancia de Aves mediante criterios de selección como el AIC para definir qué variable eliminar del modelo (y cuál dejar), pero todo indica que Tiempo saldría del mejor modelo y que la intensidad del Pastoreo es la variable predictora de la abundancia de aves más importante. Por otro lado es fundamental un análisis completo de los residuos para contrastar los supuestos del modelo de Regresión ganador (normalidad, homogeneidad de varionzas, independencia).

Cuadro 2: En un ANOVA nos interesa evaluar diferencias entre los niveles de uno o más factores (además de determinar la importancia, o significancia, de los factores del análisis) y siempre debemos considerar a las variables independientes como factores (usando “as.factors()” en R), como en el ejercicio de las lapitas de Gerry. En un modelo de Regresión lineal estamos interesados en ver si existe una relación lineal entre la variable dependiente (respuesta) y alguna/s de las variables independientes del modelo. Si alguna/s de las variable/s independiente/s es/son categórica/s y además es/son nominal/es (o sea, no tiene sentido ordenarlas: como “sexo”, “color” o “ëspecie”) entonces esas variables categóricas nominales deben ser convertidas a factor en un modelo de Regresión (con la función as.factor()). El modelo de regresión estimará un parámetro para cada nivel (o categoría) de cada una de esas variables nominales. Pero si la/s variable/s categórica/s en el modelo de Regresión es/son ordinal/es (o sea, tiene sentido ordenar sus categorías) y pensamos que los “valores” de las categorías representan bien la relación entre ellas, entonces es mejor representar los diferentes “valores” como números enteros que conservan el orden original. Este es el caso de Pastoreo si consideramos que la “intensidad de pastoreo 5” es 5 veces mayor que la “intensidad 1”, la “2” es 2 veces mayor que la “1”, etc. O sea, que los números representan bien las diferencias entre las intensidades, En ese caso es mejor NO convertir Pastoreo como factor (es mejor dejarlo como número entero, como hicimos en el ejemplo de las aves). En esta situación los grados de libertad de Pastoreo dan igual a 1 porque, al considerarlo a la variable Pastoreo como numérica (“escalar”, no como “factor”) los GL se calculan como el número de parámetros (2: intercepto y pendiente) menos 1 (i.e. 2-1). Para el caso de una RLM con una o más variables independientes cuantitativas y uno o más factores estaríamos hablando de un ANCOVA, ver Wiki haciendo click.

Noten que el \(\mathcal{F}_{global}\) calculado a partir de la tabla de la salida de anova() es idéntico al reportado en la salida de summary() (\(\mathcal{F}_{global} = 23.87\)), como es esperable. Basado en el valor \(p\) de la prueba de hipótesis global (ver saliida de summary()), rechazo \(H_0\) con una probabilidad de error \(<<\) 5%. Entonces, hay evidencia significativa para afirmar que existe una relación lineal entre la variable dependiente y al menos una de las variables independientes.

anova(mod1)

SCReg = 1605.8+1398.0  # Calculo la Suma de Cuadrados explicada. 
CMReg = SCReg/2        # Calculo los Cuadrados Medios de la parte explicada.
Fvalue = CMReg / 62.91 # Calculo el F global de la prueba.
Fvalue
#> [1] 23.87379

summary(mod1)
#> 
#> Call:
#> lm(formula = abund ~ tiempo + pastoreo, data = aves)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -19.4688  -4.7551   0.3169   4.4733  19.3304 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 34.48115    2.41674  14.268  < 2e-16 ***
#> tiempo      -0.04898    0.05415  -0.905     0.37    
#> pastoreo    -4.43984    0.94182  -4.714  1.8e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 7.931 on 53 degrees of freedom
#> Multiple R-squared:  0.4739, Adjusted R-squared:  0.4541 
#> F-statistic: 23.87 on 2 and 53 DF,  p-value: 4.049e-08

¿Cuáles son las diferencias entre el \(R^2\) y el \(R^2\)-ajustado?.

  • \(R^2 = 0.4739\) (significa que el modelo explica el 47% de la variabilidad en la abundancia de aves.

  • \(R^2 \ ajustado = 0.4541\) (como explico debajo, es menor al \(R^2\)).

El \(R^2\), o coeficiente de determinación, se utiliza para explicar el grado en que las variables independientes (variables predictoras) explican la variación de la/s variable/s dependiente/s (variables respuesta). Toma valores de 0 a 1. Por ejemplo, un \(R^2\) = 0.9, indica que el 90% de la variación de la variable dependiente se explica por las variables independientes en el modelo. En general, un \(R^2\) más alto indica un mejor ajuste del modelo. El \(R^2\) tiene un problema inherente: las variables independientes adicionales harán que el \(R^2\) se mantenga o aumente (esto se debe a la forma en que se calcula el \(R^2\)). Por lo tanto, aunque las variables independientes adicionales no muestren ninguna relación con la variable dependiente, el \(R^2\) aumentará. El \(R^2\)-ajustado también indica qué tan bueno es el ajuste del modelo para explicar la variabilidad de la variable dependiente, pero es corregido por el número de términos (variables independientes) en un modelo. Si se añaden más y más variables independientes innecesarias a un modelo, el \(R^2\)-ajustado disminuirá. Si se añaden más variables independientes valiosas, el \(R^2\)-ajustado aumentará. El \(R^2\)-ajustado siempre será menor o igual que el \(R^2\). Esencialmente, el \(R^2\)-ajustado examina si las variables independientes adicionales están contribuyendo al modelo.

Intervalo de confianza de los parámetros

Ahora calcularemos los Intervalos de Confianza de 95% para los coeficientes correspondientes para el modelo mod1 usando la función confint(). Escriban en la consola: ?confint(), para saber qué clase de objeto requiere la fucnión, y qué argumentos usar para especificar el nivel de significancia y los parámetros del modelo. Para interpretar los resultados no olviden fijarse si el IC contiene al valor 0 o no (y recuerden para quién es el IC que están construyendo, ¿qué significa que el coeficiente de un modelo lineal sea igual a 0?).

ICcoef.mod1 <- confint(mod1, parm = c("tiempo", "pastoreo"), level = 0.95)
ICcoef.mod1
#>               2.5 %     97.5 %
#> tiempo   -0.1575948  0.0596298
#> pastoreo -6.3288913 -2.5507855

El IC del coeficiente de tiempo indica que su estimación no difiere de 0 (la pendiente del cambio de la variable dependiente por unidad de cambio de la variable Tiempo no es diferenciabl de 0). Entonces, el efecto del tiempo en la abundancia de aves puede despreciarse. Contrariamente, el IC del coeficiente de pastoreo no comprende al cero, indicando que la pendiente del cambio de la variable dependiente por unidad de cambio de la variable Pastoreo es distinta de cero.

Cuadro 3: ¿Cuándo comprobamos los supuestos del modelo? Aquí no me explayaré en la respuesta a esta pregunta, pero pueden encontrar una explicación aquí. Mi sugerencia es hacer la evaluación de los residuos del modelo inicial (el modelo lineal con todas las variables independientes e interacciones), así como del modelo finalmente seleccionado. De hecho, Zuur et al. (2009) sugieren hacer la comprobación iterativamente siguiendo un protocolo que parte del modelo lineal con todos los términos fijos, pasando sucesivamente por modelos lineales mixtos, generalizados, y/o aditivos hasta satisfacer los supuestos. Si nuestras observaciones son independientes y preferimos ajustar un modelo lineal (antes de optar por modelos más complejos) una posibilidad es transformar las variables predictoras y/o la variable respuesta, dependiendo de las causas de los patrones observados en los residuos del modelo ajustado inicialmente. Las pruebas y diagnósticos recomendados, así como las transformaciones potencialmente existosas en diferentes contextos constituyen material de otro tutorial.

Métodos de selección de variables en RLM: Regresión por Pasos

Por simplicidad en este tutorial ejemplifico la selección de variables independientes con el método de Regresión por Pasos (“Stepwise Regression”), aunque existen métodos alternativos que pueden dar mejores resultados. Ahora evaluaremos al modelo que considera a todas las variables independientes con los tres métodos de selección de variables independientes (“Forward”, “Backward”, “Both”) para RLM usando el criterio de Akaike (AIC).Recomiendo la función step(), del paquete stats. Para comenzar, la misma sugerencia de siempre cuando se introduce una nueva función: ?step(). Uno de los argumentos más importantes de la función es el que define el método a utilizar (direction). Notar debajo que obviamos las interacciones entre las variables independientes, aunque formalmente deberían ser consideradas en el modelo completo. En efecto, esta es una limitación importante del método de selección de variables de Regresión por Pasos. Como se menciona acá, inevitablemente violaremos el “Principio de Jerarquía” al implementar la Regresión por Pasos para un modelo que incluye las interecciones, ya que durante el proceso nada impedirá que se seleccione un modelo con un término para una interacción pero con alguno, o todas, los términos principales participantes de la interacción ausentes.

Método Forward: comienza sin predictores en el modelo y añade iterativamente los predictores más contributivos y se detiene cuando la mejora deja de ser estadísticamente significativa. Si eligen “forward” es recomendable construir un modelo nulo para decirle a la función que comience por lo más simple, que es un modelo lineal que sólo contiene al intercepto (es una especie de modelo neutro con el cual puedo comparar mis modelos). Constrúyanlo así:

nulo.lm <- lm(abund ~ 1., data=aves) # modelo con sólo el intercepto.

Otro parámetro es “scope”, útil para especificar el modelo más extenso permitido (que será el que contiene a todas las variables independientes que consideró en su estudio la investigadora). En definitiva lo que hará el algoritmo del método forward es comenzar por el modelo nulo y calcularle el AIC, y en el primer paso de selección escogerá alternativamente una variable independiente por vez y calculará el AIC para cada modelo con el intercepto + 1 variable independiente. Comparará a todos los modelos con el intercepto + sólo una variable y se quedará con el de menor AIC. En el segundo pasa sumará una variable más al modelo escogido en el paso anterior, probando con cada una de las restantes, y evaluará si alguno de los modelos con el intercepto y dos variables independientes es mejor que el modelo seleccionado antes (con una variable más el intercepto). Y así sucesivamente hasta llegar al modelo completo indicado en “scope” (a menos que frene antes porque la adición de las variables restantes tiene un AIC mayor al del modelo seleccionado al momento).

# Método Forward:

nulo.lm <- lm(abund ~ 1., data=aves) # modelo con sólo el intercepto.

step(object = nulo.lm,   # "lm" o "glm". Lo usa como modelo inicial.
     scope = abund ~ area + tiempo + dist.1 + dist.2 +pastoreo +
       altura, # útil para especificar el modelo más extenso permitido.
     direction = "forward" # allow it remove predictors but not add them.
     )
#> Start:  AIC=266.82
#> abund ~ 1
#> 
#>            Df Sum of Sq    RSS    AIC
#> + pastoreo  1   2952.35 3385.6 233.71
#> + tiempo    1   1605.83 4732.1 252.46
#> + altura    1    943.52 5394.4 259.80
#> + area      1    415.27 5922.7 265.03
#> + dist.1    1    353.33 5984.6 265.61
#> <none>                  6337.9 266.82
#> + dist.2    1     48.14 6289.8 268.39
#> 
#> Step:  AIC=233.71
#> abund ~ pastoreo
#> 
#>          Df Sum of Sq    RSS    AIC
#> <none>                3385.6 233.71
#> + altura  1    88.519 3297.1 234.22
#> + tiempo  1    51.473 3334.1 234.85
#> + dist.2  1    29.361 3356.2 235.22
#> + dist.1  1    25.649 3359.9 235.28
#> + area    1    13.651 3371.9 235.48
#> 
#> Call:
#> lm(formula = abund ~ pastoreo, data = aves)
#> 
#> Coefficients:
#> (Intercept)     pastoreo  
#>      34.369       -4.981

Método Backward: el algoritmo del método “Backward” funciona de manera similar pero comienza por el modelo de regresión completo, incluyendo todos los posibles predictores. Luego, en cada “paso” va probando todas las formas posibles de eliminar una de las variables, y se queda con el mejor modelo en términos del valor del AIC (el de menor AIC). Este se convierte en nuestro nuevo modelo de regresión; y a continuación vuelve a comparar todas las posibles eliminaciones de una variable independiente, eligiendo la opción con menor AIC. Se detiene cuando se tiene un modelo en el que todos los predictores son estadísticamente significativos y un valor de AIC menor que cualquiera de los otros modelos posibles que podría producir al eliminar uno de sus predictores. Obviamente, cuando usen este método en la función step deberán especificar en object al modelo de regresión completo (que deberán construir y correr en R antes de ejecutar la función).

# Método Backward:

# Modelo completo (con todas las vars. independientes):
aves.lm <- lm(abund ~ area + tiempo + dist.1 + dist.2 + pastoreo + altura, data=aves)

step(object = aves.lm,   # "lm" o "glm". Lo usa como modelo inicial.
     direction = "backward")
#> Start:  AIC=238.67
#> abund ~ area + tiempo + dist.1 + dist.2 + pastoreo + altura
#> 
#>            Df Sum of Sq    RSS    AIC
#> - area      1      2.29 3096.5 236.71
#> - dist.1    1     31.24 3125.5 237.23
#> - dist.2    1     73.94 3168.2 237.99
#> - tiempo    1     87.54 3181.8 238.23
#> <none>                  3094.2 238.67
#> - altura    1    150.73 3245.0 239.33
#> - pastoreo  1    612.49 3706.7 246.78
#> 
#> Step:  AIC=236.71
#> abund ~ tiempo + dist.1 + dist.2 + pastoreo + altura
#> 
#>            Df Sum of Sq    RSS    AIC
#> - dist.1    1     32.61 3129.1 235.30
#> - dist.2    1     76.94 3173.5 236.09
#> - tiempo    1     86.04 3182.6 236.25
#> <none>                  3096.5 236.71
#> - altura    1    182.64 3279.2 237.92
#> - pastoreo  1    675.82 3772.3 245.77
#> 
#> Step:  AIC=235.3
#> abund ~ tiempo + dist.2 + pastoreo + altura
#> 
#>            Df Sum of Sq    RSS    AIC
#> - tiempo    1     82.05 3211.2 234.75
#> - dist.2    1    111.28 3240.4 235.25
#> <none>                  3129.1 235.30
#> - altura    1    163.85 3293.0 236.16
#> - pastoreo  1    809.22 3938.3 246.18
#> 
#> Step:  AIC=234.75
#> abund ~ dist.2 + pastoreo + altura
#> 
#>            Df Sum of Sq    RSS    AIC
#> - dist.2    1     85.89 3297.1 234.22
#> <none>                  3211.2 234.75
#> - altura    1    145.05 3356.2 235.22
#> - pastoreo  1   1888.71 5099.9 258.65
#> 
#> Step:  AIC=234.22
#> abund ~ pastoreo + altura
#> 
#>            Df Sum of Sq    RSS    AIC
#> - altura    1     88.52 3385.6 233.71
#> <none>                  3297.1 234.22
#> - pastoreo  1   2097.34 5394.4 259.80
#> 
#> Step:  AIC=233.71
#> abund ~ pastoreo
#> 
#>            Df Sum of Sq    RSS    AIC
#> <none>                  3385.6 233.71
#> - pastoreo  1    2952.3 6337.9 266.82
#> 
#> Call:
#> lm(formula = abund ~ pastoreo, data = aves)
#> 
#> Coefficients:
#> (Intercept)     pastoreo  
#>      34.369       -4.981

Método Combinado: El algoritmo del método “both” es una combinación del Forward y el Backward. Comienza con todos los predictores del modelo (modelo completo) y en el primer paso compara su AIC con cada una de las versiones de ese modelo sin una de sus variables. Similarmente, al modelo con menor AIC resultante lo compara en el paso 2 con cada una de las versiones de ese modelo sin una de sus variables pero incluyendo además entre los modelos competidores al modelo completo. Al modelo seleccionado en el paso 2 lo compara en el tercer paso con sus versiones sin una de sus variables, pero además con el modelo seleccionado en el paso 2 con la adición de la variable eliminada en el paso 1 y con el modelo seleccionado en el paso 2 con la adición de la variable eliminada en el paso 2. El algoritmo elimina y añade variables de esa forma en cada paso hasta quedarse con el modelo con menor AIC.

# Método Combinado:
step(object = aves.lm,   # "lm" o "glm". Lo usa como modelo inicial.
     scope = abund ~ area + tiempo + dist.1 + dist.2 + pastoreo
     + altura, # útil para especificar el modelo más extenso permitido.
     direction = "both")
#> Start:  AIC=238.67
#> abund ~ area + tiempo + dist.1 + dist.2 + pastoreo + altura
#> 
#>            Df Sum of Sq    RSS    AIC
#> - area      1      2.29 3096.5 236.71
#> - dist.1    1     31.24 3125.5 237.23
#> - dist.2    1     73.94 3168.2 237.99
#> - tiempo    1     87.54 3181.8 238.23
#> <none>                  3094.2 238.67
#> - altura    1    150.73 3245.0 239.33
#> - pastoreo  1    612.49 3706.7 246.78
#> 
#> Step:  AIC=236.71
#> abund ~ tiempo + dist.1 + dist.2 + pastoreo + altura
#> 
#>            Df Sum of Sq    RSS    AIC
#> - dist.1    1     32.61 3129.1 235.30
#> - dist.2    1     76.94 3173.5 236.09
#> - tiempo    1     86.04 3182.6 236.25
#> <none>                  3096.5 236.71
#> - altura    1    182.64 3279.2 237.92
#> + area      1      2.29 3094.2 238.67
#> - pastoreo  1    675.82 3772.3 245.77
#> 
#> Step:  AIC=235.3
#> abund ~ tiempo + dist.2 + pastoreo + altura
#> 
#>            Df Sum of Sq    RSS    AIC
#> - tiempo    1     82.05 3211.2 234.75
#> - dist.2    1    111.28 3240.4 235.25
#> <none>                  3129.1 235.30
#> - altura    1    163.85 3293.0 236.16
#> + dist.1    1     32.61 3096.5 236.71
#> + area      1      3.65 3125.5 237.23
#> - pastoreo  1    809.22 3938.3 246.18
#> 
#> Step:  AIC=234.75
#> abund ~ dist.2 + pastoreo + altura
#> 
#>            Df Sum of Sq    RSS    AIC
#> - dist.2    1     85.89 3297.1 234.22
#> <none>                  3211.2 234.75
#> - altura    1    145.05 3356.2 235.22
#> + tiempo    1     82.05 3129.1 235.30
#> + dist.1    1     28.61 3182.6 236.25
#> + area      1      0.22 3211.0 236.74
#> - pastoreo  1   1888.71 5099.9 258.65
#> 
#> Step:  AIC=234.22
#> abund ~ pastoreo + altura
#> 
#>            Df Sum of Sq    RSS    AIC
#> - altura    1     88.52 3385.6 233.71
#> <none>                  3297.1 234.22
#> + dist.2    1     85.89 3211.2 234.75
#> + dist.1    1     57.37 3239.7 235.24
#> + tiempo    1     56.66 3240.4 235.25
#> + area      1      0.82 3296.2 236.21
#> - pastoreo  1   2097.34 5394.4 259.80
#> 
#> Step:  AIC=233.71
#> abund ~ pastoreo
#> 
#>            Df Sum of Sq    RSS    AIC
#> <none>                  3385.6 233.71
#> + altura    1     88.52 3297.1 234.22
#> + tiempo    1     51.47 3334.1 234.85
#> + dist.2    1     29.36 3356.2 235.22
#> + dist.1    1     25.65 3359.9 235.28
#> + area      1     13.65 3371.9 235.48
#> - pastoreo  1   2952.35 6337.9 266.82
#> 
#> Call:
#> lm(formula = abund ~ pastoreo, data = aves)
#> 
#> Coefficients:
#> (Intercept)     pastoreo  
#>      34.369       -4.981

CONCLUSIÓN: Los tres métodos seleccionaron al mismo modelo, el modelo con Pastoreo como única variable independiente que contribuye de forma importante a la abundancia de aves.