Saltar al contenido principal

# 16.6: ANOVA como modelo lineal

$$\newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$$

$$\newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$$

$$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$

( \newcommand{\kernel}{\mathrm{null}\,}\) $$\newcommand{\range}{\mathrm{range}\,}$$

$$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$

$$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}[1]{\| #1 \|}$$

$$\newcommand{\inner}[2]{\langle #1, #2 \rangle}$$

$$\newcommand{\Span}{\mathrm{span}}$$

$$\newcommand{\id}{\mathrm{id}}$$

$$\newcommand{\Span}{\mathrm{span}}$$

$$\newcommand{\kernel}{\mathrm{null}\,}$$

$$\newcommand{\range}{\mathrm{range}\,}$$

$$\newcommand{\RealPart}{\mathrm{Re}}$$

$$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$

$$\newcommand{\Argument}{\mathrm{Arg}}$$

$$\newcommand{\norm}[1]{\| #1 \|}$$

$$\newcommand{\inner}[2]{\langle #1, #2 \rangle}$$

$$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\AA}{\unicode[.8,0]{x212B}}$$

$$\newcommand{\vectorA}[1]{\vec{#1}} % arrow$$

$$\newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow$$

$$\newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$$

$$\newcommand{\vectorC}[1]{\textbf{#1}}$$

$$\newcommand{\vectorD}[1]{\overrightarrow{#1}}$$

$$\newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}}$$

$$\newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}}$$

$$\newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$$

$$\newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$$

$$\newcommand{\avec}{\mathbf a}$$ $$\newcommand{\bvec}{\mathbf b}$$ $$\newcommand{\cvec}{\mathbf c}$$ $$\newcommand{\dvec}{\mathbf d}$$ $$\newcommand{\dtil}{\widetilde{\mathbf d}}$$ $$\newcommand{\evec}{\mathbf e}$$ $$\newcommand{\fvec}{\mathbf f}$$ $$\newcommand{\nvec}{\mathbf n}$$ $$\newcommand{\pvec}{\mathbf p}$$ $$\newcommand{\qvec}{\mathbf q}$$ $$\newcommand{\svec}{\mathbf s}$$ $$\newcommand{\tvec}{\mathbf t}$$ $$\newcommand{\uvec}{\mathbf u}$$ $$\newcommand{\vvec}{\mathbf v}$$ $$\newcommand{\wvec}{\mathbf w}$$ $$\newcommand{\xvec}{\mathbf x}$$ $$\newcommand{\yvec}{\mathbf y}$$ $$\newcommand{\zvec}{\mathbf z}$$ $$\newcommand{\rvec}{\mathbf r}$$ $$\newcommand{\mvec}{\mathbf m}$$ $$\newcommand{\zerovec}{\mathbf 0}$$ $$\newcommand{\onevec}{\mathbf 1}$$ $$\newcommand{\real}{\mathbb R}$$ $$\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}$$ $$\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}$$ $$\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}$$ $$\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}$$ $$\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}$$ $$\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}$$ $$\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}$$ $$\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}$$ $$\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}$$ $$\newcommand{\laspan}[1]{\text{Span}\{#1\}}$$ $$\newcommand{\bcal}{\cal B}$$ $$\newcommand{\ccal}{\cal C}$$ $$\newcommand{\scal}{\cal S}$$ $$\newcommand{\wcal}{\cal W}$$ $$\newcommand{\ecal}{\cal E}$$ $$\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}$$ $$\newcommand{\gray}[1]{\color{gray}{#1}}$$ $$\newcommand{\lgray}[1]{\color{lightgray}{#1}}$$ $$\newcommand{\rank}{\operatorname{rank}}$$ $$\newcommand{\row}{\text{Row}}$$ $$\newcommand{\col}{\text{Col}}$$ $$\renewcommand{\row}{\text{Row}}$$ $$\newcommand{\nul}{\text{Nul}}$$ $$\newcommand{\var}{\text{Var}}$$ $$\newcommand{\corr}{\text{corr}}$$ $$\newcommand{\len}[1]{\left|#1\right|}$$ $$\newcommand{\bbar}{\overline{\bvec}}$$ $$\newcommand{\bhat}{\widehat{\bvec}}$$ $$\newcommand{\bperp}{\bvec^\perp}$$ $$\newcommand{\xhat}{\widehat{\xvec}}$$ $$\newcommand{\vhat}{\widehat{\vvec}}$$ $$\newcommand{\uhat}{\widehat{\uvec}}$$ $$\newcommand{\what}{\widehat{\wvec}}$$ $$\newcommand{\Sighat}{\widehat{\Sigma}}$$ $$\newcommand{\lt}{<}$$ $$\newcommand{\gt}{>}$$ $$\newcommand{\amp}{&}$$ $$\definecolor{fillinmathshade}{gray}{0.9}$$

Una de las cosas más importantes a entender sobre el ANOVA y la regresión es que básicamente son lo mismo. En la superficie de esto, no pensarías que esto es cierto: después de todo, la forma en que los he descrito hasta ahora sugiere que ANOVA se ocupa principalmente de probar las diferencias de grupo, y la regresión se ocupa principalmente de comprender las correlaciones entre variables. Y hasta donde va, eso es perfectamente cierto. Pero cuando miras bajo el capó, por así decirlo, la mecánica subyacente del ANOVA y la regresión son muy similares. De hecho, si lo piensas, ya has visto evidencias de ello. Tanto el ANOVA como la regresión dependen en gran medida de sumas de cuadrados (SS), ambos hacen uso de pruebas F, y así sucesivamente. Mirando hacia atrás, es difícil escapar de la sensación de que los Capítulos 14 y 15 fueron un poco repetitivos.

La razón de esto es que ANOVA y regresión son ambos tipos de modelos lineales. En el caso de la regresión, esto es algo obvio. La ecuación de regresión que usamos para definir la relación entre predictores y resultados es la ecuación para una línea recta, por lo que obviamente es un modelo lineal. Y si esa no fuera una pista lo suficientemente grande, el simple hecho de que el comando para ejecutar una regresión sea lm () también es una especie de pista. Cuando usamos una fórmula R como resultado ~ predictor1 + predictor2 con lo que realmente estamos trabajando es el modelo lineal algo más feo:

Y p =b 1 X 1p +b 2 X 2p +b 0 +p

donde Y p es el valor de resultado para la p-ésima observación (por ejemplo, p-ésima persona), X 1p es el valor del primer predictor para la p-ésima observación, X 2p es el valor del segundo predictor para la p-ésima observación, el b 1, b 2 y b 0 términos son nuestros coeficientes de regresión, y p es el p-ésimo residual. Si ignoramos los residuales p y solo nos enfocamos en la propia línea de regresión, obtenemos la siguiente fórmula:

$$\ \hat{Y_p}$$= b 1 X 1p +b 2 X 2p +b 0

donde$$\ \hat{Y_p}$$ está el valor de Y que la línea de regresión predice para la persona p, a diferencia del valor observado realmente Y p. Lo que no es inmediatamente obvio es que también podemos escribir ANOVA como modelo lineal. Sin embargo, en realidad es bastante sencillo hacer esto. Empecemos con un ejemplo realmente sencillo: reescribir un ANOVA factorial 2×2 como modelo lineal.

## Algunos datos

Para concretar las cosas, supongamos que nuestra variable de resultado es la calificación que recibe un alumno en mi clase, una variable ratio-escala correspondiente a una nota de 0% a 100%. Hay dos variables predictoras de interés: si el estudiante acudió o no a clases magistrales (la variable asistir), y si el estudiante realmente leyó o no el libro de texto (la variable de lectura). Diremos que asistir = 1 si el alumno asistió a clase, y asistir = 0 si no lo hizo. De igual manera, diremos que la lectura = 1 si el alumno leyó el libro de texto, y la lectura = 0 si no lo hizo.

Bien, hasta ahora eso es bastante simple. Lo siguiente que tenemos que hacer es envolver algunas matemáticas alrededor de esto (¡lo siento!). Para los efectos de este ejemplo, que Y p denote la calificación del alumno p-ésimo en la clase. Esta no es exactamente la misma notación que usamos anteriormente en este capítulo: anteriormente, hemos usado la notación Y rci para referirnos a la i-ésima persona en el grupo r-ésimo para el predictor 1 (el factor de fila) y el grupo c-ésimo para el predictor 2 (el factor de columna). Esta notación extendida fue muy útil para describir cómo se calculan los valores SS, pero es un dolor en el contexto actual, así que cambiaré la notación aquí. Ahora, la notación Yp es visualmente más simple que Y rci, ¡pero tiene la deficiencia de que en realidad no realiza un seguimiento de las membresías del grupo! Es decir, si te dijera que Y 0,0,3 =35, inmediatamente sabrías que estamos hablando de un estudiante (el 3er tal estudiante, de hecho) que no asistió a las clases (es decir, asiste = 0) y no leyó el libro de texto (es decir, leyendo = 0), y que terminó reprobando la clase ( grado = 35). Pero si te digo que Y p =35 todo lo que sabes es que el p-ésimo alumno no obtuvo una buena nota. Hemos perdido información clave aquí. Por supuesto, no hace falta mucho pensar en cómo solucionar esto: lo que haremos en cambio es introducir dos nuevas variables X 1p y X 2p que hagan un seguimiento de esta información. En el caso de nuestro estudiante hipotético, sabemos que X 1p =0 (es decir, asistir = 0) y X 2p =0 (es decir, lectura = 0). Entonces los datos podrían verse así:

knitr::kable(tibble::tribble(
~V1,   ~V2,     ~V3,     ~V4,
"1",  "90",     "1",     "1",
"2",  "87",     "1",     "1",
"3",  "75",     "0",     "1",
"4",  "60",     "1",     "0",
"5",  "35",     "0",     "0",
"6",  "50",     "0",     "0",
"7",  "65",     "1",     "0",
"8",  "70",     "0",     "1"),
col.names= c("person $p$", "grade $Y_p$", "attendance $X_{1p}$", "reading $X_{2p}$"), align = 'cccc')

persona p grado Y p asistencia X 1p lectura X 2p
5 35 0 0
6 50 0 0
4 60 1 0
7 65 1 0
8 70 0 1
3 75 0 1
2 87 1 1
1 90 1 1

Esto no es nada particularmente especial, claro: ¡es exactamente el formato en el que esperamos ver nuestros datos! En otras palabras, si sus datos se han almacenado como un marco de datos en R, entonces probablemente esté esperando ver algo que se parezca al marco de datos rtfm.1:

load("./rbook-master/data/rtfm.rdata")
rtfm.1
##   grade attend reading
## 1    90      1       1
## 2    87      1       1
## 3    75      0       1
## 4    60      1       0
## 5    35      0       0
## 6    50      0       0
## 7    65      1       0
## 8    70      0       1

Bueno, algo así como. Sospecho que algunos lectores probablemente estén frunciendo un poco el ceño en este punto. Anteriormente en el libro enfaticé la importancia de convertir variables de escala nominal como atender y leer a factores, en lugar de codificarlos como variables numéricas. El marco de datos rtfm.1 no hace esto, pero el marco de datos rtfm.2 sí, por lo que podría estar esperando ver datos como estos:

rtfm.2
##   grade attend reading
## 1    90    yes     yes
## 2    87    yes     yes
## 3    75     no     yes
## 4    60    yes      no
## 5    35     no      no
## 6    50     no      no
## 7    65    yes      no
## 8    70     no     yes

No obstante, para los efectos de esta sección es importante que podamos alternar de un lado a otro entre estas dos formas diferentes de pensar sobre los datos. Después de todo, nuestro objetivo en esta sección es mirar algunas de las matemáticas que sustentan el ANOVA, y si queremos hacerlo necesitamos poder ver la representación numérica de los datos (en rtfm.1) así como la representación factorial más significativa (en rtfm.2). En cualquier caso, podemos usar la función xtabs () para confirmar que este conjunto de datos corresponde a un diseño equilibrado

 xtabs( ~ attend + reading, rtfm.2 )


##       reading
## attend no yes
##    no   2   2
##    yes  2   2

Por cada posible combinación de las variables de asistencia y lectura, tenemos exactamente dos alumnos. Si nos interesa calcular la calificación media para cada una de estas celdas, podemos usar la función aggregate ():

aggregate( grade ~ attend + reading, rtfm.2, mean )
##   attend reading grade
## 1     no      no  42.5
## 2    yes      no  62.5
## 3     no     yes  72.5
## 4    yes     yes  88.5

Al mirar esta mesa, uno tiene la fuerte impresión de que leer el texto y asistir a la clase importan mucho.

## ANOVA con factores binarios como modelo de regresión

Bien, volvamos a hablar de las matemáticas. Ahora tenemos nuestros datos expresados en términos de tres variables numéricas: la variable continua Y, y las dos variables binarias X1 y X2. Lo que quiero que reconozcas es que nuestro ANOVA factorial de 22 es exactamente equivalente al modelo de regresión

Y p =b 1 X 1p +b 2 X 2p +b 0 +p

¡Esta es, por supuesto, exactamente la misma ecuación que usé antes para describir un modelo de regresión de dos predictores! La única diferencia es que X1 y X2 son ahora variables binarias (es decir, los valores sólo pueden ser 0 o 1), mientras que en un análisis de regresión esperamos que X1 y X2 sean continuos. Hay un par de formas en las que podría tratar de convencerte de esto. Una posibilidad sería hacer un largo ejercicio matemático, demostrando que los dos son idénticos. Sin embargo, voy a salir en una extremidad y adivinar que la mayoría de los lectores de este libro encontrarán que eso es molesto en lugar de útil. En cambio, explicaré las ideas básicas, y luego confiaré en R para demostrar que los análisis ANOVA y los análisis de regresión no son solo similares, son idénticos a todos los efectos. 239 Empecemos por ejecutar esto como un ANOVA. Para ello, usaremos el marco de datos rtfm.2, ya que ese es el que hice lo correcto de codificar asistir y leer como factores, y usaré la función aov () para hacer el análisis. Esto es lo que obtenemos...

anova.model <- aov( grade ~ attend + reading, data = rtfm.2 )
summary( anova.model ) 
##             Df Sum Sq Mean Sq F value  Pr(>F)
## attend       1    648     648   21.60 0.00559 **
## reading      1   1568    1568   52.27 0.00079 ***
## Residuals    5    150      30
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Entonces, al leer los números clave de la tabla ANOVA y la tabla de medios que presentamos anteriormente, podemos ver que los estudiantes obtuvieron una calificación superior si asistían a clase (F 1,5 =26.1, p=.0056) y si leyeron el libro de texto (F 1,5 =52.3, p=.0008). Tomemos nota de esos valores p y esas estadísticas de F.

library(effects)
Effect( c("attend","reading"), anova.model )
##
## attend   no  yes
##    no  43.5 71.5
##    yes 61.5 89.5

Ahora pensemos en el mismo análisis desde una perspectiva de regresión lineal. En el conjunto de datos rtfm.1, hemos codificado atender y leer como si fueran predictores numéricos. En este caso, esto es perfectamente aceptable. Realmente hay un sentido en el que un estudiante que acude a clase (es decir, asiste = 1) de hecho ha hecho “más asistencia” que un estudiante que no (es decir, asiste = 0). Por lo que no es nada irrazonable incluirlo como predictor en un modelo de regresión. Es un poco inusual, porque el predictor sólo toma dos valores posibles, pero no viola ninguno de los supuestos de regresión lineal. Y es fácil de interpretar. Si el coeficiente de regresión para asistir es mayor a 0, significa que los estudiantes que asisten a clases obtienen calificaciones más altas; si es menor que cero, entonces los estudiantes que asisten a clases obtienen calificaciones más bajas. Lo mismo es cierto para nuestra variable de lectura.

Espera un segundo... ¿por qué es esto cierto? Es algo que es intuitivamente obvio para todos los que han tomado algunas clases de estadísticas y se sienten cómodos con las matemáticas, pero no es claro para todos los demás en el primer pase. Para ver por qué esto es cierto, ayuda mirar de cerca a algunos estudiantes específicos. Empecemos considerando a los alumnos 6º y 7º en nuestro conjunto de datos (es decir, p=6 y p=7). Ninguno ha leído el libro de texto, así que en ambos casos podemos establecer la lectura = 0. O, para decir lo mismo en nuestra notación matemática, observamos X 2,6 =0 y X 2,7 =0. Sin embargo, el estudiante número 7 sí acudió a clases magistrales (es decir, asiste = 1, X 1,7 =1) mientras que el estudiante número 6 no lo hizo (es decir, asiste = 0, X 1,6 =0). Ahora veamos qué sucede cuando insertamos estos números en la fórmula general para nuestra línea de regresión. Para el estudiante número 6, la regresión predice que

\begin{aligned} Y_6 &= b_1X_{1,6} + b_2X_{2,6} + b_0 \\& = (b_1 \times 0) + (b_2 \times 0) + b_0 \\&=b_0 \end{aligned}

Por lo que estamos esperando que este alumno obtenga una calificación correspondiente al valor del término de intercepción b0. ¿Qué pasa con el estudiante 7? Esta vez, cuando insertamos los números en la fórmula para la línea de regresión, obtenemos lo siguiente:

\begin{aligned} Y_6 &= b_1X_{1,7} + b_2X_{2,7} + b_0 \\& = (b_1 \times 1) + (b_2 \times 0) + b_0 \\&= b_1 + b_0 \end{aligned}

Debido a que este estudiante asistió a clase, la calificación prevista es igual al término de intercepción b 0 más el coeficiente asociado a la variable asistir, b 1. Entonces, si b 1 es mayor que cero, estamos esperando que los alumnos que acuden a clases obtengan calificaciones más altas que aquellos que no lo hacen.Si este coeficiente es negativo, estamos esperando lo contrario: los estudiantes que aparecen en clase terminen desempeñándose mucho peor. De hecho, podemos empujar esto un poco más allá. ¿Qué pasa con el estudiante número 1, que acudió a clase (X 1,1 =1) y leyó el libro de texto (X 2,1 =1)? Si enchufamos estos números en la regresión, obtenemos

\begin{aligned} Y_6 &= b_1X_{1,1} + b_2X_{2,1} + b_0 \\& = (b_1 \times 1) + (b_2 \times 1) + b_0 \\&= b_1 + b_2 + b_0 \end{aligned}

Entonces, si asumimos que asistir a clase te ayuda a obtener una buena nota (es decir, b1>0) y si asumimos que leer el libro de texto también te ayuda a obtener una buena calificación (es decir, b 2 >0), entonces nuestra expectativa es que el estudiante 1 obtenga una calificación que sea superior a la del estudiante 6 y al estudiante 7.

Y en este punto, no te sorprenderá en absoluto al saber que el modelo de regresión predice que el estudiante 3, que leyó el libro pero no asistió a conferencias, obtendrá una calificación de b 2 +b 0. No te aburriré con otra fórmula de regresión más. En cambio, lo que voy a hacer es mostrarte la siguiente tabla de calificaciones esperadas:

knitr::kable(tibble::tribble(
~V1,           ~V2,   ~V3,

"attended - no","$b_0$","$b_0 + b_2$",
"attended - yes", "$b_0 + b_1$", "$b_0 + b_1 + b_2$"),
col.names = c("","read textbook? no", "read textbook? yes"))
leer libro de texto? no leer libro de texto? si
atendido - no b 0 b 0 +b 2
atendido - sí b 0 +b 1 b 0 +b 1 +b 2

Como puedes ver, el término de intercepción b 0 actúa como una especie de calificación “base” que esperarías de aquellos alumnos que no se toman el tiempo para asistir a clase o leer el libro de texto. De igual manera, b 1 representa el impulso que se espera que obtengas si vienes a clase, y b 2 representa el impulso que viene de leer el libro de texto. De hecho, si esto fuera un ANOVA muy bien podrías querer caracterizar a b 1 como el efecto principal de la asistencia, ¡y b 2 como el efecto principal de la lectura! De hecho, para un simple ANOVA 2×2 así es exactamente como se desarrolla.

Bien, ahora que realmente estamos empezando a ver por qué ANOVA y regresión son básicamente lo mismo, vamos a ejecutar nuestra regresión usando los datos rtfm.1 y la función lm () para convencernos de que esto es realmente cierto. Ejecutar la regresión de la manera habitual nos da el siguiente resultado: 240

 regression.model <- lm( grade ~ attend + reading, data = rtfm.1 )
summary( regression.model )
##
## Call:
##
## Residuals:
##    1    2    3    4    5    6    7    8
##  0.5 -2.5  3.5 -1.5 -8.5  6.5  3.5 -1.5
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   43.500      3.354  12.969 4.86e-05 ***
## attend        18.000      3.873   4.648  0.00559 **
## reading       28.000      3.873   7.230  0.00079 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.477 on 5 degrees of freedom
## Multiple R-squared:  0.9366, Adjusted R-squared:  0.9112
## F-statistic: 36.93 on 2 and 5 DF,  p-value: 0.001012

Aquí hay algunas cosas interesantes a tener en cuenta. En primer lugar, observe que el término de interceptación es de 43.5, que se acerca a la media “grupal” de 42.5 observada para aquellos dos alumnos que no leyeron el texto ni asistieron a clase. Además, es idéntico al grupo predicho, significa que sacamos de nuestro ANOVA usando la función Effects ()! En segundo lugar, observe que tenemos el coeficiente de regresión de b 1 =18.0 para la variable de asistencia, sugiriendo que aquellos estudiantes que asistieron a clase obtuvieron un 18% más alto que los que no lo hicieron. entonces nuestra expectativa sería que aquellos estudiantes que acudieron a clase pero no leyeron el libro de texto lo harían obtener una calificación de b 0 +b 1, que es igual a 43.5+18.0=61.5. Nuevamente, esto es similar a la media observada del grupo de 62.5, e idéntica a la media esperada del grupo que sacamos de nuestro ANOVA. Puedes verificar por ti mismo que sucede lo mismo cuando miramos a los alumnos que leen el libro de texto.

En realidad, podemos empujar un poco más para establecer la equivalencia de nuestro ANOVA y nuestra regresión. Observe los valores p asociados con la variable attend y la variable de lectura en el resultado de regresión. Son idénticos a los que encontramos antes al ejecutar el ANOVA. Esto puede parecer un poco sorprendente, ya que la prueba utilizada al ejecutar nuestro modelo de regresión calcula un estadístico t y el ANOVA calcula un estadístico F. Sin embargo, si recuerdas todo el camino de regreso al Capítulo 9, mencioné que existe una relación entre la distribución t y la distribución F: si tienes alguna cantidad que se distribuye de acuerdo con una distribución t con k grados de libertad y la cuadras, entonces esta nueva cantidad al cuadrado sigue una Distribución F cuyos grados de libertad son 1 y k, lo podemos verificar con respecto a la estadística t en nuestro modelo de regresión. Para la variable attend obtenemos un valor t de 4.648. Si cuadramos este número terminamos con 21.604, que es idéntico al estadístico F correspondiente en nuestro ANOVA.

Por último, una última cosa que debes saber. Debido a que R entiende el hecho de que tanto ANOVA como regresión son ejemplos de modelos lineales, le permite extraer la tabla ANOVA clásica de su modelo de regresión usando la función anova (). Todo lo que tienes que hacer es esto:

 anova( regression.model )
## Analysis of Variance Table
##
##           Df Sum Sq Mean Sq F value    Pr(>F)
## attend     1    648     648  21.600 0.0055943 **
## reading    1   1568    1568  52.267 0.0007899 ***
## Residuals  5    150      30
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Cambiar la categoría de línea base

En este punto, probablemente estés convencido de que el ANOVA y la regresión son en realidad idénticos entre sí. Entonces hay una última cosa que debería mostrarte. ¿Qué pasa si utilizo los datos de rtfm.2 para ejecutar la regresión? En rtfm.2, codificamos las variables de asistencia y lectura como factores y no como variables numéricas. ¿Importa esto? Resulta que no lo hace Las únicas diferencias son superficiales:

regression.model.2 <- lm( grade ~ attend + reading, data = rtfm.2 )
summary( regression.model.2 )
##
## Call:
##
## Residuals:
##    1    2    3    4    5    6    7    8
##  0.5 -2.5  3.5 -1.5 -8.5  6.5  3.5 -1.5
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   43.500      3.354  12.969 4.86e-05 ***
## attendyes     18.000      3.873   4.648  0.00559 **
## readingyes    28.000      3.873   7.230  0.00079 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.477 on 5 degrees of freedom
## Multiple R-squared:  0.9366, Adjusted R-squared:  0.9112
## F-statistic: 36.93 on 2 and 5 DF,  p-value: 0.001012

Lo único que es diferente es que R etiqueta las dos variables de manera diferente: la salida ahora se refiere a attendyes y readingyes. Probablemente puedas adivinar lo que esto significa. Cuando R se refiere a leeryes se trata de indicar que está asumiendo que “yes = 1” y “no = 0”. Esto es importante. Supongamos que quisiéramos decir que “sí = 0” y “no = 1”. Todavía podríamos ejecutar esto como un modelo de regresión, pero ahora todos nuestros coeficientes irán en sentido contrario, porque el efecto de leerno estaría refiriéndose a las consecuencias de no leer el libro de texto. Para mostrarte cómo funciona esto, podemos usar la función relevel () en R para cambiar qué nivel de la variable de lectura se establece en “0”. Así es como funciona. Primero, hagamos que R imprima el factor de lectura tal como está actualmente:

rtfm.2$reading ## [1] yes yes yes no no no no yes ## Levels: no yes Observe que el orden en el que R imprime los niveles es “no” y luego “sí”. Ahora vamos a aplicar la función relevel ():  relevel( x = rtfm.2$reading, ref = "yes" )
## [1] yes yes yes no  no  no  no  yes
## Levels: yes no

R ahora enumera “sí” antes que “no”. Esto significa que R ahora tratará al “sí” como el nivel de “referencia” (a veces llamado nivel basal) cuando lo incluyas en un ANOVA. Así que ahora vamos a crear un nuevo marco de datos con nuestros factores recodificados...

rtfm.3 <- rtfm.2                                        # copy the old data frame
rtfm.3$reading <- relevel( rtfm.2$reading, ref="yes" )  # re-level the reading factor
rtfm.3$attend <- relevel( rtfm.2$attend, ref="yes" )    # re-level the attend factor

Finalmente, volvamos a ejecutar nuestra regresión, esta vez usando los datos re-codificados:

regression.model.3 <- lm( grade ~ attend + reading, data = rtfm.3 )
summary( regression.model.3 )

##
## Call:
##
## Residuals:
##    1    2    3    4    5    6    7    8
##  0.5 -2.5  3.5 -1.5 -8.5  6.5  3.5 -1.5
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   89.500      3.354  26.684 1.38e-06 ***
## attendno     -18.000      3.873  -4.648  0.00559 **
## readingno    -28.000      3.873  -7.230  0.00079 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.477 on 5 degrees of freedom
## Multiple R-squared:  0.9366, Adjusted R-squared:  0.9112
## F-statistic: 36.93 on 2 and 5 DF,  p-value: 0.001012

Como puede ver, ahora hay algunos cambios. Lo más obvio es que los efectos attendno y readingno son ambos negativos, aunque son de la misma magnitud que antes: si no lees el libro de texto, por ejemplo, debes esperar que tu calificación baje 28% en relación con alguien que lo hizo. Las estadísticas t también han invertido el signo. Los valores p siguen siendo los mismos, por supuesto. La intercepción también ha cambiado. En nuestra regresión original, la línea de base correspondía a un estudiante que no asistió a clase y no leyó el libro de texto, por lo que obtuvimos un valor de 43.5 como nota basal esperada. No obstante, ahora que hemos recodificado nuestras variables, la línea de base corresponde a un estudiante que ha leído el libro de texto y ha asistido a clase, y para ese estudiante esperaríamos una calificación de 89.5.

## codificar factores no binarios como contrastes

En este punto, te he mostrado cómo podemos ver un ANOVA 2×2 en un modelo lineal. Y es bastante fácil ver como esto generaliza a un ANOVA de 2×2×2 o un ANOVA de 2×2×2×2... es lo mismo, realmente: solo agregas una nueva variable binaria para cada uno de tus factores. Donde empieza a ponerse más complicado es cuando consideramos factores que tienen más de dos niveles. Consideremos, por ejemplo, el ANOVA 3×2 que ejecutamos anteriormente en este capítulo usando los datos de clin.trial. ¿Cómo podemos convertir el factor de fármaco de tres niveles en una forma numérica apropiada para una regresión?

La respuesta a esta pregunta es bastante simple, en realidad. Todo lo que tenemos que hacer es darnos cuenta de que un factor de tres niveles puede ser redescrito como dos variables binarias. Supongamos, por ejemplo, que iba a crear una nueva variable binaria llamada druganxifree. Siempre que la variable droga sea igual a “anxifree” establecemos druganxifree = 1. De lo contrario, establecemos druganxifree = 0. Esta variable establece un contraste, en este caso entre ansifree y las otras dos drogas. Por sí mismo, por supuesto, el contraste libre de drogaxias no es suficiente para capturar completamente toda la información en nuestra variable de medicamento. Necesitamos un segundo contraste, uno que nos permita distinguir entre joyzepam y el placebo. Para ello, podemos crear un segundo contraste binario, llamado farmjoyzepam, que equivale a 1 si el medicamento es joyzepam, y 0 si no lo es. Tomados en conjunto, estos dos contrastes nos permiten discriminar perfectamente entre las tres drogas posibles. La siguiente tabla ilustra esto:

knitr::kable(tibble::tribble(
~V1,              ~V2,              ~V3,
"placebo",              "0",              "0",
"anxifree",              "1",              "0",
"joyzepam",              "0",              "1"
), col.names = c(      "drug", "druganxifree", "drugjoyzepam"))
droga druganxifree drogajoyzepam
placebo 0 0
ansioso 1 0
joyzepam 0 1

Si el fármaco administrado a un paciente es un placebo, entonces ambas variables de contraste serán iguales a 0. Si el medicamento es Anxifree, entonces la variable droganxifree será igual a 1, y farmjoyzepam será 0. Lo contrario es cierto para Joyzepam: drogajoyzepam es 1, y droganxilibre es 0.

Crear variables de contraste manualmente no es demasiado difícil de hacer usando comandos R básicos. Por ejemplo, así es como crearíamos la variable druganxifree:

druganxifree <- as.numeric( clin.trial$drug == "anxifree" ) druganxifree  ## [1] 0 0 0 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 El clin.trial$drug == parte “anxifree” del comando devuelve un vector lógico que tiene un valor de TRUE si el medicamento es Anxifree, y un valor de FALSO si el medicamento es Joyzepam o el placebo. La función.numeric () solo convierte TRUE a 1 y FALSE a 0. Obviamente, este comando crea la variable druganxifree dentro del espacio de trabajo. Si quisieras agregarlo al marco de datos clin.trial, usarías un commmand como este en su lugar:

clin.trial$druganxifree <- as.numeric( clin.trial$drug == "anxifree" )

Entonces podrías repetir esto para los otros contrastes que quisieras usar. Sin embargo, es un poco tedioso hacer esto una y otra vez por cada contraste que quieras crear. Para hacerlo un poco más fácil, el paquete lsr contiene una función simple llamada expandFactors () que convertirá cada factor en un marco de datos en un conjunto de variables de contraste. 241 Podemos usarlo para crear un nuevo marco de datos, clin.trial.2 que contenga los mismos datos que clin.trial, pero con los dos factores representados en términos de las variables de contraste:

 clin.trial.2 <- expandFactors( clin.trial )
##    (Intercept) druganxifree drugjoyzepam therapyCBT mood.gain druganxifree
## 1            1            0            0          0       0.5            0
## 2            1            0            0          0       0.3            0
## 3            1            0            0          0       0.1            0
## 4            1            1            0          0       0.6            1
## 5            1            1            0          0       0.4            1
## 6            1            1            0          0       0.2            1
## 7            1            0            1          0       1.4            0
## 8            1            0            1          0       1.7            0
## 9            1            0            1          0       1.3            0
## 10           1            0            0          1       0.6            0
## 11           1            0            0          1       0.9            0
## 12           1            0            0          1       0.3            0
## 13           1            1            0          1       1.1            1
## 14           1            1            0          1       0.8            1
## 15           1            1            0          1       1.2            1
## 16           1            0            1          1       1.8            0
## 17           1            0            1          1       1.3            0
## 18           1            0            1          1       1.4            0
## attr(,"assign")
## [1] 0 1 1 2 3 4
## attr(,"contrasts")
## attr(,"contrasts")$drug ## [1] "contr.treatment" ## ## attr(,"contrasts")$therapy
## [1] "contr.treatment"
 clin.trial.2
##    druganxifree drugjoyzepam therapyCBT mood.gain druganxifree
## 1             0            0          0       0.5            0
## 2             0            0          0       0.3            0
## 3             0            0          0       0.1            0
## 4             1            0          0       0.6            1
## 5             1            0          0       0.4            1
## 6             1            0          0       0.2            1
## 7             0            1          0       1.4            0
## 8             0            1          0       1.7            0
## 9             0            1          0       1.3            0
## 10            0            0          1       0.6            0
## 11            0            0          1       0.9            0
## 12            0            0          1       0.3            0
## 13            1            0          1       1.1            1
## 14            1            0          1       0.8            1
## 15            1            0          1       1.2            1
## 16            0            1          1       1.8            0
## 17            0            1          1       1.3            0
## 18            0            1          1       1.4            0

No es tan bonito como los datos originales de clin.trial, pero definitivamente está diciendo lo mismo. Ahora hemos recodificado nuestro factor de tres niveles en términos de dos variables binarias, y ya hemos visto que el ANOVA y la regresión se comportan de la misma manera para las variables binarias. No obstante, existen algunas complejidades adicionales que surgen en este caso, que discutiremos en la siguiente sección.

## equivalencia entre ANOVA y regresión para factores no binarios

Ahora tenemos dos versiones diferentes del mismo conjunto de datos: nuestro marco de datos original clin.trial en el que la variable de fármaco se expresa como un solo factor de tres niveles, y el conjunto de datos expandido clin.trial.2 en el que se expande en dos contrastes binarios. Una vez más, lo que queremos demostrar es que nuestro ANOVA factorial 3×2 original es equivalente a un modelo de regresión aplicado a las variables de contraste. Empecemos por volver a ejecutar el ANOVA:

 drug.anova <- aov( mood.gain ~ drug + therapy, clin.trial )
summary( drug.anova )
##             Df Sum Sq Mean Sq F value   Pr(>F)
## drug         2  3.453  1.7267  26.149 1.87e-05 ***
## therapy      1  0.467  0.4672   7.076   0.0187 *
## Residuals   14  0.924  0.0660
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Obviamente, aquí no hay sorpresas. Ese es exactamente el mismo ANOVA que ejecutamos antes, excepto por el hecho de que arbitrariamente he decidido cambiar el nombre de la variable de salida como drug.anova por alguna estúpida razón. 242 A continuación, hagamos una regresión, usando druganxifree, drug joyzepam y TherapyCBT como predictores. Esto es lo que obtenemos:

drug.regression <- lm( mood.gain ~ druganxifree + drugjoyzepam + therapyCBT, clin.trial.2 )
summary( drug.regression )
##
## Call:
## lm(formula = mood.gain ~ druganxifree + drugjoyzepam + therapyCBT,
##     data = clin.trial.2)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -0.3556 -0.1806  0.0000  0.1972  0.3778
##
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)    0.2889     0.1211   2.385   0.0318 *
## druganxifree   0.2667     0.1484   1.797   0.0939 .
## drugjoyzepam   1.0333     0.1484   6.965  6.6e-06 ***
## therapyCBT     0.3222     0.1211   2.660   0.0187 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.257 on 14 degrees of freedom
## Multiple R-squared:  0.8092, Adjusted R-squared:  0.7683
## F-statistic: 19.79 on 3 and 14 DF,  p-value: 2.64e-05

Hm. Esta no es la misma salida que obtuvimos la última vez. No es sorprendente que la salida de regresión imprima los resultados de cada uno de los tres predictores por separado, tal como lo hacía cada dos veces que usamos lm (). Por un lado, podemos ver que el valor p para la variable terapiaCBT es exactamente el mismo que el del factor de terapia en nuestro ANOVA original, por lo que podemos estar seguros de que el modelo de regresión está haciendo lo mismo que el ANOVA. Por otro lado, este modelo de regresión está probando el contraste libre de druganxix y el contraste de fármaco-joyzepam por separado, como si fueran dos variables completamente no relacionadas. Por supuesto, no es sorprendente, porque la pobre función lm () no tiene forma de saber que drogajoyzepam y druganxifree son en realidad los dos contrastes diferentes que usamos para codificar nuestro factor de fármaco de tres niveles. Hasta donde se sabe, drogajoyzepam y droganxifree no están más relacionados entre sí que drogajoyzepam y terapiaCBT. Sin embargo, usted y yo sabemos mejor. En esta etapa no nos interesa en absoluto determinar si estos dos contrastes son individualmente significativos. Sólo queremos saber si hay un efecto “general” de la droga. Es decir, lo que queremos que haga R es ejecutar algún tipo de prueba “ómnibus”, una en la que se agrupen los dos contrastes “relacionados con las drogas” para el propósito de la prueba. ¿Te suena familiar? Esta es exactamente la situación que discutimos en la Sección 16.5, y es precisamente esta situación la que la prueba F está construida para manejar. Todo lo que necesitamos hacer es especificar nuestro modelo nulo, que en este caso incluiría el predictor TherapyCBT, y omitir ambas variables relacionadas con el fármaco, y luego ejecutarlo a través de la función anova ():

nodrug.regression <- lm( mood.gain ~ therapyCBT, clin.trial.2 )
anova( nodrug.regression, drug.regression )
## Analysis of Variance Table
##
## Model 1: mood.gain ~ therapyCBT
## Model 2: mood.gain ~ druganxifree + drugjoyzepam + therapyCBT
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)
## 1     16 4.3778
## 2     14 0.9244  2    3.4533 26.149 1.872e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ah, eso es mejor. Nuestro estadístico F es 26.1, los grados de libertad son 2 y 14, y el valor p es 0.000019. Los números son idénticos a los que obtuvimos para el efecto principal del fármaco en nuestro ANOVA original. Una vez más, vemos que ANOVA y regresión son esencialmente los mismos: ambos son modelos lineales, y la maquinaria estadística subyacente para ANOVA es idéntica a la maquinaria utilizada en regresión. No se debe subestimar la importancia de este hecho. A lo largo del resto de este capítulo vamos a confiar mucho en esta idea.

Por fin, por fin puedo dar una definición de grados de libertad con los que estoy contento. Los grados de libertad se definen en términos del número de parámetros que se tienen que estimar en un modelo. Para un modelo de regresión o un ANOVA, el número de parámetros corresponde al número de coeficientes de regresión (es decir, valores b), incluyendo la intercepción. Teniendo en cuenta que cualquier prueba F es siempre una comparación entre dos modelos, el primer df es la diferencia en el número de parámetros. Por ejemplo, la comparación del modelo anterior, el modelo nulo (mood.gain ~ therapyCBT) tiene dos parámetros: hay un coeficiente de regresión para la variable TherapyCBT, y un segundo para la intercepción. El modelo alternativo (mood.gain ~ druganxifree + drug joyzepam + therapyCBT) tiene cuatro parámetros: un coeficiente de regresión para cada uno de los tres contrastes y uno más para el intercepto. Entonces los grados de libertad asociados a la diferencia entre estos dos modelos es df 1 =4−2=2.

¿Qué pasa con el caso cuando no parece haber un modelo nulo? Por ejemplo, podrías estar pensando en la prueba F que aparece en la parte inferior de la salida de regresión. Originalmente lo describí como una prueba del modelo de regresión en su conjunto. Sin embargo, eso sigue siendo una comparación entre dos modelos. El modelo nulo es el modelo trivial que solo incluye una intercepción, la cual se escribe como resultado ~ 1 en R, y el modelo alternativo es el modelo de regresión completa. El modelo nulo en este caso contiene 1 coeficiente de regresión, para el término de intercepción. El modelo alternativo contiene coeficientes de regresión K+1, uno para cada una de las K variables predictoras y otro más para la intercepción. Entonces el valor df que ves en esta prueba F es igual a df 1 =K+1−1=K.

## posdata

Hay una última cosa que quiero mencionar en esta sección. En el ejemplo anterior, utilicé la función aov () para ejecutar un ANOVA usando los datos clin.trial que codifica la variable de fármaco como un solo factor. También utilicé la función lm () para ejecutar una regresión usando los datos de clin.trial en los que tenemos dos contrastes separados que describen el fármaco. Sin embargo, también es posible utilizar la función lm () en los datos originales. Es decir, podrías usar un comando como este:

drug.lm <- lm( mood.gain ~ drug + therapy, clin.trial )

El hecho de que la droga sea un factor de tres niveles no importa. Siempre y cuando se haya declarado que la variable fármaco es un factor, R la traducirá automáticamente en dos variables de contraste binario, y realizará el análisis apropiado. Después de todo, como he estado diciendo a lo largo de esta sección, ANOVA y regresión son modelos lineales, y lm () es la función que maneja los modelos lineales. De hecho, la función aov () en realidad no hace mucho del trabajo cuando ejecutas un ANOVA usándolo: internamente, R simplemente pasa todo el trabajo duro directamente a lm (). No obstante, quiero recalcar nuevamente que es crítico que sus variables factoriales sean declaradas como tales. Si el fármaco se declarara como una variable numérica, entonces R estaría feliz de tratarlo como uno solo. Después de todo, podría ser que droga se refiere al número de medicamentos que uno ha tomado en el pasado, o algo que es genuinamente numérico. R no va a dudar de ti aquí. Asume que tus factores son factores y tus números son números. No cometas el error de codificar tus factores como números, o R ejecutará el análisis incorrecto. Esto no es un defecto en R: es su responsabilidad como analista asegurarse de que está especificando el modelo adecuado para sus datos. Realmente no se puede confiar en el software con este tipo de cosas.

Bien, a un lado las advertencias, en realidad es un poco ordenado ejecutar tu ANOVA usando la función lm () de la manera que hice arriba. Debido a que has llamado a la función lm (), el summary () que R saca está formateado como una regresión. Para ahorrar espacio no te mostraré la salida aquí, pero puedes verificar esto fácilmente escribiendo

 summary( drug.lm )
##
## Call:
## lm(formula = mood.gain ~ drug + therapy, data = clin.trial)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -0.3556 -0.1806  0.0000  0.1972  0.3778
##
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)    0.2889     0.1211   2.385   0.0318 *
## druganxifree   0.2667     0.1484   1.797   0.0939 .
## drugjoyzepam   1.0333     0.1484   6.965  6.6e-06 ***
## therapyCBT     0.3222     0.1211   2.660   0.0187 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.257 on 14 degrees of freedom
## Multiple R-squared:  0.8092, Adjusted R-squared:  0.7683
## F-statistic: 19.79 on 3 and 14 DF,  p-value: 2.64e-05

Sin embargo, debido a que las variables de fármaco y terapia fueron ambos factores, la función anova () en realidad sabe qué contrasta agruparse para los fines de ejecutar las pruebas F, por lo que se puede extraer la tabla ANOVA clásica. Nuevamente, no voy a reproducir aquí la salida ya que es idéntica a la tabla ANOVA que mostré al inicio de la sección, pero vale la pena probar el siguiente comando

 anova( drug.lm )
## Analysis of Variance Table
##
## Response: mood.gain
##           Df Sum Sq Mean Sq F value    Pr(>F)
## drug       2 3.4533 1.72667 26.1490 1.872e-05 ***
## therapy    1 0.4672 0.46722  7.0757   0.01866 *
## Residuals 14 0.9244 0.06603
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sólo para ver por ti mismo. Sin embargo, este comportamiento de la función anova () solo ocurre cuando las variables predictoras son factores. Si intentamos un comando como anova (drug.regression), la salida continuará tratando druganxifree y drug joyzepam como si fueran dos factores binarios distintos. Esto se debe a que en el modelo fármaco.regresión incluimos todos los contrastes como variables “crudas”, por lo que R no tenía idea de cuáles pertenecían juntas. Sin embargo, cuando ejecutamos el modelo drug.lm, le dimos a R las variables factoriales originales, por lo que sí sabe qué contrastes van juntos. El comportamiento de la función anova () refleja eso.

This page titled 16.6: ANOVA como modelo lineal is shared under a CC BY-SA 4.0 license and was authored, remixed, and/or curated by Danielle Navarro via source content that was edited to the style and standards of the LibreTexts platform.