ANCOVA como GLM

Maestría en Ciencias Actuariales Aplicadas

1 Introducción

En la clase anterior se estudió ANOVA como una técnica para comparar medias. Esa descripción es correcta, pero incompleta. En esta sesión se dará un paso conceptual más importante: ANOVA será reinterpretado como un modelo lineal.

Esta reinterpretación permitirá comprender que:

• la regresión, ANOVA y ANCOVA no son técnicas aisladas;
• las tres comparten la misma estructura matemática;
• ANCOVA incorpora simultáneamente predictores categóricos y cuantitativos;
• los efectos fijos y aleatorios responden preguntas inferenciales distintas;
• el Modelo Lineal General ofrece un lenguaje común para organizar todos estos análisis.

Idea central del modelo

\[ \text{Datos}=\text{Modelo}+\text{Error} \]

ANOVA, regresión y ANCOVA se distinguen principalmente por la forma en que se construye el componente denominado modelo.

1.1 Objetivos que se prentenden alcanzar

Al finalizar la sesión, debemos ser capaces de:

1. Explicar ANOVA como un caso particular del Modelo Lineal General.
2. Mostrar la equivalencia entre ANOVA y una regresión con variables indicadoras.
3. Diferenciar predictores categóricos y cuantitativos dentro de una misma matriz de diseño.
4. Formular ANCOVA como la combinación de ANOVA y regresión.
5. Interpretar medias observadas y medias ajustadas.
6. Verificar los principales supuestos de ANCOVA (clase de supuestos).
7. Distinguir correctamente efectos fijos, aleatorios y mixtos.
8. Reconocer cuándo debe emplearse un modelo lineal mixto.
9. Comparar modelos completos y reducidos.
10. Establecer el puente hacia el estudio formal del Modelo Lineal General (próxima clase).

2 Preparación del entorno de trabajo

El documento utiliza paquetes disponibles en CRAN.

install.packages(
  c(
    "ggplot2",
    "dplyr",
    "tidyr",
    "car",
    "emmeans",
    "lmtest",
    "lme4",
    "lmerTest"
  )
)

3 ¿Qué es realmente ANOVA?

ANOVA suele presentarse como una prueba para determinar si varias medias son iguales. Sin embargo, desde la perspectiva del modelamiento, la pregunta es más general:

¿La incorporación de un factor permite explicar una parte relevante de la variabilidad observada en la respuesta?

Supóngase una variable respuesta \(Y\) y un factor con tres niveles. El modelo de ANOVA de un factor puede escribirse como:

\[ Y_{ij}=\mu+\tau_i+\varepsilon_{ij} \]

donde:

• \(Y_{ij}\) es la observación \(j\) del grupo \(i\);
• \(\mu\) es la media general;
• \(\tau_i\) es el efecto correspondiente al grupo \(i\);
• \(\varepsilon_{ij}\) es la variación que el modelo no explica.

La expresión anterior ya tiene la estructura:

\[ \text{Dato}=\text{Modelo}+\text{Error} \]

El modelo no contiene únicamente una prueba. Contiene una representación de cómo se supone que se originan los datos.

3.1 Modelo reducido y modelo completo

El modelo reducido sostiene que todos los individuos comparten la misma media:

\[ Y_{ij}=\mu+\varepsilon_{ij} \]

El modelo completo incorpora diferencias entre grupos:

\[ Y_{ij}=\mu+\tau_i+\varepsilon_{ij} \]

La prueba F del ANOVA evalúa si la reducción del error lograda por el modelo completo es suficientemente grande en relación con el error que todavía permanece.

\[ F= \frac{ \text{variabilidad explicada por el factor}/gl_{\text{factor}} }{ \text{variabilidad residual}/gl_{\text{error}} } \]

Distinción de efectos

En esta primera formulación se asumirá que el factor es fijo: los métodos A, B y C son exactamente los niveles que interesa comparar. La inferencia se dirige a esos niveles específicos.

Si los niveles hubieran sido seleccionados aleatoriamente de una población más amplia de niveles, el factor sería aleatorio y la pregunta ya no sería cuál nivel tiene mayor media, sino cuánta variabilidad existe entre niveles.

4 Construcción de un ejemplo reproducible

Se estudiará el desempeño final de estudiantes asignados a tres métodos de enseñanza. También se registró su calificación inicial y la escuela a la que pertenecen.

El diseño permitirá ilustrar tres componentes:

• Método de enseñanza: factor fijo.
• Calificación inicial: covariable cuantitativa.
• Escuela: factor aleatorio, porque las escuelas representan una muestra de una población mayor.

set.seed(2026)

n_escuelas <- 12
n_por_metodo_escuela <- 7

escuelas <- factor(
  paste0("E", sprintf("%02d", 1:n_escuelas))
)

datos <- expand.grid(
  escuela = escuelas,
  metodo = factor(
    c("A", "B", "C"),
    levels = c("A", "B", "C")
  ),
  estudiante = 1:n_por_metodo_escuela
)

efecto_escuela <- rnorm(
  n_escuelas,
  mean = 0,
  sd = 3.5
)

names(efecto_escuela) <- escuelas

media_pretest <- c(
  A = 60,
  B = 67,
  C = 74
)

efecto_metodo <- c(
  A = 0,
  B = 4,
  C = 7
)

datos$pretest <- rnorm(
  nrow(datos),
  mean = media_pretest[datos$metodo],
  sd = 8
)

datos$postest <- 20 +
  0.65 * datos$pretest +
  efecto_metodo[datos$metodo] +
  efecto_escuela[datos$escuela] +
  rnorm(
    nrow(datos),
    mean = 0,
    sd = 5
  )

datos$pretest_c <- datos$pretest - mean(datos$pretest)

head(datos)

4.1 Descripción del conjunto de datos simulado

El conjunto de datos representa un estudio educativo hipotético en el que se desea comparar el desempeño de estudiantes que reciben tres métodos de enseñanza diferentes: A, B y C. Los estudiantes pertenecen a distintas escuelas y cuentan con una medición académica antes y después de la aplicación del método de enseñanza.

La simulación considera que el resultado final de un estudiante puede estar relacionado con tres componentes principales:

1. Su nivel académico inicial.
2. El método de enseñanza recibido.
3. Las características particulares de la escuela a la que pertenece.

El objetivo será determinar si existen diferencias entre los métodos de enseñanza después de controlar estadísticamente el nivel inicial de los estudiantes. Además, se considerará que los estudiantes pertenecientes a una misma escuela pueden presentar resultados similares debido a que comparten docentes, instalaciones, recursos, contexto social y otras condiciones institucionales.

Por estas razones, el conjunto de datos permitirá desarrollar un análisis de covarianza con efectos mixtos, en el que el método de enseñanza y el pretest serán tratados como efectos fijos, mientras que la escuela será incorporada como un efecto aleatorio.

4.2 Estructura del estudio

Se simulan datos correspondientes a:

• 12 escuelas.
• 3 métodos de enseñanza.
• 7 estudiantes por cada combinación de escuela y método.

Por lo tanto, el número total de observaciones es:

[ 12 = 252]

Cada fila del conjunto de datos representa a un estudiante perteneciente a una escuela y asignado a uno de los tres métodos de enseñanza.

4.3 Variables del conjunto de datos

4.3.1 Escuela

La variable escuela identifica la institución educativa a la que pertenece cada estudiante. Las escuelas se representan mediante los códigos E01, E02, …, E12.

Esta variable es categórica y se almacena como un factor en R.

La escuela se utilizará como un efecto aleatorio, debido a que las doce escuelas pueden considerarse una muestra de una población más amplia de instituciones educativas. No se pretende comparar individualmente una escuela contra otra, sino reconocer que existen diferencias institucionales que pueden influir en el desempeño de los estudiantes.

Cada escuela recibe un efecto aleatorio generado mediante una distribución normal:

[ u_j N(0,3.5^2)]

donde (u_j) representa el efecto particular de la escuela (j).

Un valor positivo indica que la escuela tiende a presentar resultados superiores al promedio, mientras que un valor negativo indica resultados inferiores, manteniendo constantes el método y el nivel inicial de los estudiantes.

La inclusión de esta variable permite reconocer que las observaciones no son completamente independientes, ya que los estudiantes de una misma escuela comparten características institucionales.

4.3.2 Método de enseñanza

La variable metodo identifica el método educativo recibido por cada estudiante. Esta variable tiene tres categorías:

• Método A.
• Método B.
• Método C.

El método A se utiliza como categoría de referencia.

En la simulación se establecen los siguientes efectos:

[ _A = 0]

[ _B = 4]

[ _C = 7]

Esto significa que, después de controlar el pretest y el efecto de la escuela, se espera que:

• el método B produzca aproximadamente 4 puntos más que el método A;
• el método C produzca aproximadamente 7 puntos más que el método A.

El método se utilizará como un efecto fijo, debido a que el interés principal del análisis es comparar específicamente estos tres métodos de enseñanza.

La pregunta que se buscará responder es:

¿Existen diferencias estadísticamente significativas en el postest entre los métodos A, B y C, después de controlar el nivel inicial de los estudiantes?

4.3.3 Estudiante

La variable estudiante contiene los números del 1 al 7 para identificar a los estudiantes dentro de cada combinación de escuela y método.

Este número no constituye un identificador único para todo el conjunto de datos, debido a que se repite entre escuelas y métodos. La identificación completa de cada observación está determinada por la combinación de:

• escuela;
• método;
• número de estudiante.

Esta variable se utiliza principalmente para organizar la estructura del conjunto de datos. En el análisis propuesto no será incorporada directamente como predictor, debido a que cada estudiante cuenta con una única observación de pretest y postest.

4.3.4 Pretest

La variable pretest representa la calificación obtenida por cada estudiante antes de recibir el método de enseñanza.

Esta variable se utiliza como una medida del nivel académico inicial. En la simulación, los grupos presentan diferentes medias iniciales:

[ _{,A}=60]

[ _{,B}=67]

[ _{,C}=74]

Las puntuaciones individuales se generan alrededor de estas medias con una desviación estándar de 8 puntos.

Esto significa que los grupos no comienzan en condiciones equivalentes. Los estudiantes asignados al método C presentan, en promedio, un nivel inicial mayor que los estudiantes de los métodos A y B.

Este desequilibrio se introduce deliberadamente para mostrar una situación en la que un ANOVA aplicado únicamente al postest podría producir conclusiones engañosas. Una diferencia observada en el postest podría deberse parcialmente al método de enseñanza, pero también podría deberse a diferencias que ya existían antes de aplicar los métodos.

El pretest será utilizado como una covariable cuantitativa. Su función será ajustar la comparación entre los métodos de enseñanza según el nivel inicial de cada estudiante.

La pregunta asociada a esta variable es:

¿Cómo cambia el postest cuando aumenta la calificación inicial del estudiante?

En la simulación se establece una pendiente de 0.65, lo que significa que por cada punto adicional en el pretest se espera un aumento promedio de 0.65 puntos en el postest, manteniendo constantes el método y la escuela.

4.3.5 Postest

La variable postest representa la calificación obtenida por cada estudiante después de recibir el método de enseñanza.

Esta será la variable dependiente o variable respuesta del análisis.

El postest se genera mediante el siguiente modelo:

[ Y_{ijk} =======

20+ 0.65X_{ijk}+ k+ u_j+ {ijk} ]

donde:

• (Y_{ijk}) es el postest del estudiante;
• 20. es el intercepto general;
• (X_{ijk}) es el pretest;
• (0.65) es el efecto esperado del pretest;
• (_k) es el efecto fijo del método;
• (u_j) es el efecto aleatorio de la escuela;
• (_{ijk}) es el error residual individual.

El término residual se genera mediante:

[ _{ijk}N(0,5^2)]

Este componente representa la variabilidad individual que no puede explicarse mediante el pretest, el método de enseñanza o la escuela.

Puede incluir diferencias relacionadas con la motivación, asistencia, hábitos de estudio, condiciones familiares, errores de medición y otras características no observadas.

El propósito principal será explicar la variación del postest mediante el nivel inicial, el método de enseñanza y la escuela.

4.3.6 Pretest centrado

La variable pretest_c corresponde al pretest centrado alrededor de su media general. Se calcula mediante:

[ _i = _i-]

Cuando pretest_c es igual a cero, el estudiante tiene una puntuación inicial igual a la media general del conjunto de datos.

Los valores negativos indican que el estudiante se encuentra por debajo de la media inicial, mientras que los valores positivos indican que se encuentra por encima.

El centrado del pretest no modifica la relación entre el pretest y el postest, tampoco cambia las diferencias estimadas entre los métodos. Su principal beneficio consiste en mejorar la interpretación del intercepto.

Cuando se utiliza el pretest sin centrar, el intercepto representa el postest esperado para un estudiante con un pretest igual a cero. Esta situación puede carecer de sentido práctico si una puntuación de cero no se encuentra dentro del rango observado.

Al utilizar pretest_c, el intercepto representa el postest esperado para un estudiante con un nivel inicial promedio, perteneciente al método de referencia y a una escuela cuyo efecto aleatorio es igual a cero.

4.4 Propósito del análisis

Con este conjunto de datos se pretende mostrar la transición entre el ANOVA, el ANCOVA y los modelos mixtos.

En primer lugar, podría realizarse un ANOVA simple para comparar las medias del postest entre los métodos:

modelo_anova <- aov(postest ~ metodo, data = datos)
summary(modelo_anova)

##              Df Sum Sq Mean Sq F value              Pr(>F)    
## metodo        2   9934    4967    76.6 <0.0000000000000002 ***
## Residuals   249  16141      65                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sin embargo, este modelo no considera que los estudiantes comienzan con diferentes niveles de desempeño. Debido a que el método C contiene estudiantes con un pretest promedio mayor, una parte de su ventaja en el postest podría deberse a esa diferencia inicial y no exclusivamente al método utilizado.

Posteriormente, se ajustará un modelo ANCOVA:

modelo_ancova <- lm(
  postest ~ pretest_c + metodo,
  data = datos
)

summary(modelo_ancova)

## 
## Call:
## lm(formula = postest ~ pretest_c + metodo, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.666  -3.901   0.414   3.941  17.799 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  63.0092     0.7568   83.26 < 0.0000000000000002 ***
## pretest_c     0.6758     0.0498   13.57 < 0.0000000000000002 ***
## metodoB       2.6199     1.0105    2.59                 0.01 *  
## metodoC       5.7351     1.1808    4.86            0.0000021 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.11 on 248 degrees of freedom
## Multiple R-squared:  0.645,  Adjusted R-squared:  0.64 
## F-statistic:  150 on 3 and 248 DF,  p-value: <0.0000000000000002

Este modelo permitirá comparar los métodos después de controlar estadísticamente el pretest.

Finalmente, se incorporará la estructura jerárquica de los datos mediante un modelo mixto:

modelo_mixto <- lmer(
  postest ~ pretest_c + metodo + (1 | escuela),
  data = datos
)

summary(modelo_mixto)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: postest ~ pretest_c + metodo + (1 | escuela)
##    Data: datos
## 
## REML criterion at convergence: 1571
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -2.221 -0.650 -0.015  0.660  3.587 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  escuela  (Intercept) 10.8     3.29    
##  Residual             27.3     5.22    
## Number of obs: 252, groups:  escuela, 12
## 
## Fixed effects:
##             Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)  63.0307     1.1514  18.8258   54.74 < 0.0000000000000002 ***
## pretest_c     0.6787     0.0437 240.0435   15.51 < 0.0000000000000002 ***
## metodoB       2.5981     0.8666 237.4488    3.00                0.003 ** 
## metodoC       5.6924     1.0194 238.2003    5.58          0.000000064 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) prtst_ metodB
## pretest_c  0.273              
## metodoB   -0.426 -0.368       
## metodoC   -0.444 -0.612  0.593

En este modelo:

• pretest_c representa el efecto fijo del nivel inicial;
• metodo representa el efecto fijo del tratamiento educativo;
• (1 | escuela) representa un intercepto aleatorio para cada escuela;
• postest representa la variable dependiente.

4.5 Objetivos específicos

El análisis permitirá responder las siguientes preguntas:

1. ¿Existe una relación significativa entre el pretest y el postest?

2. ¿Existen diferencias entre los métodos de enseñanza después de controlar el nivel inicial de los estudiantes?

3. ¿Cuál método presenta el mayor postest ajustado?

4. ¿Cuánta variabilidad del postest se debe a diferencias entre escuelas?

5. ¿Los estudiantes pertenecientes a una misma escuela presentan resultados más similares entre sí que los estudiantes de escuelas diferentes?

6. ¿Cambian las conclusiones sobre los métodos cuando se incorpora el pretest como covariable?

7. ¿Cambian los errores estándar y las pruebas de significancia cuando se reconoce la agrupación de estudiantes dentro de las escuelas?

4.6 Interpretación conceptual del modelo

El conjunto de datos combina variables de diferentes niveles.

En el nivel individual se encuentran:

• el pretest;
• el postest;
• el método recibido;
• la variabilidad residual de cada estudiante.

En el nivel institucional se encuentra:

• la escuela.

Esta estructura puede representarse como:

[ ]

Por ello, el modelo mixto resulta más apropiado que un modelo lineal tradicional cuando se desea reconocer que los estudiantes de una misma institución comparten condiciones comunes.

El análisis no se limitará a preguntar cuál grupo presenta la media de postest más alta. Se buscará estimar cuál método produce mejores resultados cuando los estudiantes son comparados bajo un nivel inicial equivalente y cuando se controla la variabilidad atribuible a las escuelas.

4.7 Descripción inicial

datos |>
  group_by(metodo) |>
  summarise(
    n = n(),
    media_pretest = mean(pretest),
    sd_pretest = sd(pretest),
    media_postest = mean(postest),
    sd_postest = sd(postest),
    .groups = "drop"
  )

En los datos simulados, los grupos no inician con el mismo nivel promedio. El grupo C presenta un pretest mayor. Por tanto, una comparación exclusiva de las medias finales mezclará dos fuentes de explicación:

1. diferencias atribuibles al método;
2. diferencias que ya existían antes del método.

5 ANOVA como descomposición de la variabilidad

El siguiente gráfico muestra:

• observaciones individuales;
• media de cada método;
• media general.

media_general <- mean(datos$postest)

ggplot(
  datos,
  aes(
    x = metodo,
    y = postest
  )
) +
  geom_jitter(
    width = 0.12,
    alpha = 0.35,
    size = 1.7
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 4,
    shape = 18
  ) +
  geom_hline(
    yintercept = media_general,
    linetype = 2,
    linewidth = 0.9
  ) +
  labs(
    title = "ANOVA: variación entre grupos y dentro de los grupos",
    subtitle = "Los rombos representan medias de grupo y la línea discontinua la media general",
    x = "Método de enseñanza",
    y = "Calificación final"
  ) +
  theme_minimal(base_size = 13)

La variabilidad total puede descomponerse como:

\[ SS_{\text{total}} = SS_{\text{entre grupos}} + SS_{\text{dentro de grupos}} \]

En lenguaje de modelos:

\[ SS_{\text{total}} = SS_{\text{modelo}} + SS_{\text{error}} \]

5.1 ANOVA tradicional en R

modelo_aov <- aov(
  postest ~ metodo,
  data = datos
)

summary(modelo_aov)

##              Df Sum Sq Mean Sq F value              Pr(>F)    
## metodo        2   9934    4967    76.6 <0.0000000000000002 ***
## Residuals   249  16141      65                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.2 El mismo ANOVA mediante lm()

modelo_anova <- lm(
  postest ~ metodo,
  data = datos
)

anova(modelo_anova)

summary(modelo_anova)

## 
## Call:
## lm(formula = postest ~ metodo, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.948  -5.624   0.524   5.806  18.271 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)   58.153      0.878   66.20 < 0.0000000000000002 ***
## metodoB        7.546      1.242    6.07         0.0000000046 ***
## metodoC       15.378      1.242   12.38 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.05 on 249 degrees of freedom
## Multiple R-squared:  0.381,  Adjusted R-squared:  0.376 
## F-statistic: 76.6 on 2 and 249 DF,  p-value: <0.0000000000000002

Las funciones aov() y lm() producen el mismo ajuste porque el ANOVA de efectos fijos es un modelo lineal con un predictor categórico.

c(
  RSS_aov = deviance(modelo_aov),
  RSS_lm = deviance(modelo_anova)
)

## RSS_aov  RSS_lm 
##   16141   16141

all.equal(
  fitted(modelo_aov),
  fitted(modelo_anova)
)

## [1] TRUE

all.equal(
  residuals(modelo_aov),
  residuals(modelo_anova)
)

## [1] TRUE

Conclusión

ANOVA no es matemáticamente ajeno a la regresión. Es una regresión en la que los predictores representan categorías.

6 ¿Cómo entra un factor dentro de una regresión?

Una computadora no interpreta directamente etiquetas como A, B y C. El factor se representa mediante columnas numéricas en una matriz de diseño.

X_tratamiento <- model.matrix(
  ~ metodo,
  data = datos
)

head(X_tratamiento, 10)

##    (Intercept) metodoB metodoC
## 1            1       0       0
## 2            1       0       0
## 3            1       0       0
## 4            1       0       0
## 5            1       0       0
## 6            1       0       0
## 7            1       0       0
## 8            1       0       0
## 9            1       0       0
## 10           1       0       0

Con codificación de tratamiento, A funciona como categoría de referencia. El modelo es:

\[ Y_i = \beta_0 + \beta_1D_{Bi} + \beta_2D_{Ci} + \varepsilon_i \]

donde:

\[ D_{Bi}= \begin{cases} 1,& \text{si el individuo pertenece a B}\\ 0,& \text{en otro caso} \end{cases} \]

y análogamente para \(D_C\).

6.1 Interpretación de los coeficientes

coef(modelo_anova)

## (Intercept)     metodoB     metodoC 
##      58.153       7.546      15.378

La interpretación es:

• \(\beta_0\): media estimada del método A;
• \(\beta_1\): diferencia entre B y A;
• \(\beta_2\): diferencia entre C y A.

Se puede comprobar con las medias observadas.

medias_metodo <- tapply(
  datos$postest,
  datos$metodo,
  mean
)

medias_metodo

##     A     B     C 
## 58.15 65.70 73.53

c(
  media_A = coef(modelo_anova)[1],
  diferencia_B_A = medias_metodo["B"] - medias_metodo["A"],
  diferencia_C_A = medias_metodo["C"] - medias_metodo["A"]
)

## media_A.(Intercept)    diferencia_B_A.B    diferencia_C_A.C 
##              58.153               7.546              15.378

7 La ecuación común: el Modelo Lineal General

Regresión, ANOVA y ANCOVA pueden representarse mediante:

\[ \mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

donde:

• \(\mathbf{Y}\) es el vector de respuestas;
• \(\mathbf{X}\) es la matriz de diseño;
• \(\boldsymbol{\beta}\) contiene los parámetros;
• \(\boldsymbol{\varepsilon}\) contiene los errores.

Para el Modelo Lineal General clásico se supone:

\[ \boldsymbol{\varepsilon} \sim N(\mathbf{0},\sigma^2\mathbf{I}) \]

La estimación por mínimos cuadrados es:

\[ \widehat{\boldsymbol{\beta}} = (\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top}\mathbf{Y} \]

No es necesario realizar manualmente toda la operación para utilizar el modelo, pero sí comprender que las diferencias entre ANOVA, regresión y ANCOVA se encuentran principalmente en las columnas incluidas en \(\mathbf{X}\).

7.1 Qué cambia cuando se incorpora una covariable

La covariable puede no ser ortogonal al factor. Esto ocurre, por ejemplo, cuando las medias del pretest son distintas entre métodos.

metodo_h <- datos$metodo
contrasts(metodo_h) <- contr.helmert(3)

X_ancova <- model.matrix(
  ~ metodo_h + pretest_c,
  data = datos
)

round(
  crossprod(X_ancova),
  2
)

##             (Intercept) metodo_h1 metodo_h2 pretest_c
## (Intercept)         252       0.0         0       0.0
## metodo_h1             0     168.0         0     612.4
## metodo_h2             0       0.0       504    1785.0
## pretest_c             0     612.4      1785   23615.9

Los productos cruzados entre las columnas del factor y la covariable ya no son necesariamente cero.

Puente conceptual

En un ANOVA balanceado y ortogonal, cada efecto puede recibir una porción claramente separable de la variación.

En ANCOVA, el factor y la covariable pueden compartir información. Por eso se habla de un efecto del método ajustado por el pretest y se vuelve importante especificar qué modelo se está comparando.

8 Regresión: un modelo con predictor cuantitativo

La regresión simple sustituye el predictor categórico por una variable cuantitativa:

\[ Y_i = \beta_0 + \beta_1X_i + \varepsilon_i \]

modelo_regresion <- lm(
  postest ~ pretest_c,
  data = datos
)

summary(modelo_regresion)

## 
## Call:
## lm(formula = postest ~ pretest_c, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.892  -3.691  -0.229   3.988  18.392 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)  65.7942     0.4013   163.9 <0.0000000000000002 ***
## pretest_c     0.8212     0.0415    19.8 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.37 on 250 degrees of freedom
## Multiple R-squared:  0.611,  Adjusted R-squared:  0.609 
## F-statistic:  392 on 1 and 250 DF,  p-value: <0.0000000000000002

ggplot(
  datos,
  aes(
    x = pretest,
    y = postest
  )
) +
  geom_point(
    alpha = 0.4
  ) +
  geom_smooth(
    method = "lm",
    se = TRUE
  ) +
  labs(
    title = "Regresión lineal: un predictor cuantitativo",
    x = "Calificación inicial",
    y = "Calificación final"
  ) +
  theme_minimal(base_size = 13)

La pendiente indica cuánto cambia, en promedio, el postest por cada punto adicional en el pretest.

Una precisión importante

Que el pretest varíe entre estudiantes no lo convierte en un “efecto aleatorio” en el sentido de los modelos mixtos.

En esta ANCOVA, la pendiente del pretest es un efecto fijo: se estima un único coeficiente común para la población de interés.

Un efecto aleatorio implica que el coeficiente correspondiente se modela como una realización de una distribución, por ejemplo, una variación aleatoria entre escuelas.

9 ANCOVA: ANOVA y regresión en un solo modelo

ANCOVA incorpora ambos tipos de predictores:

\[ Y_{ij} = \beta_0 + \tau_i + \beta_1X_{ij} + \varepsilon_{ij} \]

En este ejemplo:

• metodo es categórico;
• pretest_c es cuantitativo;
• postest es la respuesta.

modelo_ancova <- lm(
  postest ~ metodo + pretest_c,
  data = datos
)

summary(modelo_ancova)

## 
## Call:
## lm(formula = postest ~ metodo + pretest_c, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.666  -3.901   0.414   3.941  17.799 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  63.0092     0.7568   83.26 < 0.0000000000000002 ***
## metodoB       2.6199     1.0105    2.59                 0.01 *  
## metodoC       5.7351     1.1808    4.86            0.0000021 ***
## pretest_c     0.6758     0.0498   13.57 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.11 on 248 degrees of freedom
## Multiple R-squared:  0.645,  Adjusted R-squared:  0.64 
## F-statistic:  150 on 3 and 248 DF,  p-value: <0.0000000000000002

La matriz de diseño contiene simultáneamente:

• una columna de intercepto;
• columnas que codifican el método;
• una columna para el pretest.

head(
  model.matrix(modelo_ancova),
  10
)

##    (Intercept) metodoB metodoC pretest_c
## 1            1       0       0    -9.270
## 2            1       0       0    -9.305
## 3            1       0       0   -27.873
## 4            1       0       0     3.278
## 5            1       0       0    -2.567
## 6            1       0       0    -5.758
## 7            1       0       0   -13.936
## 8            1       0       0    -1.980
## 9            1       0       0   -10.128
## 10           1       0       0    -8.816

9.1 Interpretación del modelo

Con la covariable centrada:

• el intercepto representa el postest esperado del método A para un estudiante con pretest igual a la media general;
• metodoB representa la diferencia ajustada entre B y A;
• metodoC representa la diferencia ajustada entre C y A;
• pretest_c representa el cambio esperado en el postest por cada punto adicional de pretest, manteniendo constante el método.

9.2 Comparación de modelos completos y reducidos

9.2.1 ¿El método agrega explicación después del pretest?

Modelo reducido:

\[ Y=\beta_0+\beta_1X+\varepsilon \]

Modelo completo:

\[ Y=\beta_0+\tau+\beta_1X+\varepsilon \]

anova(
  modelo_regresion,
  modelo_ancova
)

9.2.2 ¿El pretest agrega explicación después del método?

anova(
  modelo_anova,
  modelo_ancova
)

Esta lógica de comparar modelos completos y reducidos será central en el estudio formal del Modelo Lineal General.

10 Medias observadas y medias ajustadas

Las medias observadas resumen lo ocurrido en cada grupo sin igualar el pretest.

Las medias ajustadas responden:

¿Qué postest se esperaría en cada método si los grupos se compararan en el mismo valor de pretest?

emm_metodo <- emmeans(
  modelo_ancova,
  ~ metodo
)

emm_metodo

##  metodo emmean    SE  df lower.CL upper.CL
##  A        63.0 0.757 248     61.5     64.5
##  B        65.6 0.667 248     64.3     66.9
##  C        68.7 0.754 248     67.3     70.2
## 
## Confidence level used: 0.95

medias_crudas <- datos |>
  group_by(metodo) |>
  summarise(
    media = mean(postest),
    .groups = "drop"
  ) |>
  mutate(
    tipo = "Media observada"
  )

medias_ajustadas <- as.data.frame(
  emm_metodo
) |>
  transmute(
    metodo,
    media = emmean,
    tipo = "Media ajustada"
  )

bind_rows(
  medias_crudas,
  medias_ajustadas
) |>
  ggplot(
    aes(
      x = metodo,
      y = media,
      group = tipo,
      linetype = tipo,
      shape = tipo
    )
  ) +
  geom_point(
    size = 3
  ) +
  geom_line() +
  labs(
    title = "Medias observadas y medias ajustadas",
    subtitle = "El ajuste controla las diferencias iniciales en el pretest",
    x = "Método",
    y = "Calificación final estimada",
    linetype = NULL,
    shape = NULL
  ) +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = "bottom"
  )

10.1 Comparaciones múltiples ajustadas

pairs(
  emm_metodo,
  adjust = "tukey"
)

##  contrast estimate   SE  df t.ratio p.value
##  A - B       -2.62 1.01 248  -2.593  0.0272
##  A - C       -5.74 1.18 248  -4.857 <0.0001
##  B - C       -3.12 1.00 248  -3.100  0.0061
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Las comparaciones deben interpretarse sobre las medias ajustadas, no únicamente sobre las medias crudas.

11 Representación gráfica de ANCOVA

ggplot(
  datos,
  aes(
    x = pretest,
    y = postest,
    linetype = metodo,
    shape = metodo
  )
) +
  geom_point(
    alpha = 0.45
  ) +
  geom_smooth(
    method = "lm",
    formula = y ~ x,
    se = FALSE
  ) +
  labs(
    title = "ANCOVA con pendientes aproximadamente paralelas",
    subtitle = "El modelo aditivo supone una pendiente común para la covariable",
    x = "Calificación inicial",
    y = "Calificación final",
    linetype = "Método",
    shape = "Método"
  ) +
  theme_minimal(base_size = 13)

En un ANCOVA tradicional, las diferencias verticales entre las rectas representan diferencias ajustadas entre métodos.

12 Supuesto distintivo: homogeneidad de pendientes

El ANCOVA aditivo supone que el efecto del pretest es el mismo en todos los métodos.

Para evaluarlo se incorpora una interacción:

\[ Y = \beta_0 + \tau + \beta_1X + \tau X + \varepsilon \]

modelo_pendientes <- lm(
  postest ~ metodo * pretest_c,
  data = datos
)

anova(
  modelo_ancova,
  modelo_pendientes
)

Si la interacción no resulta significativa, el modelo con pendiente común suele ser una representación parsimoniosa.

Si la interacción es significativa, el efecto del método cambia con el nivel de la covariable.

No debe tratarse una interacción significativa como una simple molestia.

La heterogeneidad de pendientes es un resultado sustantivo: indica que las diferencias entre métodos dependen del nivel inicial del estudiante.

13 Ejemplo donde el ANCOVA tradicional no es suficiente

Se simulará un segundo conjunto de datos en el que cada método presenta una pendiente diferente.

set.seed(2035)

datos_heterogeneos <- data.frame(
  metodo = factor(
    rep(
      c("A", "B", "C"),
      each = 45
    ),
    levels = c("A", "B", "C")
  )
)

datos_heterogeneos$pretest_c <- rnorm(
  nrow(datos_heterogeneos),
  mean = 0,
  sd = 10
)

interceptos <- c(
  A = 64,
  B = 67,
  C = 70
)

pendientes <- c(
  A = 0.30,
  B = 0.65,
  C = 1.00
)

datos_heterogeneos$postest <- interceptos[datos_heterogeneos$metodo] +
  pendientes[datos_heterogeneos$metodo] *
  datos_heterogeneos$pretest_c +
  rnorm(
    nrow(datos_heterogeneos),
    mean = 0,
    sd = 4.5
  )

ggplot(
  datos_heterogeneos,
  aes(
    x = pretest_c,
    y = postest,
    linetype = metodo,
    shape = metodo
  )
) +
  geom_point(
    alpha = 0.45
  ) +
  geom_smooth(
    method = "lm",
    se = FALSE
  ) +
  labs(
    title = "Interacción entre método y covariable",
    subtitle = "Las diferencias entre métodos cambian según el pretest",
    x = "Pretest centrado",
    y = "Postest",
    linetype = "Método",
    shape = "Método"
  ) +
  theme_minimal(base_size = 13)

modelo_heterogeneo_aditivo <- lm(
  postest ~ metodo + pretest_c,
  data = datos_heterogeneos
)

modelo_heterogeneo_interaccion <- lm(
  postest ~ metodo * pretest_c,
  data = datos_heterogeneos
)

anova(
  modelo_heterogeneo_aditivo,
  modelo_heterogeneo_interaccion
)

summary(
  modelo_heterogeneo_interaccion
)

## 
## Call:
## lm(formula = postest ~ metodo * pretest_c, data = datos_heterogeneos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.287  -2.407   0.133   3.184  11.278 
## 
## Coefficients:
##                   Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)        64.5499     0.6833   94.47 < 0.0000000000000002 ***
## metodoB             2.3362     0.9708    2.41               0.0175 *  
## metodoC             5.5363     0.9715    5.70        0.00000007798 ***
## pretest_c           0.2218     0.0786    2.82               0.0055 ** 
## metodoB:pretest_c   0.4761     0.1000    4.76        0.00000508537 ***
## metodoC:pretest_c   0.7593     0.1080    7.03        0.00000000011 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.58 on 129 degrees of freedom
## Multiple R-squared:  0.737,  Adjusted R-squared:  0.727 
## F-statistic: 72.3 on 5 and 129 DF,  p-value: <0.0000000000000002

13.1 Pendientes por método

emtrends(
  modelo_heterogeneo_interaccion,
  ~ metodo,
  var = "pretest_c"
)

##  metodo pretest_c.trend     SE  df lower.CL upper.CL
##  A                0.222 0.0786 129   0.0662    0.377
##  B                0.698 0.0618 129   0.5757    0.820
##  C                0.981 0.0740 129   0.8347    1.128
## 
## Confidence level used: 0.95

13.2 Comparación de métodos en valores relevantes

emmeans(
  modelo_heterogeneo_interaccion,
  ~ metodo | pretest_c,
  at = list(
    pretest_c = c(-10, 0, 10)
  )
)

## pretest_c = -10:
##  metodo emmean    SE  df lower.CL upper.CL
##  A        62.3 1.050 129     60.2     64.4
##  B        59.9 0.987 129     58.0     61.9
##  C        60.3 1.080 129     58.1     62.4
## 
## pretest_c =   0:
##  metodo emmean    SE  df lower.CL upper.CL
##  A        64.5 0.683 129     63.2     65.9
##  B        66.9 0.690 129     65.5     68.2
##  C        70.1 0.691 129     68.7     71.5
## 
## pretest_c =  10:
##  metodo emmean    SE  df lower.CL upper.CL
##  A        66.8 1.030 129     64.7     68.8
##  B        73.9 0.861 129     72.2     75.6
##  C        79.9 0.935 129     78.0     81.8
## 
## Confidence level used: 0.95

Cuando existe interacción, no tiene sentido afirmar que un método es superior en todos los niveles del pretest sin examinar valores específicos.

14 Sumas de cuadrados y no ortogonalidad

En diseños no ortogonales, las sumas de cuadrados secuenciales dependen del orden de los términos.

anova(
  lm(
    postest ~ metodo + pretest_c,
    data = datos
  )
)

anova(
  lm(
    postest ~ pretest_c + metodo,
    data = datos
  )
)

14.1 Sumas de cuadrados tipo II

Para un modelo sin interacción, las sumas de cuadrados tipo II permiten evaluar cada efecto principal después de controlar los demás efectos principales.

tabla_tipo2 <- car::Anova(
  modelo_ancova,
  type = 2
)

tabla_tipo2

La selección del tipo de suma de cuadrados no debe hacerse como una regla mecánica. Debe responder a:

• la pregunta de investigación;
• la presencia de interacciones;
• el balance del diseño;
• la estructura de comparación entre modelos.

15 Tamaño del efecto ajustado

El eta cuadrado parcial para el método puede calcularse como:

\[ \eta_p^2 = \frac{ SS_{\text{método}} }{ SS_{\text{método}}+SS_{\text{error}} } \]

ss_metodo <- tabla_tipo2[
  "metodo",
  "Sum Sq"
]

ss_error <- sum(
  residuals(modelo_ancova)^2
)

eta_parcial_metodo <- ss_metodo /
  (ss_metodo + ss_error)

eta_parcial_metodo

## [1] 0.08723

El tamaño del efecto debe interpretarse junto con:

• intervalos de confianza;
• relevancia práctica;
• precisión de la medición;
• diseño del estudio;
• magnitud de las diferencias ajustadas. # Efectos fijos, aleatorios y mixtos

La distinción no depende únicamente de si una variable aparece como texto o como número. Depende de la pregunta inferencial y de cómo se conceptualizan sus niveles.

15.1 Efecto fijo

Un factor es fijo cuando interesa comparar específicamente los niveles observados.

Ejemplo:

• método A;
• método B;
• método C.

La pregunta es:

¿Cuáles métodos presentan diferencias?

Los parámetros se estiman directamente:

\[ \tau_A,\tau_B,\tau_C \]

15.2 Efecto aleatorio

Un factor es aleatorio cuando sus niveles representan una muestra de una población más amplia.

Ejemplo:

• 12 escuelas seleccionadas de una población de escuelas.

La pregunta es:

¿Cuánta variabilidad existe entre escuelas?

El efecto de escuela se modela como:

\[ b_j\sim N(0,\sigma_b^2) \]

No interesa necesariamente ordenar las escuelas. Interesa estimar \(\sigma_b^2\).

15.3 Modelo mixto

El modelo contiene simultáneamente:

• efectos fijos;
• efectos aleatorios.

Para este ejemplo:

\[ Y_{ijk} = \beta_0 + \tau_i + \beta_1X_{ijk} + b_j + \varepsilon_{ijk} \]

donde:

• método es fijo;
• pretest tiene pendiente fija;
• escuela tiene intercepto aleatorio.

modelo_mixto <- lmerTest::lmer(
  postest ~ metodo + pretest_c + (1 | escuela),
  data = datos,
  REML = TRUE
)

summary(
  modelo_mixto
)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: postest ~ metodo + pretest_c + (1 | escuela)
##    Data: datos
## 
## REML criterion at convergence: 1571
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -2.221 -0.650 -0.015  0.660  3.587 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  escuela  (Intercept) 10.8     3.29    
##  Residual             27.3     5.22    
## Number of obs: 252, groups:  escuela, 12
## 
## Fixed effects:
##             Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)  63.0307     1.1514  18.8258   54.74 < 0.0000000000000002 ***
## metodoB       2.5981     0.8666 237.4488    3.00                0.003 ** 
## metodoC       5.6924     1.0194 238.2003    5.58          0.000000064 ***
## pretest_c     0.6787     0.0437 240.0435   15.51 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) metodB metodC
## metodoB   -0.426              
## metodoC   -0.444  0.593       
## pretest_c  0.273 -0.368 -0.612

15.4 Componentes de la varianza

componentes <- as.data.frame(
  lme4::VarCorr(modelo_mixto)
)

componentes

var_escuela <- componentes$vcov[
  componentes$grp == "escuela"
]

var_residual <- componentes$vcov[
  componentes$grp == "Residual"
]

ICC <- var_escuela /
  (var_escuela + var_residual)

ICC

## [1] 0.2842

El coeficiente de correlación intraclase representa la proporción de la variabilidad total atribuible a diferencias entre escuelas.

Regla práctica

• Método específico de enseñanza: efecto fijo.
• Escuelas seleccionadas de una población: efecto aleatorio.
• Método fijo dentro de escuelas aleatorias: modelo mixto.
• Pretest: covariable con pendiente fija, salvo que se permita que su pendiente cambie aleatoriamente entre escuelas.

15.5 Pendiente aleatoria

Si se considera que la relación entre pretest y postest cambia entre escuelas, puede ajustarse:

modelo_pendiente_aleatoria <- lmerTest::lmer(
  postest ~ metodo + pretest_c +
    (1 + pretest_c | escuela),
  data = datos,
  REML = TRUE
)

Este modelo requiere suficiente información dentro de cada escuela. Una estructura aleatoria compleja no debe agregarse sin sustento teórico ni información suficiente.

16 ANOVA fijo, ANOVA aleatorio y modelo mixto

Situación	Interés principal	Parámetro relevante	Ejemplo
ANOVA fijo	Comparar niveles específicos	Medias o contrastes	Tres métodos seleccionados
ANOVA aleatorio	Generalizar a una población de niveles	Componente de varianza	Escuelas elegidas al azar
Modelo mixto	Comparar efectos específicos y modelar agrupación	Coeficientes y varianzas	Métodos dentro de escuelas

La elección del denominador de una prueba F depende de la estructura del modelo. En modelos aleatorios y mixtos no siempre corresponde utilizar automáticamente el cuadrado medio residual del ANOVA fijo.

17 Aplicación correcta de ANCOVA

Antes de utilizar ANCOVA deben responderse las siguientes preguntas.

17.1 ¿La covariable fue medida antes del tratamiento?

Una covariable posterior al tratamiento puede ser un mediador. Ajustarla puede eliminar una parte del efecto que se desea estudiar.

17.2 ¿Existe relación lineal con la respuesta?

Si la relación es curva, puede requerirse:

• término cuadrático;
• transformación;
• spline;
• modelo aditivo.

17.3 ¿Las pendientes son homogéneas?

Debe evaluarse la interacción:

postest ~ metodo * pretest

17.4 ¿Existe solapamiento?

ggplot(
  datos,
  aes(
    x = metodo,
    y = pretest
  )
) +
  geom_boxplot(
    alpha = 0.35
  ) +
  geom_jitter(
    width = 0.12,
    alpha = 0.3
  ) +
  labs(
    title = "Solapamiento de la covariable entre grupos",
    x = "Método",
    y = "Pretest"
  ) +
  theme_minimal(base_size = 13)

Si los grupos ocupan rangos completamente distintos de la covariable, las medias ajustadas dependerán de extrapolaciones.

17.5 ¿La asignación fue aleatoria?

ANCOVA puede mejorar precisión y controlar diferencias observadas, pero no convierte automáticamente un estudio observacional en un experimento aleatorizado.

17.6 ¿La covariable se mide con precisión?

El error de medición puede atenuar la pendiente y dejar confusión residual.

18 Guía de interpretación del ejemplo

Una conclusión académica debe separar los siguientes elementos.

18.1 Efecto de la covariable

Manteniendo constante el método, cada punto adicional en el pretest se asoció con un cambio estimado de \(\widehat{\beta}_1\) puntos en el postest.

18.2 Efecto fijo del método

Después de controlar el pretest, se evaluó si las medias ajustadas de los tres métodos eran iguales.

18.3 Efecto aleatorio de escuela

Se estimó la variabilidad atribuible a diferencias entre escuelas mediante un intercepto aleatorio.

18.4 Interacción

La interacción método por pretest permitió determinar si la pendiente de la covariable era común a los métodos.

19 Puente hacia la siguiente clase: el Modelo Lineal General

En esta sesión se han utilizado varias fórmulas:

19.0.1 ANOVA

\[ Y=\beta_0+\text{factor}+\varepsilon \]

19.0.2 Regresión

\[ Y=\beta_0+\beta_1X+\varepsilon \]

19.0.3 ANCOVA

\[ Y=\beta_0+\text{factor}+\beta_1X+\varepsilon \]

19.0.4 Modelo mixto

\[ Y=\mathbf{X}\boldsymbol{\beta} + \mathbf{Z}\mathbf{b} + \boldsymbol{\varepsilon} \]

Todos comparten una misma lógica:

1. especificar el modelo;
2. estimar parámetros;
3. comparar modelos;
4. revisar residuos;
5. interpretar predicciones y contrastes.

Base para la siguiente clase

El Modelo Lineal General no es una técnica adicional que aparece después de ANOVA y ANCOVA. Es el marco matemático que permite entenderlas conjuntamente.

20 Referencias conceptuales

Rutherford, A. (2011). ANOVA and ANCOVA: A GLM Approach (2.ª ed.). John Wiley & Sons.

Pérez López, C. (2007). Econometría básica: técnicas y herramientas. Pearson Educación.

Walpole, R. E., Myers, R. H., Myers, S. L. y Ye, K. (2012). Probabilidad y estadística para ingeniería y ciencias (9.ª ed.). Pearson Educación.

```