Capítulo 15: Regresión múltiple

Clase de un enfoque de Modelo Lineal Múltiple General, ANOVA y ANCOVA

1 Orientación de la clase

Esta clase desarrolla la regresión lineal múltiple como una ampliación de la regresión lineal simple y, simultáneamente, como una instancia del Modelo Lineal General. La idea central es que regresión, ANOVA y ANCOVA comparten una misma estructura:

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

La regresión utiliza predictores cuantitativos y también puede incorporar predictores cualitativos mediante variables indicadoras. Cuando el modelo contiene únicamente factores, puede representar un ANOVA. Cuando combina factores y covariables cuantitativas, puede representar un ANCOVA.

2 Del modelo simple al modelo múltiple

En regresión lineal simple se utiliza una sola variable independiente. En regresión múltiple se estudia una variable dependiente \(Y\) a partir de \(p\) variables independientes:

\[ Y_i = \beta_0+\beta_1X_{1i}+\beta_2X_{2i}+\cdots+\beta_pX_{pi}+\varepsilon_i. \]

La ecuación de regresión poblacional es:

\[ E(Y\mid X_1,\ldots,X_p) = \beta_0+\beta_1X_1+\cdots+\beta_pX_p. \]

La ecuación estimada mediante una muestra es:

\[ \widehat{Y}_i = b_0+b_1X_{1i}+\cdots+b_pX_{pi}. \]

El residual de la observación \(i\) es:

\[ e_i=Y_i-\widehat{Y}_i. \]

2.1 Interpretación condicional de los coeficientes

En un modelo múltiple, \(b_j\) representa el cambio esperado en \(Y\) ante un aumento de una unidad en \(X_j\), manteniendo constantes las demás variables incluidas en el modelo.

Esta última frase es indispensable. El coeficiente múltiple no describe necesariamente la relación marginal entre \(X_j\) y \(Y\). Describe una relación parcial o condicional, después de controlar estadísticamente los otros predictores.

Error frecuente: interpretar un coeficiente como causal. Un modelo de regresión puede estimar asociaciones ajustadas, pero la causalidad requiere un diseño, supuestos y una estrategia de identificación adecuados.

2.2 Lineal no significa necesariamente “recta simple”

La palabra lineal se refiere a que el modelo es lineal en los parámetros \(\beta\). Por ejemplo, los modelos siguientes siguen siendo modelos lineales:

\[ Y=\beta_0+\beta_1X+\beta_2X^2+\varepsilon, \]

\[ Y=\beta_0+\beta_1X+\beta_2Z+\beta_3XZ+\varepsilon. \]

El primero permite curvatura y el segundo permite interacción. En ambos casos, los parámetros aparecen de forma lineal.

3 Caso de estudio: costo anual de siniestros

Se construirá una cartera simulada de pólizas de automóvil. El objetivo es explicar y predecir el costo anual de siniestros a partir de características del vehículo, exposición y perfil del conductor.

Las variables serán:

• costo_siniestros: costo anual observado.
• valor_vehiculo_10k: valor del vehículo en unidades de 10 000.
• kilometraje_1000: kilometraje anual en kilometros.
• edad: edad del conductor.
• experiencia: años aproximados de experiencia de conducción.
• zona: zona de riesgo, baja, media o alta.
• tipo_vehiculo: sedán, SUV o pickup.

set.seed(1501)

n <- 220

zona <- factor(
  sample(
    c("Baja", "Media", "Alta"),
    size = n,
    replace = TRUE,
    prob = c(0.35, 0.40, 0.25)
  ),
  levels = c("Baja", "Media", "Alta")
)

tipo_vehiculo <- factor(
  sample(
    c("Sedan", "SUV", "Pickup"),
    size = n,
    replace = TRUE,
    prob = c(0.50, 0.30, 0.20)
  ),
  levels = c("Sedan", "SUV", "Pickup")
)

edad <- round(runif(n, min = 18, max = 70))

edad_licencia <- sample(16:22, size = n, replace = TRUE)
experiencia <- pmax(
  0,
  edad - edad_licencia + round(rnorm(n, mean = 0, sd = 2))
)

valor_vehiculo <- round(
  exp(rnorm(n, mean = log(90000), sd = 0.45)),
  digits = -2
)
valor_vehiculo <- pmin(pmax(valor_vehiculo, 25000), 300000)

kilometraje <- round(runif(n, min = 5000, max = 35000))

efecto_zona <- c(Baja = 0, Media = 450, Alta = 1000)[as.character(zona)]
efecto_tipo <- c(Sedan = 0, SUV = 300, Pickup = 550)[as.character(tipo_vehiculo)]

# La zona alta presenta además una pendiente mayor respecto del kilometraje.
efecto_interaccion <- ifelse(
  zona == "Alta",
  25 * ((kilometraje / 1000) - mean(kilometraje / 1000)),
  0
)

costo_siniestros <- (
  600 +
  220 * (valor_vehiculo / 10000) +
  28 * (kilometraje / 1000) -
  7 * edad +
  efecto_zona +
  efecto_tipo +
  efecto_interaccion +
  rnorm(n, mean = 0, sd = 700)
)

costo_siniestros <- round(pmax(costo_siniestros, 200), digits = 2)

polizas <- data.frame(
  costo_siniestros = costo_siniestros,
  valor_vehiculo_10k = valor_vehiculo / 10000,
  kilometraje_1000 = kilometraje / 1000,
  edad = edad,
  experiencia = experiencia,
  zona = zona,
  tipo_vehiculo = tipo_vehiculo
)

head(polizas)

3.1 Exploración inicial

summary(polizas)

##  costo_siniestros valor_vehiculo_10k kilometraje_1000      edad        experiencia       zona   
##  Min.   : 431.2   Min.   : 2.500     Min.   : 5.12    Min.   :18.00   Min.   : 0.00   Baja :84  
##  1st Qu.:2713.4   1st Qu.: 6.657     1st Qu.:12.55    1st Qu.:33.75   1st Qu.:15.00   Media:73  
##  Median :3401.0   Median : 8.875     Median :20.35    Median :47.00   Median :28.00   Alta :63  
##  Mean   :3625.7   Mean   : 9.893     Mean   :20.02    Mean   :45.66   Mean   :26.80             
##  3rd Qu.:4353.1   3rd Qu.:12.315     3rd Qu.:27.10    3rd Qu.:58.00   3rd Qu.:39.25             
##  Max.   :9825.6   Max.   :30.000     Max.   :34.94    Max.   :70.00   Max.   :56.00             
##  tipo_vehiculo
##  Sedan :116   
##  SUV   : 60   
##  Pickup: 44   
##               
##               
## 

variables_numericas <- polizas[
  c(
    "costo_siniestros",
    "valor_vehiculo_10k",
    "kilometraje_1000",
    "edad",
    "experiencia"
  )
]

round(cor(variables_numericas), 3)

##                    costo_siniestros valor_vehiculo_10k kilometraje_1000   edad experiencia
## costo_siniestros              1.000              0.762            0.211 -0.040      -0.044
## valor_vehiculo_10k            0.762              1.000            0.008 -0.052      -0.062
## kilometraje_1000              0.211              0.008            1.000  0.065       0.062
## edad                         -0.040             -0.052            0.065  1.000       0.982
## experiencia                  -0.044             -0.062            0.062  0.982       1.000

par(mfrow = c(1, 3))

plot(
  polizas$valor_vehiculo_10k,
  polizas$costo_siniestros,
  xlab = "Valor del vehículo, unidades de 10 000",
  ylab = "Costo anual de siniestros",
  main = "Costo y valor del vehículo",
  pch = 19
)

plot(
  polizas$kilometraje_1000,
  polizas$costo_siniestros,
  xlab = "Kilometraje anual, miles",
  ylab = "Costo anual de siniestros",
  main = "Costo y kilometraje",
  pch = 19
)

boxplot(
  costo_siniestros ~ zona,
  data = polizas,
  xlab = "Zona",
  ylab = "Costo anual de siniestros",
  main = "Costo por zona"
)

par(mfrow = c(1, 1))

Comprobación 1. Antes de ajustar un modelo, indique cuáles predictores parecen tener asociación con el costo. ¿La gráfica permite concluir que una asociación es causal?

Respuesta orientativa

El valor del vehículo, el kilometraje y la zona deberían mostrar asociaciones visibles. Sin embargo, una gráfica bivariada no controla otros factores y no demuestra causalidad. La asociación observada puede cambiar cuando se incorporan predictores correlacionados.

4 Estimación por mínimos cuadrados

El método de mínimos cuadrados ordinarios selecciona los coeficientes que minimizan:

\[ SCE=\sum_{i=1}^{n}(Y_i-\widehat{Y}_i)^2 =\sum_{i=1}^{n}e_i^2. \]

En notación matricial:

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

y, cuando \(\mathbf{X}'\mathbf{X}\) es invertible,

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

R realiza este cálculo con lm().

4.1 Primer modelo múltiple

modelo_1 <- lm(
  costo_siniestros ~ valor_vehiculo_10k + kilometraje_1000 + edad,
  data = polizas
)

summary(modelo_1)

## 
## Call:
## lm(formula = costo_siniestros ~ valor_vehiculo_10k + kilometraje_1000 + 
##     edad, data = polizas)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2673.59  -565.75   -19.45   567.69  2340.14 
## 
## Coefficients:
##                    Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)         845.557    259.009   3.265              0.00127 ** 
## valor_vehiculo_10k  220.464     12.151  18.143 < 0.0000000000000002 ***
## kilometraje_1000     32.795      6.668   4.918           0.00000173 ***
## edad                 -1.257      3.960  -0.317              0.75120    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 854.9 on 216 degrees of freedom
## Multiple R-squared:  0.6226, Adjusted R-squared:  0.6173 
## F-statistic: 118.8 on 3 and 216 DF,  p-value: < 0.00000000000000022

4.2 Tabla de coeficientes legible

tabla_coeficientes <- function(modelo, digitos = 4) {
  tabla <- as.data.frame(coef(summary(modelo)))
  tabla$termino <- rownames(tabla)
  rownames(tabla) <- NULL

  tabla <- tabla[
    c("termino", "Estimate", "Std. Error", "t value", "Pr(>|t|)")
  ]

  names(tabla) <- c(
    "Término",
    "Estimación",
    "Error estándar",
    "Estadístico t",
    "Valor p"
  )

  columnas_numericas <- 2:ncol(tabla)
  tabla[columnas_numericas] <- lapply(
    tabla[columnas_numericas],
    round,
    digits = digitos
  )

  tabla
}

knitr::kable(
  tabla_coeficientes(modelo_1),
  caption = "Coeficientes del modelo múltiple inicial"
)

Coeficientes del modelo múltiple inicial

Término	Estimación	Error estándar	Estadístico t	Valor p
(Intercept)	845.5570	259.0090	3.2646	0.0013
valor_vehiculo_10k	220.4641	12.1513	18.1433	0.0000
kilometraje_1000	32.7954	6.6685	4.9180	0.0000
edad	-1.2571	3.9597	-0.3175	0.7512

4.3 Interpretación sustantiva

Para interpretar un coeficiente:

1. Identifique la variable dependiente y sus unidades.
2. Identifique las unidades del predictor.
3. Indique el cambio esperado en la respuesta.
4. Añada “manteniendo constantes las demás variables”.
5. Distinga la estimación puntual de su incertidumbre.

Ejemplo de redacción:

Manteniendo constantes el kilometraje y la edad, un aumento de 10 000 unidades monetarias en el valor del vehículo se asocia con un cambio promedio de \(b_1\) unidades monetarias en el costo anual de siniestros.

El intercepto representa la respuesta esperada cuando todos los predictores cuantitativos valen cero y los factores se encuentran en sus categorías de referencia. Puede carecer de interpretación práctica cuando el valor cero está fuera del rango observado.

4.4 Centrado para mejorar la interpretación

Centrar una variable significa restarle una constante, usualmente su media:

\[ X_c=X-\overline{X}. \]

El centrado no cambia las predicciones ni el ajuste global. Sí cambia la interpretación del intercepto y de los efectos principales en modelos con interacción.

polizas$valor_c <- (
  polizas$valor_vehiculo_10k -
  mean(polizas$valor_vehiculo_10k)
)

polizas$km_c <- (
  polizas$kilometraje_1000 -
  mean(polizas$kilometraje_1000)
)

polizas$edad_c <- polizas$edad - mean(polizas$edad)

modelo_1_centrado <- lm(
  costo_siniestros ~ valor_c + km_c + edad_c,
  data = polizas
)

coef(modelo_1)

##        (Intercept) valor_vehiculo_10k   kilometraje_1000               edad 
##         845.557029         220.464138          32.795431          -1.257063

coef(modelo_1_centrado)

## (Intercept)     valor_c        km_c      edad_c 
## 3625.727455  220.464138   32.795431   -1.257063

all.equal(
  fitted(modelo_1),
  fitted(modelo_1_centrado),
  tolerance = 1e-10
)

## [1] TRUE

Comprobación 2. ¿Por qué cambió el intercepto después de centrar las variables, pero no cambiaron los valores ajustados?

Solución

El centrado cambia el punto de origen de los predictores. En el modelo centrado, el intercepto corresponde al costo esperado para una póliza con valores promedio en los predictores. La superficie ajustada es la misma, únicamente se reparametriza.

5 Descomposición de la variabilidad y bondad de ajuste

La variabilidad total puede descomponerse en:

\[ STC=SCR+SCE, \]

donde:

• \(STC=\sum(Y_i-\overline{Y})^2\), suma total de cuadrados.
• \(SCR=\sum(\widehat{Y}_i-\overline{Y})^2\), suma de cuadrados de la regresión.
• \(SCE=\sum(Y_i-\widehat{Y}_i)^2\), suma de cuadrados del error.

El coeficiente de determinación múltiple es:

\[ R^2=\frac{SCR}{STC}=1-\frac{SCE}{STC}. \]

Representa la proporción de variabilidad muestral de \(Y\) explicada por el modelo.

5.1 Cálculo manual de \(R^2\)

y <- polizas$costo_siniestros
y_ajustada <- fitted(modelo_1)
residuales <- residuals(modelo_1)

STC <- sum((y - mean(y))^2)
SCR <- sum((y_ajustada - mean(y))^2)
SCE <- sum(residuales^2)

c(STC = STC, SCR = SCR, SCE = SCE)

##       STC       SCR       SCE 
## 418216587 260362449 157854138

c(STC = STC, SCR_mas_SCE = SCR + SCE)

##         STC SCR_mas_SCE 
##   418216587   418216587

R2_manual <- SCR / STC
R2_R <- summary(modelo_1)$r.squared

c(R2_manual = R2_manual, R2_R = R2_R)

## R2_manual      R2_R 
## 0.6225541 0.6225541

5.2 \(R^2\) ajustado

Agregar predictores nunca reduce \(R^2\), incluso cuando el nuevo predictor aporta poco. Por ello se utiliza:

\[ R^2_{\text{ajustado}} = 1-(1-R^2)\frac{n-1}{n-p-1}. \]

Este indicador penaliza la inclusión de predictores adicionales.

summary(modelo_1)$adj.r.squared

## [1] 0.6173118

Un \(R^2\) alto no garantiza que el modelo sea correcto, causal, estable o útil fuera de la muestra. Un \(R^2\) bajo tampoco invalida automáticamente un efecto científicamente importante.

6 Inferencia: prueba F global y pruebas t individuales

6.1 Prueba F global

La prueba \(F\) global evalúa:

\[ H_0:\beta_1=\beta_2=\cdots=\beta_p=0, \]

contra:

\[ H_a:\text{al menos uno de los coeficientes de pendiente es distinto de cero}. \]

El estadístico es:

\[ F=\frac{CMR}{CME} = \frac{SCR/p}{SCE/(n-p-1)}. \]

anova(modelo_1)

summary(modelo_1)$fstatistic

##    value    numdf    dendf 
## 118.7558   3.0000 216.0000

6.2 Pruebas t individuales

Para cada predictor:

\[ H_0:\beta_j=0, \qquad H_a:\beta_j\neq 0, \]

con:

\[ t=\frac{b_j}{SE(b_j)}. \]

tabla_coeficientes(modelo_1)

confint(modelo_1, level = 0.95)

##                        2.5 %      97.5 %
## (Intercept)        335.04835 1356.065710
## valor_vehiculo_10k 196.51393  244.414349
## kilometraje_1000    19.65179   45.939069
## edad                -9.06174    6.547614

6.3 ¿Puede la prueba F ser significativa y algunos coeficientes no?

Sí. La prueba \(F\) evalúa el conjunto de predictores. Las pruebas \(t\) evalúan coeficientes individuales condicionados a los demás predictores. La discrepancia puede surgir por:

• predictores con efectos pequeños,
• multicolinealidad,
• tamaño de muestra insuficiente,
• variables redundantes,
• pruebas de hipótesis diferentes.

Comprobación 3. Explique por qué no debe concluirse que “el modelo no sirve” únicamente porque uno de sus coeficientes tenga un valor \(p\) mayor que 0.05.

Solución

El valor \(p\) individual responde si existe evidencia de que un coeficiente parcial difiere de cero, dado el resto del modelo. No evalúa por sí solo el desempeño predictivo, la relevancia sustantiva, la calidad del diseño, la precisión de otros coeficientes ni la significancia global.

7 Modelos completos y reducidos: conexión con ANOVA y ANCOVA

Una idea central del enfoque del Modelo Lineal General es evaluar un efecto comparando:

• un modelo completo, que contiene el efecto de interés;
• un modelo reducido, que omite ese efecto.

La pregunta es:

¿Cuánto se reduce el error al agregar el conjunto de términos?

Suponga que el modelo reducido tiene suma de cuadrados del error \(SCE_R\) y el completo \(SCE_C\). La reducción atribuible a los términos agregados es:

\[ SCE_R-SCE_C. \]

La prueba parcial \(F\) es:

\[ F= \frac{(SCE_R-SCE_C)/q} {SCE_C/(n-p_C-1)}, \]

donde \(q\) es el número de parámetros agregados.

7.1 6.1 Aporte conjunto de la zona

modelo_reducido <- lm(
  costo_siniestros ~ valor_c + km_c + edad_c,
  data = polizas
)

modelo_completo <- lm(
  costo_siniestros ~ valor_c + km_c + edad_c + zona,
  data = polizas
)

anova(modelo_reducido, modelo_completo)

anova(modelo_reducido, modelo_completo) realiza una prueba \(F\) parcial. En este ejemplo, la zona requiere dos coeficientes porque tiene tres niveles y uno funciona como categoría de referencia.

7.2 6.2 Cambio en \(R^2\)

R2_reducido <- summary(modelo_reducido)$r.squared
R2_completo <- summary(modelo_completo)$r.squared

delta_R2 <- R2_completo - R2_reducido

R2_parcial <- (
  deviance(modelo_reducido) - deviance(modelo_completo)
) / deviance(modelo_reducido)

c(
  R2_reducido = R2_reducido,
  R2_completo = R2_completo,
  cambio_R2 = delta_R2,
  R2_parcial_zona = R2_parcial
)

##     R2_reducido     R2_completo       cambio_R2 R2_parcial_zona 
##       0.6225541       0.7094446       0.0868905       0.2302065

El cambio en \(R^2\) mide cuánto aumenta la proporción total explicada. El \(R^2\) parcial mide la proporción del error del modelo reducido que desaparece al agregar el efecto.

7.3 Prueba de un término individual mediante modelos anidados

drop1(modelo_completo, test = "F")

drop1() elimina un término por vez y compara el modelo resultante con el modelo completo. Para un predictor numérico con un grado de libertad, la prueba parcial \(F\) es equivalente al cuadrado de la prueba \(t\):

\[ F=t^2. \]

t_edad <- coef(summary(modelo_completo))["edad_c", "t value"]

modelo_sin_edad <- update(modelo_completo, . ~ . - edad_c)
F_edad <- anova(modelo_sin_edad, modelo_completo)$F[2]

c(t_cuadrado = t_edad^2, F_parcial = F_edad)

## t_cuadrado  F_parcial 
##   1.315186   1.315186

8 Sumas de cuadrados y orden de entrada

En diseños no ortogonales o con predictores correlacionados, la variabilidad explicada por dos predictores puede superponerse. Por ello es necesario indicar qué hipótesis se está probando.

8.1 Sumas de cuadrados secuenciales

La tabla producida por anova(modelo) en R asigna la variabilidad en el orden de entrada de los términos. Estas son sumas de cuadrados secuenciales o Tipo I.

modelo_orden_1 <- lm(
  costo_siniestros ~ edad + experiencia + zona,
  data = polizas
)

modelo_orden_2 <- lm(
  costo_siniestros ~ experiencia + edad + zona,
  data = polizas
)

anova(modelo_orden_1)

anova(modelo_orden_2)

Como edad y experiencia están correlacionadas, el predictor que entra primero puede recibir una parte importante de la variabilidad compartida.

8.2 Pruebas parciales

drop1(modelo_orden_1, test = "F")

Estas pruebas preguntan qué aporta cada término después de controlar por todos los demás términos incluidos.

Regla conceptual: no elija un tipo de suma de cuadrados únicamente porque produzca un valor \(p\) conveniente. Defina primero la hipótesis científica, la jerarquía del modelo y el orden lógico de ajuste.

9 Variables cualitativas y codificación

La regresión puede incorporar variables cualitativas. Si un factor tiene \(k\) categorías, se requieren normalmente \(k-1\) columnas para evitar sobreparametrización.

9.1 Codificación de referencia

R utiliza por defecto contrastes de tratamiento. La primera categoría del factor es la referencia.

contrasts(polizas$zona)

##       Media Alta
## Baja      0    0
## Media     1    0
## Alta      0    1

head(
  model.matrix(
    ~ zona,
    data = polizas
  )
)

##   (Intercept) zonaMedia zonaAlta
## 1           1         1        0
## 2           1         1        0
## 3           1         0        1
## 4           1         1        0
## 5           1         0        0
## 6           1         1        0

Para zona con referencia "Baja":

• el intercepto representa la media o respuesta ajustada de la zona baja;
• zonaMedia representa la diferencia entre media y baja;
• zonaAlta representa la diferencia entre alta y baja.

9.2 Modelo con predictores cualitativos

modelo_2 <- lm(
  costo_siniestros ~
    valor_c +
    km_c +
    edad_c +
    zona +
    tipo_vehiculo,
  data = polizas
)

knitr::kable(
  tabla_coeficientes(modelo_2),
  caption = "Modelo con predictores cuantitativos y cualitativos"
)

Modelo con predictores cuantitativos y cualitativos

Término	Estimación	Error estándar	Estadístico t	Valor p
(Intercept)	3032.2671	94.0683	32.2348	0.0000
valor_c	216.0369	10.3912	20.7903	0.0000
km_c	34.2100	5.6786	6.0244	0.0000
edad_c	-6.0754	3.4027	-1.7854	0.0756
zonaMedia	329.1864	117.0294	2.8129	0.0054
zonaAlta	973.0689	121.7584	7.9918	0.0000
tipo_vehiculoSUV	366.3744	116.4605	3.1459	0.0019
tipo_vehiculoPickup	528.2923	131.4540	4.0188	0.0001

9.3 Cambiar la categoría de referencia

polizas$zona_ref_alta <- relevel(polizas$zona, ref = "Alta")

modelo_ref_alta <- lm(
  costo_siniestros ~ valor_c + km_c + edad_c + zona_ref_alta,
  data = polizas
)

coef(modelo_completo)

## (Intercept)     valor_c        km_c      edad_c   zonaMedia    zonaAlta 
## 3252.864068  212.848576   33.111506   -4.026183  264.631107  995.426575

coef(modelo_ref_alta)

##        (Intercept)            valor_c               km_c             edad_c  zona_ref_altaBaja 
##        4248.290643         212.848576          33.111506          -4.026183        -995.426575 
## zona_ref_altaMedia 
##        -730.795468

all.equal(
  fitted(modelo_completo),
  fitted(modelo_ref_alta),
  tolerance = 1e-10
)

## [1] TRUE

Cambiar la referencia modifica la parametrización y las comparaciones individuales, pero no cambia los valores ajustados, los residuales, \(R^2\) ni la prueba global del factor.

9.4 Codificación de efectos

Con contrastes de suma, los coeficientes comparan niveles respecto de una media general parametrizada.

polizas$zona_suma <- polizas$zona
contrasts(polizas$zona_suma) <- contr.sum(3)

modelo_suma <- lm(
  costo_siniestros ~ valor_c + km_c + edad_c + zona_suma,
  data = polizas
)

model.matrix(modelo_suma)[1:8, ]

##   (Intercept)   valor_c       km_c     edad_c zona_suma1 zona_suma2
## 1           1 -2.672636   7.114064  22.340909          0          1
## 2           1 -1.312636  13.917064 -24.659091          0          1
## 3           1 20.107364  -7.214936 -20.659091         -1         -1
## 4           1 -3.892636   7.070064 -12.659091          0          1
## 5           1 -2.162636   2.150064   6.340909          1          0
## 6           1  1.157364 -14.900936   5.340909          0          1
## 7           1 -2.122636   3.066064  15.340909          0          1
## 8           1  2.747364  -2.600936 -15.659091          1          0

coef(modelo_suma)

## (Intercept)     valor_c        km_c      edad_c  zona_suma1  zona_suma2 
## 3672.883295  212.848576   33.111506   -4.026183 -420.019228 -155.388120

all.equal(
  fitted(modelo_completo),
  fitted(modelo_suma),
  tolerance = 1e-10
)

## [1] TRUE

10 ANOVA como regresión

Un ANOVA de un factor puede escribirse como una regresión que utiliza únicamente un predictor cualitativo:

\[ Y_i=\beta_0+\beta_1D_{1i}+\cdots+\beta_{k-1}D_{k-1,i}+\varepsilon_i. \]

10.1 Comparación de medias por zona

modelo_anova <- lm(
  costo_siniestros ~ zona,
  data = polizas
)

anova(modelo_anova)

summary(modelo_anova)

## 
## Call:
## lm(formula = costo_siniestros ~ zona, data = polizas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3945.3  -770.2  -111.0   718.8  5449.2 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)   3199.6      141.6  22.591 < 0.0000000000000002 ***
## zonaMedia      268.5      207.7   1.293                0.197    
## zonaAlta      1176.8      216.3   5.439          0.000000144 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1298 on 217 degrees of freedom
## Multiple R-squared:  0.1257, Adjusted R-squared:  0.1176 
## F-statistic:  15.6 on 2 and 217 DF,  p-value: 0.0000004682

aggregate(
  costo_siniestros ~ zona,
  data = polizas,
  FUN = mean
)

La prueba \(F\) de anova(modelo_anova) evalúa la igualdad de las medias poblacionales de las zonas.

10.2 Equivalencia entre prueba global y comparación de modelos

modelo_intercepto <- lm(
  costo_siniestros ~ 1,
  data = polizas
)

anova(modelo_intercepto, modelo_anova)

El modelo reducido contiene solo una media general. El modelo completo permite medias diferentes por zona. La reducción del error es la variabilidad entre grupos.

11 ANCOVA como regresión múltiple

El ANCOVA combina:

• factores, que representan grupos o condiciones;
• covariables cuantitativas, que explican variabilidad adicional.

Un modelo básico es:

\[ Y_i = \beta_0+ \beta_1D_{1i}+\cdots+ \beta_{k-1}D_{k-1,i}+ \gamma X_i+ \varepsilon_i. \]

El efecto del factor se evalúa después de controlar la covariable.

11.1 Modelo ANCOVA para las zonas

modelo_ancova <- lm(
  costo_siniestros ~ zona + valor_c + km_c + edad_c,
  data = polizas
)

summary(modelo_ancova)

## 
## Call:
## lm(formula = costo_siniestros ~ zona + valor_c + km_c + edad_c, 
##     data = polizas)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2271.94  -525.81     3.67   522.46  1860.62 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) 3252.864     82.535  39.412 < 0.0000000000000002 ***
## zonaMedia    264.631    120.931   2.188               0.0297 *  
## zonaAlta     995.427    126.661   7.859    0.000000000000187 ***
## valor_c      212.849     10.757  19.788 < 0.0000000000000002 ***
## km_c          33.112      5.881   5.630    0.000000056252382 ***
## edad_c        -4.026      3.511  -1.147               0.2527    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 753.5 on 214 degrees of freedom
## Multiple R-squared:  0.7094, Adjusted R-squared:  0.7027 
## F-statistic: 104.5 on 5 and 214 DF,  p-value: < 0.00000000000000022

drop1(modelo_ancova, test = "F")

La prueba del factor zona compara respuestas ajustadas a valores comunes de las covariables.

11.2 Medias ajustadas mediante predicción

Como las covariables fueron centradas, el valor cero representa su media. Se pueden obtener respuestas ajustadas por zona:

escenario_zonas <- data.frame(
  zona = factor(
    levels(polizas$zona),
    levels = levels(polizas$zona)
  ),
  valor_c = 0,
  km_c = 0,
  edad_c = 0
)

medias_ajustadas <- predict(
  modelo_ancova,
  newdata = escenario_zonas,
  interval = "confidence",
  level = 0.95
)

resultado_ajustado <- cbind(
  escenario_zonas["zona"],
  as.data.frame(medias_ajustadas)
)

knitr::kable(
  resultado_ajustado,
  digits = 2,
  caption = "Costos medios ajustados por zona"
)

Costos medios ajustados por zona

zona	fit	lwr	upr
Baja	3252.86	3090.18	3415.55
Media	3517.50	3343.44	3691.55
Alta	4248.29	4060.17	4436.41

Las medias ajustadas son predicciones del modelo para valores específicos de las covariables. Debe indicarse siempre en qué valores se realizó el ajuste.

12 Homogeneidad de pendientes e interacción

El ANCOVA tradicional supone que la pendiente de la covariable es común en todos los grupos. Este supuesto puede expresarse como ausencia de interacción entre el factor y la covariable.

Modelo de pendientes homogéneas:

\[ Y=\beta_0+\beta_1D+\beta_2X+\varepsilon. \]

Modelo de pendientes heterogéneas:

\[ Y=\beta_0+\beta_1D+\beta_2X+\beta_3DX+\varepsilon. \]

12.1 Evaluación de la interacción zona por kilometraje

modelo_homogeneo <- lm(
  costo_siniestros ~
    zona +
    valor_c +
    km_c +
    edad_c,
  data = polizas
)

modelo_heterogeneo <- lm(
  costo_siniestros ~
    zona * km_c +
    valor_c +
    edad_c,
  data = polizas
)

anova(modelo_homogeneo, modelo_heterogeneo)

Si la interacción es relevante, la diferencia entre zonas depende del kilometraje. En ese caso no existe una única diferencia de zona válida para todos los niveles de la covariable.

12.2 Interpretación jerárquica

En un modelo con zona * km_c, R incluye:

• zona,
• km_c,
• zona:km_c.

Los efectos principales se interpretan cuando la otra variable vale cero. Debido al centrado, cero corresponde al kilometraje promedio.

knitr::kable(
  tabla_coeficientes(modelo_heterogeneo),
  caption = "Modelo con pendientes diferentes por zona"
)

Modelo con pendientes diferentes por zona

Término	Estimación	Error estándar	Estadístico t	Valor p
(Intercept)	3252.6906	82.5858	39.3856	0.0000
zonaMedia	267.8050	120.9811	2.2136	0.0279
zonaAlta	996.5458	126.7289	7.8636	0.0000
km_c	30.6309	9.1391	3.3516	0.0010
valor_c	213.8184	10.7829	19.8294	0.0000
edad_c	-3.8724	3.5189	-1.1004	0.2724
zonaMedia:km_c	-4.0299	13.4448	-0.2997	0.7647
zonaAlta:km_c	16.9471	15.2905	1.1083	0.2690

12.3 Predicciones por zona y kilometraje

rejilla <- expand.grid(
  km_c = seq(
    min(polizas$km_c),
    max(polizas$km_c),
    length.out = 80
  ),
  zona = factor(
    levels(polizas$zona),
    levels = levels(polizas$zona)
  ),
  valor_c = 0,
  edad_c = 0
)

rejilla$prediccion <- predict(
  modelo_heterogeneo,
  newdata = rejilla
)

limites_y <- range(
  polizas$costo_siniestros,
  rejilla$prediccion
)

plot(
  NA,
  xlim = range(polizas$km_c),
  ylim = limites_y,
  xlab = "Kilometraje centrado, miles",
  ylab = "Costo anual de siniestros",
  main = "Pendientes estimadas por zona"
)

for (nivel in levels(polizas$zona)) {
  datos_nivel <- rejilla[rejilla$zona == nivel, ]
  lines(
    datos_nivel$km_c,
    datos_nivel$prediccion,
    lwd = 2,
    lty = which(levels(polizas$zona) == nivel)
  )
}

legend(
  "topleft",
  legend = levels(polizas$zona),
  lty = seq_along(levels(polizas$zona)),
  lwd = 2,
  bty = "n"
)

Comprobación 4. En el modelo con interacción, ¿qué representa el coeficiente de zonaAlta:km_c?

Solución

Representa cuánto difiere la pendiente del kilometraje en la zona alta respecto de la pendiente de la zona de referencia, que es la zona baja. La pendiente para la zona alta se obtiene sumando el coeficiente de km_c y el coeficiente de zonaAlta:km_c.

13 Multicolinealidad

La multicolinealidad aparece cuando un predictor puede explicarse de forma importante mediante otros predictores.

Consecuencias posibles:

• errores estándar elevados;
• coeficientes sensibles a pequeñas modificaciones;
• signos inesperados;
• pruebas individuales no significativas, aunque la prueba global sí lo sea;
• dificultad para separar contribuciones individuales.

La multicolinealidad no implica necesariamente malas predicciones dentro del rango observado, pero sí limita la interpretación individual de los coeficientes.

13.1 Edad y experiencia

cor(polizas$edad, polizas$experiencia)

## [1] 0.9820489

modelo_colineal <- lm(
  costo_siniestros ~
    valor_vehiculo_10k +
    kilometraje_1000 +
    edad +
    experiencia,
  data = polizas
)

summary(modelo_colineal)

## 
## Call:
## lm(formula = costo_siniestros ~ valor_vehiculo_10k + kilometraje_1000 + 
##     edad + experiencia, data = polizas)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2665.44  -569.13     6.02   560.03  2355.65 
## 
## Coefficients:
##                    Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)         986.504    455.092   2.168               0.0313 *  
## valor_vehiculo_10k  220.726     12.195  18.099 < 0.0000000000000002 ***
## kilometraje_1000     32.809      6.682   4.910            0.0000018 ***
## edad                 -9.027     20.986  -0.430               0.6675    
## experiencia           7.872     20.880   0.377               0.7065    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 856.6 on 215 degrees of freedom
## Multiple R-squared:  0.6228, Adjusted R-squared:  0.6158 
## F-statistic: 88.75 on 4 and 215 DF,  p-value: < 0.00000000000000022

13.2 Cálculo manual del VIF

Para el predictor \(X_j\):

\[ VIF_j=\frac{1}{1-R_j^2}, \]

donde \(R_j^2\) proviene de predecir \(X_j\) con los demás predictores.

vif_manual <- function(modelo) {
  X <- model.matrix(modelo)

  # Se elimina el intercepto, si existe.
  if ("(Intercept)" %in% colnames(X)) {
    X <- X[, colnames(X) != "(Intercept)", drop = FALSE]
  }

  if (ncol(X) < 2) {
    stop("Se necesitan al menos dos columnas predictoras.")
  }

  resultado <- numeric(ncol(X))
  names(resultado) <- colnames(X)

  for (j in seq_len(ncol(X))) {
    respuesta <- X[, j]
    predictores <- X[, -j, drop = FALSE]

    modelo_auxiliar <- lm(respuesta ~ predictores)
    R2_j <- summary(modelo_auxiliar)$r.squared

    resultado[j] <- 1 / (1 - R2_j)
  }

  resultado
}

round(vif_manual(modelo_colineal), 3)

## valor_vehiculo_10k   kilometraje_1000               edad        experiencia 
##              1.006              1.004             28.174             28.198

Reglas como VIF mayor que 5 o 10 son señales orientativas, no leyes universales. Debe examinarse el contexto, la precisión requerida y la función del predictor.

13.3 Comparar estabilidad de coeficientes

modelo_sin_experiencia <- update(
  modelo_colineal,
  . ~ . - experiencia
)

coef(modelo_colineal)

##        (Intercept) valor_vehiculo_10k   kilometraje_1000               edad        experiencia 
##         986.503868         220.726036          32.809378          -9.026702           7.872147

coef(modelo_sin_experiencia)

##        (Intercept) valor_vehiculo_10k   kilometraje_1000               edad 
##         845.557029         220.464138          32.795431          -1.257063

Posibles decisiones: conservar ambos predictores por razones teóricas, eliminar uno redundante, combinarlos, redefinir el objetivo inferencial, aumentar información muestral o utilizar regularización cuando el objetivo principal sea predictivo.

14 Supuestos del modelo

Para inferencia clásica por mínimos cuadrados se revisan los siguientes supuestos.

14.1 Linealidad de la media

\[ E(Y\mid X)=X\beta. \]

La forma funcional debe representar razonablemente la media condicional. La falta de linealidad puede requerir transformaciones, términos polinómicos o interacciones.

14.2 Media condicional cero

\[ E(\varepsilon\mid X)=0. \]

Este supuesto implica que el error no contiene un patrón sistemático relacionado con los predictores. Puede violarse por variables omitidas, simultaneidad, errores de medición o selección de muestra.

14.3 Varianza constante

\[ Var(\varepsilon_i\mid X)=\sigma^2. \]

La heterocedasticidad altera los errores estándar clásicos, aunque los coeficientes OLS pueden continuar siendo insesgados bajo otros supuestos.

14.4 Independencia

Los errores no deben presentar dependencia no modelada. La independencia depende principalmente del diseño de recolección. Datos longitudinales, agrupados o temporales requieren estructuras específicas.

14.5 Normalidad del error

La normalidad es importante para inferencia exacta en muestras pequeñas. No es necesaria para que mínimos cuadrados defina la recta o superficie que minimiza la suma de cuadrados.

14.6 Ausencia de multicolinealidad perfecta

Ninguna columna de la matriz de diseño debe ser una combinación lineal exacta de otras columnas.

Los diagnósticos no “demuestran” que los supuestos sean verdaderos. Permiten identificar patrones incompatibles con ellos y evaluar la sensibilidad de las conclusiones.

15 Diagnóstico residual

Se utilizará el modelo con factores como modelo principal.

modelo_diagnostico <- modelo_2

15.1 Gráficos estándar de R

par(mfrow = c(2, 2))
plot(modelo_diagnostico)

par(mfrow = c(1, 1))

15.2 Prueba manual de Breusch-Pagan

La prueba de Breusch-Pagan puede implementarse mediante una regresión auxiliar de los residuales cuadrados.

breusch_pagan_manual <- function(modelo) {
  e2 <- residuals(modelo)^2
  X <- model.matrix(modelo)

  modelo_auxiliar <- lm(e2 ~ X[, -1, drop = FALSE])
  LM <- length(e2) * summary(modelo_auxiliar)$r.squared
  gl <- ncol(X) - 1
  valor_p <- pchisq(LM, df = gl, lower.tail = FALSE)

  data.frame(
    estadistico_LM = LM,
    grados_libertad = gl,
    valor_p = valor_p
  )
}

breusch_pagan_manual(modelo_diagnostico)

No debe utilizarse una prueba de supuestos como un interruptor mecánico de “válido/no válido”. Examine también magnitud, gráficos y consecuencias prácticas.

16 Observaciones atípicas, leverage e influencia

16.1 Residuales estudentizados eliminados

Los residuales estudentizados eliminados recalculan la estimación de la variabilidad dejando fuera temporalmente la observación evaluada.

rstudent_valores <- rstudent(modelo_diagnostico)

summary(rstudent_valores)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.68952 -0.66002 -0.01775 -0.00022  0.55182  2.61214

which(abs(rstudent_valores) > 2)

##  25  39  76  79 107 111 120 161 172 200 205 215 
##  25  39  76  79 107 111 120 161 172 200 205 215

El umbral de 2 es una regla inicial. Cuando se examinan muchas observaciones, conviene considerar ajustes por multiplicidad y revisar el contexto.

16.2 Leverage

El leverage \(h_{ii}\) mide qué tan inusual es la combinación de valores predictores de una observación.

h <- hatvalues(modelo_diagnostico)

p_parametros <- length(coef(modelo_diagnostico))
n_observaciones <- nobs(modelo_diagnostico)

umbral_h_2 <- 2 * p_parametros / n_observaciones
umbral_h_3 <- 3 * p_parametros / n_observaciones

c(
  umbral_dos = umbral_h_2,
  umbral_tres = umbral_h_3
)

##  umbral_dos umbral_tres 
##  0.07272727  0.10909091

which(h > umbral_h_3)

## 201 
## 201

16.3 Distancia de Cook

La distancia de Cook resume cuánto podrían cambiar los valores ajustados al eliminar una observación.

D <- cooks.distance(modelo_diagnostico)
umbral_cook <- 4 / n_observaciones

which(D > umbral_cook)

##  21  25  79 107 109 120 161 172 200 201 215 
##  21  25  79 107 109 120 161 172 200 201 215

plot(
  D,
  type = "h",
  xlab = "Observación",
  ylab = "Distancia de Cook",
  main = "Influencia de las observaciones"
)
abline(h = umbral_cook, lty = 2)

16.4 Tabla consolidada de diagnóstico

diagnostico_observaciones <- data.frame(
  observacion = seq_len(n_observaciones),
  residual = residuals(modelo_diagnostico),
  residual_estudentizado = rstudent_valores,
  leverage = h,
  cook = D
)

diagnostico_observaciones <- diagnostico_observaciones[
  order(diagnostico_observaciones$cook, decreasing = TRUE),
]

knitr::kable(
  head(diagnostico_observaciones, 10),
  digits = 4,
  caption = "Diez observaciones con mayor distancia de Cook"
)

Diez observaciones con mayor distancia de Cook

	observacion	residual	residual_estudentizado	leverage	cook
215	215	-2517.5822	-3.6895	0.0580	0.0989
107	107	-2465.4706	-3.5744	0.0411	0.0648
200	200	-1912.7390	-2.7460	0.0450	0.0431
21	21	1370.6822	1.9685	0.0617	0.0314
161	161	1678.9582	2.3974	0.0425	0.0312
120	120	1693.7403	2.4139	0.0385	0.0285
172	172	1834.3373	2.6121	0.0325	0.0279
201	201	855.8091	1.2570	0.1125	0.0250
79	79	-1428.9433	-2.0355	0.0450	0.0240
25	25	-1415.7209	-2.0127	0.0417	0.0217

Una observación influyente no debe eliminarse automáticamente. Primero verifique errores de registro, pertenencia a la población, mecanismo de generación y sensibilidad de las conclusiones con y sin la observación.

17 Estimación y predicción

Para una combinación de predictores \(x_0\), pueden construirse dos intervalos diferentes.

17.1 Intervalo de confianza para la respuesta media

Estima:

\[ E(Y\mid X=x_0). \]

17.2 Intervalo de predicción para una observación individual

Estima un nuevo valor individual:

\[ Y_{\text{nuevo}}\mid X=x_0. \]

El intervalo de predicción es más amplio porque incluye:

• incertidumbre en la media estimada;
• variabilidad individual alrededor de la media.

17.3 Ejemplo

nueva_poliza <- data.frame(
  valor_c = 12 - mean(polizas$valor_vehiculo_10k),
  km_c = 18 - mean(polizas$kilometraje_1000),
  edad_c = 35 - mean(polizas$edad),
  zona = factor("Media", levels = levels(polizas$zona)),
  tipo_vehiculo = factor(
    "SUV",
    levels = levels(polizas$tipo_vehiculo)
  )
)

IC_media <- predict(
  modelo_diagnostico,
  newdata = nueva_poliza,
  interval = "confidence",
  level = 0.95
)

IP_individual <- predict(
  modelo_diagnostico,
  newdata = nueva_poliza,
  interval = "prediction",
  level = 0.95
)

IC_media

##        fit      lwr      upr
## 1 4178.718 3920.964 4436.473

IP_individual

##        fit      lwr      upr
## 1 4178.718 2729.104 5628.333

Comprobación 5. ¿Cuál intervalo debe reportarse para estimar el costo promedio de todas las pólizas con ese perfil y cuál para anticipar el costo de una nueva póliza individual?

Solución

Para el costo promedio se utiliza el intervalo de confianza de la respuesta media. Para una nueva póliza individual se utiliza el intervalo de predicción.

18 Explicación frente a predicción

La inferencia busca comprender parámetros, efectos ajustados e incertidumbre. La predicción busca desempeño sobre datos no utilizados en el ajuste.

Un modelo puede:

• explicar bien ciertos coeficientes y predecir modestamente;
• predecir bien, pero tener coeficientes difíciles de interpretar;
• presentar un \(R^2\) alto en la muestra y generalizar mal.

18.1 División entrenamiento-prueba

set.seed(1517)

indice_entrenamiento <- sample(
  seq_len(nrow(polizas)),
  size = floor(0.75 * nrow(polizas))
)

entrenamiento <- polizas[indice_entrenamiento, ]
prueba <- polizas[-indice_entrenamiento, ]

modelo_entrenamiento <- lm(
  costo_siniestros ~
    valor_vehiculo_10k +
    kilometraje_1000 +
    edad +
    zona +
    tipo_vehiculo,
  data = entrenamiento
)

pred_prueba <- predict(
  modelo_entrenamiento,
  newdata = prueba
)

RMSE <- sqrt(
  mean((prueba$costo_siniestros - pred_prueba)^2)
)

MAE <- mean(
  abs(prueba$costo_siniestros - pred_prueba)
)

c(RMSE = RMSE, MAE = MAE)

##     RMSE      MAE 
## 728.8596 567.8174

18.2 \(R^2\) fuera de muestra

SSE_prueba <- sum(
  (prueba$costo_siniestros - pred_prueba)^2
)

SST_prueba <- sum(
  (
    prueba$costo_siniestros -
    mean(entrenamiento$costo_siniestros)
  )^2
)

R2_prueba <- 1 - SSE_prueba / SST_prueba
R2_prueba

## [1] 0.5932393

La evaluación predictiva debe usar datos nuevos, validación cruzada o remuestreo. El \(R^2\) de entrenamiento es una medida de ajuste, no una garantía de generalización.

19 Estandarización y comparación de magnitudes

Los coeficientes no estandarizados dependen de las unidades. Para comparar asociaciones en escalas distintas puede estimarse un modelo con variables estandarizadas:

\[ Z=\frac{X-\overline{X}}{s_X}. \]

datos_estandarizados <- data.frame(
  costo_z = as.numeric(scale(polizas$costo_siniestros)),
  valor_z = as.numeric(scale(polizas$valor_vehiculo_10k)),
  km_z = as.numeric(scale(polizas$kilometraje_1000)),
  edad_z = as.numeric(scale(polizas$edad))
)

modelo_estandarizado <- lm(
  costo_z ~ valor_z + km_z + edad_z,
  data = datos_estandarizados
)

coef(modelo_estandarizado)

##              (Intercept)                  valor_z                     km_z                   edad_z 
## -0.000000000000000227079  0.759522869720220650258  0.206028106100576036130 -0.013317297462060586483

Un coeficiente estandarizado indica el cambio esperado en desviaciones estándar de \(Y\) por un aumento de una desviación estándar en el predictor, manteniendo los otros constantes.

Los coeficientes estandarizados no sustituyen la interpretación en unidades originales y no resuelven multicolinealidad, confusión ni problemas de diseño.

20 Construcción razonada del modelo

Una estrategia disciplinada:

1. Definir la variable dependiente.
2. Formular la pregunta sustantiva.
3. Clasificar predictores como exposición principal, controles, factores, mediadores o posibles modificadores de efecto.
4. Examinar rangos, unidades, datos faltantes y observaciones imposibles.
5. Dibujar relaciones bivariadas.
6. Especificar el modelo inicial antes de observar valores \(p\).
7. Respetar la jerarquía: si se incluye una interacción, mantener sus efectos principales.
8. Comparar modelos anidados cuando la hipótesis lo requiera.
9. Examinar supuestos e influencia.
10. Evaluar precisión, magnitud y utilidad práctica.
11. Validar predicciones si ese es el objetivo.
12. Reportar limitaciones.

20.1 Lo que no debe hacerse

• Seleccionar variables únicamente por significancia individual.
• Interpretar \(R^2\) como porcentaje de causalidad.
• eliminar automáticamente observaciones influyentes.
• Extrapolar lejos del rango observado.
• Reportar únicamente valores \(p\).
• Ignorar la codificación de factores.
• Interpretar efectos principales sin considerar una interacción relevante.
• Tratar la ausencia de significancia como demostración de ausencia de efecto.

21 Referencias utilizadas para la estructura conceptual

Anderson, D. R., Sweeney, D. J. y Williams, T. A. Estadística para negocios y economía, 11.ª edición. Capítulo 15, Regresión múltiple.

Rutherford, A. ANOVA and ANCOVA: A GLM Approach, 2.ª edición. Secciones sobre Modelo Lineal General, comparación de modelos completos y reducidos, codificación de condiciones experimentales, ANCOVA, supuestos y regresiones heterogéneas.

El material presenta una síntesis didáctica original y utiliza un caso simulado para fines docentes.