1 Introducción

Interpretar las interacciones entre variables nn análisis de regresión es un poco más sencillo que en ANOVA, así que veamos un ejemplo con una ecuación lineal donde \(y\) es una función de dos variables \(x_1\) y \(x_2\):

\[\begin{equation} y(x_1, x_2) = a + b_1 x_1 + b_2 x_2 \tag{1.1} \end{equation}\]

En un espacio lineal el producto \(x_1 \cdot x_2\) no puede existir, por lo que la ecuación:

\[\begin{equation} y(x_1, x_2) = a + b_1 x_1 + b_2 x_2 + b_3 x_1 x_2 \tag{1.2} \end{equation}\]

tampoco existe en un espacio lineal. Sin embargo si \(x_3 = x_1 \cdot x_2\), de modo que la pendiente \(b_3\) de la interacción \(x_1 \cdot x_2\) es lineal con respecto de \(x_3\) la ecuación

\[ y(x_1, x_2) = a + b_1 x_1 + b_2 x_2 + b_3 x_3 \] sí existe en un espacio lineal. Es en este contexto que se estiman las interacciones estadísticas entre variables independientes \(x\). Geométricamente, una ecuación lineal (ecuación (1.1)) sin interacciones se ve así:

Gráfica de perspectiva de una ecuación lineal sin interacciones con $a=0$, $b_1 = 1$ y $b_2 = 1$.

Figure 1.1: Gráfica de perspectiva de una ecuación lineal sin interacciones con \(a=0\), \(b_1 = 1\) y \(b_2 = 1\).

y cuando hay interacciones (1.2):

Gráfica de perspectiva de una ecuación lineal con interacciones con $a=0$, $b_1 = 1$, $b_2 = 1$ y $b3 = 0.5$.

Figure 1.2: Gráfica de perspectiva de una ecuación lineal con interacciones con \(a=0\), \(b_1 = 1\), \(b_2 = 1\) y \(b3 = 0.5\).

En las gráficas se puede apreciar cómo los valores de \(y\) aumentan más rápido cuando más grandes son los valores de \(x_1\) y \(x_2\). Este efecto, sin embargo no es universal pues es posible que la pendiente de \(b_3 < 0\):

Gráfica de perspectiva de una ecuación lineal con interacciones con $a=0$, $b_1 = 1$, $b_2 = 1$ y $b3 = - 0.5$.

Figure 1.3: Gráfica de perspectiva de una ecuación lineal con interacciones con \(a=0\), \(b_1 = 1\), \(b_2 = 1\) y \(b3 = - 0.5\).

Debido a la existencia de grupos discretos con medias “arbitrarias” dentro de cada variable independiente (\(x_1\) y \(x_2\)), las interacciones pueden tener efectos diversos dependiento de los niveles que interactúen. Por ejemplo, es posible que ciertas combinaciones de tratamientos de un estudio factorial tengan efectos negativos sobre ma media, mientras que otros incrementaron el valor promedio. En general, se considera que la interacciones son significativas cuando el efecto de un tratamiento depende de la presencia de otro tratamiento. Entonces, para interpretar correctamente las interacciones es necesario ver las diferencias entre grupos y los efectos estimados para cada combinación de tratamientos por el modelo estadístico.

2 Hipótesis estadísticas para interacciones

Si pensamos en una ANOVA de una vía, donde la variable dependiente \(y\) proviene de \(n\) tratamientos de un factor, la hipótesis nula establece que las medias de cada tratamiento \(i\) son:

\[ H_0 : \bar y_1 = \bar y_2 = \dots = \bar y_n\] En un ANOVA de dos vías aditivo, si las medias de cada tratamiento del primer factor las denotamos por \(\alpha_i\), las hipótesis nulas con respecto del efecto del primer factor independiente son:

\[ H_0 : \alpha_1 = \alpha_1 = \alpha_2 = \alpha_3 = \dots = \alpha_{n-1} = \alpha_{n} \] y las medias del segundo factor, denotadas por \(\beta_j\) son

\[H_1 : \beta_1 = \beta_2 = \beta_3 = \dots = \beta_{m-1} = \beta_m\]

En ambos casos la cantidad de efectos estimados corresponde con la cantidad de niveles que haya en cada factor. En as interacciones, el número total de efectos está determinado por el producto de la cantidad de niveles \(n \times m\). Si cada efecto estimado para las diferentes interacciones es \(\gamma_{i, j}\) tenemos que la hipótesis nula que se prueba es:

\[ H_2 : \gamma_{1, 1} = \gamma_{1, 2} = \dots = \gamma_{2, 1} = \dots = \gamma_{n, m}\]

Es por esta razon que en una tabla de ANOVA para un modelo de dos vías con interacciones se tienen tres valores de suma de cuadrados y de significancia, una para cada factor por separado y otra para cada interacción.

3 Ejemplos

Aquí volveremos a analizar la base de datos npk y construiremos el modelo probando varias combinaciones de las variables N, P y K. De entrada podemos identificar tres modelos con base en las diferentes combinaciones:

  1. yield ~ N * P
  2. yield ~ N * K
  3. yield ~ P * K

Para ahorrar código, ajustaremos los tres modelos con lm únicamente con efectos fijos (asumiremos para este análisis que no hay bloques experimentales), analizaremos sus respectivas sumas de cuadrados con \(SS-II\):

data(npk)
m1 <- lm(yield ~ N * P, data = npk)
m2 <- lm(yield ~ N * K, data = npk)
m3 <- lm(yield ~ P * K, data = npk)

car::Anova(m1, type = 2)
## Anova Table (Type II tests)
## 
## Response: yield
##           Sum Sq Df F value  Pr(>F)  
## N         189.28  1  5.7585 0.02627 *
## P           8.40  1  0.2556 0.61868  
## N:P        21.28  1  0.6474 0.43049  
## Residuals 657.40 20                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m2, type = 2)
## Anova Table (Type II tests)
## 
## Response: yield
##           Sum Sq Df F value  Pr(>F)  
## N         189.28  1  6.7752 0.01702 *
## K          95.20  1  3.4077 0.07975 .
## N:K        33.14  1  1.1860 0.28908  
## Residuals 558.75 20                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m3, type = 2)
## Anova Table (Type II tests)
## 
## Response: yield
##           Sum Sq Df F value Pr(>F)
## P           8.40  1  0.2176 0.6459
## K          95.20  1  2.4655 0.1321
## P:K         0.48  1  0.0125 0.9122
## Residuals 772.28 20

Para comenzar a comparar los modelos podemos ver cuál tiene la suma de cuadrados más pequeña para los residuales, es decir, en promedio cuál modelo predice mejor los datos pues estos están más cercanos a las verdaderas medias de los grupos. Con esto podemos ver que m2 con fómula yield ~ N * K tiene la menor suma de cuadrados de los residuales y también es el modelo en el cual la interacción tiene una \(SS\) más alta. Sin embargo, en todos los casos las interacciones no parecen afectar las medias de los grupos de las distintas combinaciones de uno y otro factor. En este caso los niveles de cada factor son 1 y 0, con lo que las comparaciones entre grupos con interacciones son sumamente sencillas, donde los grupos relevantes son:

Table 3.1: Combinaciones de tratamientos de las interacciones
N K
0 0
1 0
0 1
1 1

Si recordamos el significado de las interacciones en un modelo lineal como en la ecuación (1.1) es \(x_1 \times x_2\). Debido a que los valores de N y K son 0 y 1, los únicos valores posibles de la interacción son 0 y 1, es decir que sólo estamos haciendo la comparación entre los grupos que recibieron ambos tratamientos N y K. Viendo la tabla ANOVA de este modelo:

car::Anova(m2, type = 2)
## Anova Table (Type II tests)
## 
## Response: yield
##           Sum Sq Df F value  Pr(>F)  
## N         189.28  1  6.7752 0.01702 *
## K          95.20  1  3.4077 0.07975 .
## N:K        33.14  1  1.1860 0.28908  
## Residuals 558.75 20                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No podemos rechazar \(H_0\) de no hay diferencia entre los tratamientos que recibieron nitrógeno y potasio y los que sí los recibieron. Debido a que este modelo es el mejor pues sus residuales fueron menores, podemos proceder a evaluar el efecto de los bloques con un modelo de efectos mixtos:

library(lme4)
m4 <- lmer(yield ~ N * K + (1|block), data = npk)
car::Anova(m4, type = 2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: yield
##       Chisq Df Pr(>Chisq)    
## N   13.1780  1  0.0002833 ***
## K    6.6281  1  0.0100385 *  
## N:K  2.3069  1  0.1288009    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En este modelo, el factor de agrupación experimental ayudó a que los tratamientos experimentales explicaran más variabilidad, lo que disminuyó los valores de \(P\). Aún así podemos seguir rechazando \(H_0\) de no diferencia entre los tratamientos que sí recibieron N y K y los que no los recibieron. Con base en esto, sería sensible volver a ajustar un modelo sin la interacción:

m5 <- lmer(yield ~ N + K + (1|block), data = npk)
car::Anova(m5, type = 2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: yield
##     Chisq Df Pr(>Chisq)    
## N 12.1829  1  0.0004823 ***
## K  6.1275  1  0.0133091 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ahora sí, me parece que este es el mejor modelo, y concluimos, biológicamente, que el efecto del nitrógeno no depende de la aplciación de potasio y viceversa. Para saber si sus efectos son positivos o negativos veamos los coeficientes ajustados:

m5
## Linear mixed model fit by REML ['lmerMod']
## Formula: yield ~ N + K + (1 | block)
##    Data: npk
## REML criterion at convergence: 131.394
## Random effects:
##  Groups   Name        Std.Dev.
##  block    (Intercept) 3.644   
##  Residual             3.942   
## Number of obs: 24, groups:  block, 6
## Fixed Effects:
## (Intercept)           N1           K1  
##      54.058        5.617       -3.983

Con esto vemos que la aplicación de nitrógeno aumentó la producción (N1 = 5.617), y que la aplicación de potasio disminuyó la producción (K1 = - 3.983).

4 Aplicación del conocimiento

Dirígete a classroom para completar la actividad correspondiente.

Regresar al índice del curso