0.1 Comparación entre RStudio y Python, porque salen los mismos resultados:

Los resultados obtenidos en R y en Python coinciden porque ambos lenguajes ajustan el modelo de regresión lineal múltiple con el mismo método, el de mínimos cuadrados ordinarios, cuyo estimador β̂ = (XᵀX)⁻¹Xᵀy es único cuando las variables no son perfectamente colineales, y lo aplican sobre exactamente los mismos datos. Como la estimación de los coeficientes, la varianza residual σ̂², los errores estándar, el estadístico F global, el R² y los intervalos de confianza y de predicción (que se construyen con la misma distribución t de Student y los mismos grados de libertad) dependen únicamente de esos datos y de esas fórmulas, ambos programas deben entregar los mismos valores. Las únicas diferencias observables son de redondeo numérico, del orden de 10⁻¹² o menores, debidas a los algoritmos de álgebra lineal internos de cada implementación, y de formato de presentación, como la forma de mostrar p-valores muy pequeños (<2e-16 en R frente a 0.000 en Python). Por tanto, la coincidencia confirma que la implementación es correcta y que el análisis es reproducible independientemente del software utilizado.

1 Ejercicio 1. Publicidad y ventas

1.1 a) EDA

# --- Lectura y estructura ---
ventas <- read.table("ventas.txt", header = TRUE)

cat("Dimensiones:", dim(ventas), "\n")
## Dimensiones: 180 4
kable(head(ventas), caption = "Primeras filas del dataset (ventas)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Primeras filas del dataset (ventas)
tv radio periodico ventas
230.1 37.8 69.2 22.1
44.5 39.3 45.1 10.4
17.2 45.9 69.3 9.3
151.5 41.3 58.5 18.5
180.8 10.8 58.4 12.9
8.7 48.9 75.0 7.2 
# --- Datos faltantes ---
kable(t(colSums(is.na(ventas))), caption = "Valores faltantes por variable") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Valores faltantes por variable
tv radio periodico ventas
0 0 0 0 
# --- Descriptivos ---
desc1 <- data.frame(
  media = sapply(ventas, mean),
  sd    = sapply(ventas, sd),
  cv    = sapply(ventas, function(x) sd(x) / mean(x)),
  min   = sapply(ventas, min),
  max   = sapply(ventas, max)
)
kable(round(desc1, 4), caption = "Estadísticos descriptivos: Publicidad y Ventas") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Estadísticos descriptivos: Publicidad y Ventas
media sd cv min max
tv 145.2989 87.2465 0.6005 0.7 296.4
radio 24.0583 15.0130 0.6240 0.0 49.6
periodico 31.3228 22.0205 0.7030 0.3 114.0
ventas 14.0789 5.3873 0.3826 1.6 27.0 
library(ggcorrplot)

# --- Matriz de correlaciones (Gráfico) ---
ggcorrplot(cor(ventas),
           hc.order = FALSE,
           type = "lower",
           lab = TRUE,
           lab_size = 4.5,
           method = "square",
           colors = c("#2166ac", "#f7f7f7", "#b2182b"),
           title = "Matriz de Correlación: Publicidad y Ventas",
           ggtheme = theme_minimal())

# --- Dispersión ventas vs. cada predictor ---
par(mfrow = c(1, 3))
for (v in c("tv", "radio", "periodico")) {
  plot(ventas[[v]], ventas$ventas,
       xlab = paste(v, "(millones)"), ylab = "ventas (millones)",
       main = paste("ventas vs.", v), pch = 19, col = "steelblue")
  abline(lm(ventas$ventas ~ ventas[[v]]), col = "red", lwd = 2)
}

par(mfrow = c(1, 1))
# --- Atípicos: boxplots y conteo por regla 1.5*IQR ---
par(mfrow = c(1, 4))
for (v in names(ventas)) boxplot(ventas[[v]], main = v, col = "lightblue")

par(mfrow = c(1, 1))

atipicos1 <- t(as.data.frame(sapply(ventas, function(x) length(boxplot.stats(x)$out))))
rownames(atipicos1) <- "Nº atípicos"
kable(atipicos1, caption = "Número de valores atípicos por variable (regla 1.5·IQR)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Número de valores atípicos por variable (regla 1.5·IQR)
tv radio periodico ventas
Nº atípicos 0 0 2 0

Comentarios:

Vemos que el dataset tiene 180 observaciones donde la variable respuesta es ventas (expresado en millones) y las variables explicativas son el presupuesto destinado a la publicidad en la radio, periodico y ventas. Vemos que no hay ningun dato faltante y que la mayor media (por una diferencia mayor de 100 millones) la tiene el presupuesto destinado a la televisión con un poco mas de 145 millones. El rango de publicidad televisiva es el mayor también. Esto es lógico teniendo en cuenta que la televisión es una herramienta donde la publicidad está mas valorada que en las demás herramientas.

La matriz de correlaicones nos indica que hay correlación positiva entre todas las variables, con un grado mayor entre ventas y las variables explicativas comparando con la correlación entre las variables explicativas ensi. Lo principal que vemos en las graficas de dispersion es que todas las variables tienen una asociación lineal positiva. Esto es importante, no tendría sentido que invertir en publicidad se asocie con una disminución en ventas. La mayor pendiente lo tiene la televisión. Aparte de eso, vemos que en la radio y más en el periódico, la nube de puntos es muscho más dispersa y que con el mismo dinero invertido se obtienen unos resultados distintos, luego hay una menor correlación con ventas.

Cabe destacar que esto es analisis NO es la de la regresión múltiple y que solo coincidiría con el si los predictores no estuvieran correlaiconados entre si. Que como hemos visto no ocurre así.

En cuanto a valores atípicos, hay únicamente 2 y son las dos en la variable periódico.

1.2 b) Modelo de regresión múltiple

El modelo ajustado es

\[\text{ventas}_i = \beta_0 + \beta_1\,\text{tv}_i + \beta_2\,\text{radio}_i + \beta_3\,\text{periodico}_i + \varepsilon_i, \qquad \varepsilon_i \overset{iid}{\sim} N(0,\sigma^2).\]

mod1 <- lm(ventas ~ tv + radio + periodico, data = ventas)

# Tabla de coeficientes del modelo
coef_mod1 <- as.data.frame(summary(mod1)$coefficients)
names(coef_mod1) <- c("Estimación", "Error Estándar", "t", "p-valor")
kable(round(coef_mod1, 5), caption = "Coeficientes del modelo: ventas ~ tv + radio + periodico") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
  column_spec(5, bold = TRUE, color = ifelse(coef_mod1$`p-valor` < 0.05, "green", "red"))
Coeficientes del modelo: ventas ~ tv + radio + periodico
Estimación Error Estándar t p-valor
(Intercept) 2.76972 0.32929 8.41130 0.00000
tv 0.04716 0.00145 32.48482 0.00000
radio 0.18791 0.00894 21.02306 0.00000
periodico -0.00205 0.00609 -0.33572 0.73748 
# Estadísticos globales
gl_mod1 <- data.frame(
  `R²`          = summary(mod1)$r.squared,
  `R² ajustado` = summary(mod1)$adj.r.squared,
  `Sigma (RSE)` = summary(mod1)$sigma,
  `F global`    = summary(mod1)$fstatistic[1],
  check.names   = FALSE
)
kable(round(gl_mod1, 6), caption = "Bondad de ajuste global") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Bondad de ajuste global
R² R² ajustado Sigma (RSE) F global
value 0.90324 0.901591 1.690005 547.644

Interpretación:

  • Intercepto (β̂₀): Las ventas esperadas cuando tv, radio y periodico son 0 son 2.77. Esto indica que sin ninguna publicidad las ventas serían positivas (resultado intepretable por lógica).

  • β̂₁ (tv): manteniendo constantes radio y periodico, un aumento de 1 millón en la inversión en TV se asocia con un aumento de 47.160 € en las ventas esperadas.

  • β̂₂ (radio): manteniendo constantes televisión y periodico, un aumento de 1 millón en la inversión en TV se asocia con un aumento de 187.910 € en las ventas esperadas.

  • β̂₃ (periodico): manteniendo constantes televisión y radio, un aumento de 1 millón en la inversión en TV se asocia con una disminución de 2.050 € en las ventas esperadas. 

    NOTA IMPORTANTE: Este resultado no es estadísticamente significativo (a diferencia de las dos anteriores) ya que su p-valor es 0.737. Luego no deberíamos tomar en cuenta esta interpretación.

  • Contraste t individual (cada coeficiente \(\beta_j\)):

    \[H_0:\beta_j = 0 \quad \text{frente a} \quad H_1:\beta_j \neq 0, \qquad j = 1, 2, 3.\]

    Los p-valores de \(\hat\beta_1\) (tv) y \(\hat\beta_2\) (radio) son \(\approx 0\), por lo que se rechaza \(H_0\) para ambos. El de \(\hat\beta_3\) (periódico) es \(0.737 > 0.05\), no se rechaza \(H_0\): periódico no tiene efecto significativo.

  • Contraste F global:

    \[H_0:\beta_1 = \beta_2 = \beta_3 = 0 \quad \text{frente a} \quad H_1:\exists\, j \text{ tal que } \beta_j \neq 0.\]

    Con \(F = 547.6\) y p-valor \(\approx 0\), se rechaza \(H_0\): al menos una variable explica la variación de las ventas.

  • R² ajustado: el modelo explica aproximadamente 90.16% de la variabilidad de las ventas. El no ajustado tiene un valor de 90.32

    Las sumas de cuadrados se definen como:

    \[SST = \sum_{i=1}^{n}(y_i - \bar{y})^2, \qquad SSR = \sum_{i=1}^{n}(\hat{y}_i - \bar{y})^2, \qquad SSE = \sum_{i=1}^{n}(y_i - \hat{y}_i)^2\]

    y el coeficiente de determinación se obtiene como:

    \[R^2 = \frac{SSR}{SST} = 1 - \frac{SSE}{SST}\]

    Con los valores de este modelo (\(n = 180\), \(k = 3\) predictores, \(p = k+1 = 4\) parámetros):

    \[R^2 = \frac{4{,}692.40}{5{,}195.08} = 1 - \frac{502.68}{5{,}195.08} = 0.9032\]

    El \(R^2\) ajustado penaliza por el número de parámetros para evitar sobreajuste:

    \[R^2_{\text{adj}} = 1 - \frac{SSE/(n-p)}{SST/(n-1)}\]

    Con los valores concretos:

    \[R^2_{\text{adj}} = 1 - \frac{502.68\,/\,(180-4)}{5{,}195.08\,/\,(180-1)} = 1 - \frac{502.68/176}{5{,}195.08/179} = 1 - \frac{2.8561}{29.0228} = 0.9016\]

    El estadístico F global contrasta si algún predictor es útil:

    \[F = \frac{SSR/k}{SSE/(n-p)}\]

    Con los valores concretos (\(k=3\) predictores):

    \[F = \frac{SSR/k}{SSE/(n-p)} = \frac{4{,}692.40\,/\,3}{502.68\,/\,(180-4)} = \frac{1{,}564.13}{2.856} = 547.64\]

    Este valor de \(F\) tiene una distribución \(F_{k,\,n-p} = F_{3,\,176}\) bajo \(H_0\), dando un p-valor \(\approx 0\).

1.3 c) Intervalo de confianza del 95 % para la respuesta media

Se busca un intervalo de confianza para

\[E(\text{ventas} \mid \text{tv}=147.05,\ \text{radio}=23.3,\ \text{periodico}=31) = \beta_0 + 147.05\,\beta_1 + 23.3\,\beta_2 + 31\,\beta_3 .\]

x0 <- data.frame(tv = 147.05, radio = 23.3, periodico = 31)

ic_r <- predict(mod1, newdata = x0, interval = "confidence", level = 0.95)
ic_df <- data.frame(Estimación = ic_r[,"fit"], `Límite inf. (95%)` = ic_r[,"lwr"],
                    `Límite sup. (95%)` = ic_r[,"upr"], check.names = FALSE)
kable(round(ic_df, 5), caption = "Intervalo de confianza al 95 % para la respuesta media (tv=147.05, radio=23.3, periodico=31)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Intervalo de confianza al 95 % para la respuesta media (tv=147.05, radio=23.3, periodico=31)
Estimación Límite inf. (95%) Límite sup. (95%)
14.01963 13.77065 14.26861 
# --- Cálculo manual (verificación) ---
x0v   <- c(1, 147.05, 23.3, 31)
X     <- model.matrix(mod1)
s2    <- summary(mod1)$sigma^2
gl    <- df.residual(mod1)

est   <- sum(x0v * coef(mod1))
ee    <- sqrt(s2 * as.numeric(t(x0v) %*% solve(t(X) %*% X) %*% x0v))
tcrit <- qt(0.975, gl)

manual_ic <- data.frame(
  Estimación      = est,
  `Error Estándar` = ee,
  `t crítico`     = tcrit,
  `Límite inf.`   = est - tcrit * ee,
  `Límite sup.`   = est + tcrit * ee,
  check.names = FALSE
)
kable(round(manual_ic, 5), caption = "Verificación manual del IC al 95 %") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Verificación manual del IC al 95 %
Estimación Error Estándar t crítico Límite inf. Límite sup.
14.01963 0.12616 1.97353 13.77065 14.26861

Interpretación: con un nivel de confianza del 95 %, el ingreso medio por ventas de las regiones que invierten 147.05 millones en TV, 23.3 en radio y 31 en periódicos se estima entre 13.77 y 14.27 millones, con una estimación puntual de 14.02 millones. Es importante destcar que es un intervalo para la media de ventas. Queda además verificado manualmente. Además vemos que en Python el resultado es el mismo quitando las diferencias de redondeo.

1.4 d) Predicción para la nueva sucursal

Para una nueva sucursal con \(\text{tv}=160\), \(\text{radio}=26.7\) y \(\text{periodico}=35.1\) se busca la predicción puntual \(\widehat{\text{ventas}}\) y un intervalo de predicción del 95 % para el valor de sus ventas.

x1 <- data.frame(tv = 160, radio = 26.7, periodico = 35.1)

# Intervalo de predicción (lo que pide el apartado)
ip_r <- predict(mod1, newdata = x1, interval = "prediction", level = 0.95)
ip_df <- data.frame(Predicción = ip_r[,"fit"], `Límite inf. (95%)` = ip_r[,"lwr"],
                    `Límite sup. (95%)` = ip_r[,"upr"], check.names = FALSE)
kable(round(ip_df, 5), caption = "Intervalo de predicción al 95 % (tv=160, radio=26.7, periodico=35.1)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Intervalo de predicción al 95 % (tv=160, radio=26.7, periodico=35.1)
Predicción Límite inf. (95%) Límite sup. (95%)
15.26088 11.91571 18.60605 
# IC de la media en el mismo punto (para el apartado e)
ic_x1_r <- predict(mod1, newdata = x1, interval = "confidence", level = 0.95)

ancho_IP <- unname(ip_r[, "upr"] - ip_r[, "lwr"])
ancho_IC <- unname(ic_x1_r[, "upr"] - ic_x1_r[, "lwr"])

comp_df <- data.frame(
  Intervalo    = c("IC (confianza)", "IP (predicción)"),
  Estimación   = c(ic_x1_r[,"fit"], ip_r[,"fit"]),
  `Lím. inf.`  = c(ic_x1_r[,"lwr"], ip_r[,"lwr"]),
  `Lím. sup.`  = c(ic_x1_r[,"upr"], ip_r[,"upr"]),
  Anchura      = c(ancho_IC, ancho_IP),
  check.names  = FALSE
)
comp_df[, sapply(comp_df, is.numeric)] <- round(comp_df[, sapply(comp_df, is.numeric)], 5)
kable(comp_df, caption = "Comparación IC vs. IP en el mismo punto (tv=160, radio=26.7, periodico=35.1)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Comparación IC vs. IP en el mismo punto (tv=160, radio=26.7, periodico=35.1)
Intervalo Estimación Lím. inf. Lím. sup. Anchura
IC (confianza) 15.26088 15.00383 15.51792 0.51409
IP (predicción) 15.26088 11.91571 18.60605 6.69034 
# --- Cálculo manual (verificación) ---
x1v   <- c(1, 160, 26.7, 35.1)
X     <- model.matrix(mod1)
s2    <- summary(mod1)$sigma^2
gl    <- df.residual(mod1)
h     <- as.numeric(t(x1v) %*% solve(t(X) %*% X) %*% x1v)

est   <- sum(x1v * coef(mod1))
ee_p  <- sqrt(s2 * (1 + h))
tcrit <- qt(0.975, gl)

manual_ip <- data.frame(
  Predicción       = est,
  `Error Estándar` = ee_p,
  `t crítico`      = tcrit,
  `Límite inf.`    = est - tcrit * ee_p,
  `Límite sup.`    = est + tcrit * ee_p,
  check.names = FALSE
)
kable(round(manual_ip, 5), caption = "Verificación manual del IP al 95 %") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Verificación manual del IP al 95 %
Predicción Error Estándar t crítico Límite inf. Límite sup.
15.26088 1.69502 1.97353 11.91571 18.60605

Interpretación: La predicción puntual de ventas para la nueva sucursal es de 15.261 millones. Con un 95 % de confianza, se espera que las ventas de esa sucursal concreta se encuentren entre 11.92 y 18.61 millones. Este intervalo se refiere al valor de una observación nueva, no a la media de ventas de un grupo de regiones, por lo que es más ancho que el intervalo de confianza del apartado c). Vemos que la diferencia entre anchuras es más de 6.1 millones de euros.

NOTA: Vemos que la función predict es la misma para los dos pero interval en el de confianza toma el valor “confidence” y en la de la predicción “prediction”.

1.5 e) Diferencia conceptual entre el IP y IC.

Los intervalos de los apartados c) y d) responden a preguntas distintas.

  • Intervalo de confianza (apartado c).

    Estima \(E(Y \mid X = x_0) = x_0^\top \beta\), la media de ventas de todas las regiones con esa inversión. Es un parámetro fijo, aunque desconocido; el intervalo cuantifica solo la incertidumbre por haber estimado \(\beta\) con una muestra finita. Su varianza es \[\operatorname{Var}(\hat{Y}_0) = \sigma^2\, x_0^\top (X^\top X)^{-1} x_0 .\]

  • Intervalo de predicción (apartado d). Cubre \(Y_{\text{nuevo}} \mid X = x_0\), el valor de una concreta observación futura, que es una variable aleatoria: \(Y_{\text{nuevo}} = x_0^\top \beta + \varepsilon_{\text{nuevo}}\). El error de predicción \(Y_{\text{nuevo}} - \hat{Y}_0\) tiene dos fuentes independientes de variabilidad: la incertidumbre al estimar la recta y el error aleatorio \(\varepsilon_{\text{nuevo}}\) de la nueva observación. Por eso \[\operatorname{Var}(Y_{\text{nuevo}} - \hat{Y}_0) = \sigma^2 + \sigma^2\, x_0^\top (X^\top X)^{-1} x_0 = \sigma^2\left[1 + x_0^\top (X^\top X)^{-1} x_0\right].\]

Ambos intervalos comparten el término \(\sigma^2 x_0^\top (X^\top X)^{-1} x_0\), pero el de predicción suma además \(\sigma^2\), la variabilidad inevitable de una observación individual alrededor de la media, ya que es una variable aleatoria. Como esa variabilidad no desaparece aunque se conozca perfectamente la recta, el IP es siempre más ancho que el IC en el mismo punto. Además, al aumentar el tamaño muestral el término \(x_0^\top (X^\top X)^{-1} x_0\) tiende a 0, de modo que el IC se contrae hacia un punto, mientras que el IP se estabiliza en torno a \(\pm t\,\sigma\).

2 Ejercicio 2. Biomasa y características del suelo

2.1 a) EDA

suelo <- read.table("datos.txt", header = TRUE)

cat("Dimensiones:", dim(suelo), "\n")
## Dimensiones: 40 6
kable(head(suelo), caption = "Primeras filas del dataset (suelo)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Primeras filas del dataset (suelo)
biomasa salinidad pH K Na Zn
676 33 5.00 1441.67 35184.5 16.4524
516 35 4.75 1299.19 28170.4 13.9852
1052 32 4.20 1154.27 26455.0 15.3276
868 30 4.40 1045.15 25072.9 17.3128
436 33 5.05 1273.02 25491.7 12.2778
544 36 4.25 1346.35 20877.3 17.8225

Vemos como tenemos un dataset de 40 observaciones con 5 variables explicativas: salinidad(%), pH, potasio(ppm), Sodio(ppm) y Zinc(ppm) que se utilizan para identificar características del suelo que pueden influir sobre la biomasa Y de una especie de pasto. Cabe destacar como y ahe indicado en paentesis que las variables están medidos en unidades distintos y no comparables sin una normalización previa.

# --- Datos faltantes ---
kable(t(colSums(is.na(suelo))), caption = "Valores faltantes por variable") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Valores faltantes por variable
biomasa salinidad pH K Na Zn
0 0 0 0 0 0 
# --- Descriptivos ---
desc <- data.frame(
  media = sapply(suelo, mean),
  sd    = sapply(suelo, sd),
  cv    = sapply(suelo, function(x) sd(x) / mean(x)),
  min   = sapply(suelo, min),
  max   = sapply(suelo, max)
)
kable(round(desc, 3), caption = "Estadísticos descriptivos: Características del Suelo") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Estadísticos descriptivos: Características del Suelo
media sd cv min max
biomasa 996.900 681.755 0.684 236.00 2436.000
salinidad 30.175 3.679 0.122 24.00 38.000
pH 4.541 1.219 0.268 3.20 7.450
K 804.364 293.826 0.365 350.73 1441.670
Na 16213.448 6606.895 0.407 7886.50 35184.500
Zn 18.176 8.220 0.452 0.21 31.286

Vemos como no hay ningun valor faltante (na). Además de eso, lo mas significativo es que el mayor rango lo tiene el sodio (Na).

library(ggcorrplot)

# --- Matriz de correlaciones (Gráfico) ---
ggcorrplot(cor(suelo),
           hc.order = FALSE,
           type = "lower",
           lab = TRUE,
           lab_size = 4,
           method = "square",
           colors = c("#2166ac", "#f7f7f7", "#b2182b"),
           title = "Matriz de Correlación: Características del Suelo y Biomasa",
           ggtheme = theme_minimal())

En cuanto a las correlaciones, en este dataset vemos que varían mucho más que en la anterior. Los más signficativos (con una relación mas fuerte) son el sodio y el potasio, con una correlación del 0.9, la biomasa y el pH, con 0.82, lo cual significa que ap riori será un estimado rmuy bueno por si solo, y por último, vemos que el Zinc tuene una relación negativa bastante alta tanto con la biomasa, la salinidad y el pH. Viendo estas correlaciones, se cofirma lo que he dicho anteriormente.

2.2 Graficos entre Y X2 y X3 (Biomasa pH y K (potasio))

library(ggplot2)
library(patchwork)

g1 <- ggplot(suelo, aes(pH, biomasa)) +
  geom_point(color = "steelblue", alpha = 0.7) +
  geom_smooth(method = "lm", color = "red") +
  labs(title = "biomasa vs. pH", x = "pH", y = "biomasa")

g2 <- ggplot(suelo, aes(K, biomasa)) +
  geom_point(color = "steelblue", alpha = 0.7) +
  geom_smooth(method = "lm", color = "red") +
  labs(title = "biomasa vs. K (potasio)", x = "K (ppm)", y = "biomasa")

g3 <- ggplot(suelo, aes(pH, K)) +
  geom_point(color = "steelblue", alpha = 0.7) +
  geom_smooth(method = "lm", color = "red") +
  labs(title = "K vs. pH", x = "pH", y = "K (ppm)")

g1 + g2 + g3

# Matriz de dispersión de las tres variables de interés
pairs(suelo[, c("biomasa", "pH", "K")], pch = 19, col = "steelblue")

# Correlaciones de interés para la interpretación
cor_clave <- data.frame(
  Par              = c("biomasa vs. pH", "biomasa vs. K", "pH vs. K"),
  `Correlación (r)` = c(cor(suelo$biomasa, suelo$pH),
                         cor(suelo$biomasa, suelo$K),
                         cor(suelo$pH, suelo$K)),
  check.names = FALSE
)
cor_clave[, sapply(cor_clave, is.numeric)] <- round(cor_clave[, sapply(cor_clave, is.numeric)], 5)
kable(cor_clave, caption = "Correlaciones clave entre variables") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Correlaciones clave entre variables
Par Correlación (r)
biomasa vs. pH 0.81871
biomasa vs. K -0.14870
pH vs. K 0.09228

Interpretación:

  • biomasa vs. pH (Y vs. X2): la correlación es 0.82, lo que indica una relación positiva fuerte entre el pH del suelo y la biomasa.
  • biomasa vs. K (Y vs. X3): la correlación es -0.15, lo que indica una relación negativa débil entre el contenido del pH del suelo y la biomasa.
  • pH vs. K (X2 vs. X3): la correlación entre ambos predictores es 0.09
  • NOTA IMPORTANTE: Esto es importante ya que aunque la correlación entre los dos predictores sea bajo, el efecto de K sobre la biomasa cambiará al controlar por pH si su efecto parcial es MUY grande. Ya evremos en los siguientes apartados si se cumple esto o no se cumple.

2.3 b) Regresión lineal simple: biomasa ~ pH

Se ajusta el modelo

\[\text{biomasa}_i = \beta_0 + \beta_2\,\text{pH}_i + \varepsilon_i, \qquad \varepsilon_i \overset{iid}{\sim} N(0,\sigma^2).\]

m_simple <- lm(biomasa ~ pH, data = suelo)

# Tabla de coeficientes
coef_s <- as.data.frame(summary(m_simple)$coefficients)
names(coef_s) <- c("Estimación", "Error Estándar", "t", "p-valor")
kable(round(coef_s, 6), caption = "Coeficientes del modelo simple: biomasa ~ pH") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
  column_spec(5, bold = TRUE, color = ifelse(coef_s$`p-valor` < 0.05, "green", "red"))
Coeficientes del modelo simple: biomasa ~ pH
Estimación Error Estándar t p-valor
(Intercept) -1083.2272 244.83514 -4.424313 7.9e-05
pH 458.0517 52.11537 8.789186 0.0e+00 
# --- Coeficientes con IC 95 % ---
tabla_coef <- cbind(
  Estimación      = coef(m_simple),
  confint(m_simple),
  `p-valor`       = summary(m_simple)$coefficients[, 4]
)
kable(round(tabla_coef, 6), caption = "Coeficientes e IC 95 %: biomasa ~ pH") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Coeficientes e IC 95 %: biomasa ~ pH
Estimación 2.5 % 97.5 % p-valor
(Intercept) -1083.2272 -1578.8700 -587.5844 7.9e-05
pH 458.0517 352.5496 563.5537 0.0e+00 
# --- Ajuste global ---
f <- summary(m_simple)$fstatistic
gl_s <- data.frame(
  `R²`          = summary(m_simple)$r.squared,
  `R² ajustado` = summary(m_simple)$adj.r.squared,
  `Sigma (RSE)` = summary(m_simple)$sigma,
  `F global`    = unname(f[1]),
  `p-valor F`   = pf(f[1], f[2], f[3], lower.tail = FALSE),
  check.names   = FALSE
)
kable(round(gl_s, 6), caption = "Bondad de ajuste global: biomasa ~ pH") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Bondad de ajuste global: biomasa ~ pH
R² R² ajustado Sigma (RSE) F global p-valor F
value 0.670281 0.661605 396.5888 77.24979 0 
library(ggplot2)

ggplot(suelo, aes(pH, biomasa)) +
  geom_point(color = "steelblue", alpha = 0.7) +
  geom_smooth(method = "lm", color = "red") +
  labs(title = "Recta ajustada: biomasa ~ pH", x = "pH", y = "biomasa")

Interpretación:

  • Pendiente (β̂₂): por cada unidad adicional de pH, con un signo positiv, la biomasa esperada cambia en 458.03 unidades. Esto indica que suelos con mayor pH se asocian con mayor biomasa.

  • Significancia de la pendiente:

    \[H_0:\beta_2 = 0 \quad \text{frente a} \quad H_1:\beta_2 \neq 0.\]

    El p-valor es \(1.09 \times 10^{-10} \ll 0.05\), por lo que se rechaza \(H_0\): el pH tiene un efecto significativo sobre la biomasa.

  • Intercepto (β̂₀): biomasa esperada con pH = 0 que da -1083.23, un valor fuera del rango observado, por lo que no tiene interpretación práctica; solo ancla la recta.

  • Ajuste: R² = 0.67 (el ajustado 0.66), es decir, el pH explica aproximadamente el 67% de la variabilidad de la biomasa.

    Aplicando la fórmula general con los valores de este modelo (\(n = 40\), \(k = 1\), \(p = 2\)):

    \[R^2 = \frac{12{,}150{,}057}{18{,}126{,}800} = 1 - \frac{5{,}976{,}743}{18{,}126{,}800} = 0.6703\]

    El \(R^2\) ajustado:

    \[R^2_{\text{adj}} = 1 - \frac{SSE/(n-p)}{SST/(n-1)} = 1 - \frac{5{,}976{,}743\,/\,(40-2)}{18{,}126{,}800\,/\,(40-1)} = 1 - \frac{157{,}283.8}{464{,}789.7} = 0.6616\]

    El estadístico F global (en regresión simple coincide con \(t^2\) de la pendiente):

    \[F = \frac{SSR/k}{SSE/(n-p)} = \frac{12{,}150{,}057\,/\,1}{5{,}976{,}743\,/\,(40-2)} = \frac{12{,}150{,}057}{157{,}283.8} = 77.25\]

    Bajo \(H_0\) este estadístico sigue una \(F_{1,\,38}\), con p-valor \(\approx 1.09 \times 10^{-10}\).

  • Aviso importante: este coeficiente es un efecto marginal: no controla por K ni por las demás variables del suelo. Puede cambiar al incorporar otras variables como se hace con la k en el apartado e de este mismo ejercicio.

2.4 c) Residuos y control por pH

En el apartado b) hemos ajustado biomasa ~ pH y obtenido un coeficiente marginal: ese \(\hat\beta_2 = 458\) recoge no solo el efecto directo del pH sino todo lo que K pueda tener en común con el pH. El objetivo de este apartado es aislar el efecto de K sobre la biomasa descontando el pH, y hacerlo de forma visual y conceptual antes de llegar al modelo múltiple formal.

La idea clave es la siguiente: tanto la biomasa como K están parcialmente “contaminadas” por el pH. Si queremos ver la relación pura entre biomasa y K, debemos quitarle a ambas variables lo que el pH ya explica de ellas. Para eso se construyen dos regresiones auxiliares y se toman sus residuos:

  • Primera regresión auxiliar: biomasa ~ pH \(\Rightarrow\) residuos \(e_Y\) = parte de la biomasa no explicada por el pH.
  • Segunda regresión auxiliar: K ~ pH \(\Rightarrow\) residuos \(e_{X_3}\) = parte de K no explicada por el pH.

Una vez que ambas variables están “depuradas” del pH, la relación entre ellas refleja únicamente la información que K aporta por encima de lo que ya recoge el pH.

# (i) Regresión biomasa ~ pH  -> residuos e_Y
mod_Y  <- lm(biomasa ~ pH, data = suelo)
eY     <- resid(mod_Y)

# (ii) Regresión K ~ pH  -> residuos e_X3
mod_X3 <- lm(K ~ pH, data = suelo)
eX3    <- resid(mod_X3)

# Coeficientes de ambas regresiones auxiliares
coef_Y  <- as.data.frame(summary(mod_Y)$coefficients)
coef_X3 <- as.data.frame(summary(mod_X3)$coefficients)
names(coef_Y) <- names(coef_X3) <- c("Estimación", "Error Estándar", "t", "p-valor")

kable(round(coef_Y, 6),  caption = "Coeficientes: biomasa ~ pH") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Coeficientes: biomasa ~ pH
Estimación Error Estándar t p-valor
(Intercept) -1083.2272 244.83514 -4.424313 7.9e-05
pH 458.0517 52.11537 8.789186 0.0e+00 
kable(round(coef_X3, 6), caption = "Coeficientes: K ~ pH") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Coeficientes: K ~ pH
Estimación Error Estándar t p-valor
(Intercept) 703.31236 182.98112 3.843634 0.000448
pH 22.25183 38.94918 0.571304 0.571157 
# Comprobaciones: residuos con media 0 y ortogonales a pH
check_df <- data.frame(
  `Media eY`    = mean(eY),
  `Media eX3`   = mean(eX3),
  `cor(eY,pH)`  = cor(eY, suelo$pH),
  `cor(eX3,pH)` = cor(eX3, suelo$pH),
  check.names   = FALSE
)
kable(round(check_df, 10), caption = "Verificación: residuos ortogonales a pH y con media 0") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Verificación: residuos ortogonales a pH y con media 0
Media eY Media eX3 cor(eY,pH) cor(eX3,pH)
0 0 0 0 
kable(round(head(data.frame(eY, eX3)), 5), caption = "Primeras filas de los residuos") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Primeras filas de los residuos
eY eX3
-531.0312 627.0985
-576.5183 490.1814
211.4101 357.4999
-64.2002 243.9296
-793.9338 457.3359
-319.4925 548.4674 
# (iii) Gráfico e_Y vs e_X3
library(ggplot2)

df_res <- data.frame(eX3 = eX3, eY = eY)

ggplot(df_res, aes(eX3, eY)) +
  geom_point(color = "steelblue", alpha = 0.7) +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(title = "Residuos de biomasa ~ pH  vs.  residuos de K ~ pH",
       x = "e_X3 (K libre de pH)", y = "e_Y (biomasa libre de pH)")

# Correlación entre residuos (= correlación parcial entre biomasa y K dado pH)
cor_parcial <- data.frame(
  `Correlación parcial biomasa-K dado pH` = cor(eY, eX3),
  check.names = FALSE
)
kable(round(cor_parcial, 6), caption = "Correlación parcial entre biomasa y K controlando por pH") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Correlación parcial entre biomasa y K controlando por pH
Correlación parcial biomasa-K dado pH
-0.392211

(iv) Qué representa cada conjunto de residuos e interpretación completa:

  • \(e_Y\) (residuos de biomasa ~ pH): representan la biomasa que no está explicada por el pH. Formalmente, para cada parcela \(i\), \(e_{Y,i} = \text{biomasa}_i - \hat{\text{biomasa}}_i(\text{pH})\). Un valor positivo indica que esa parcela tiene más biomasa de lo que su pH predeciría; uno negativo, menos. Esta variable ya no lleva información del pH.

  • \(e_{X_3}\) (residuos de K ~ pH): representan la parte del potasio que no está correlacionada con el pH. Es la variación en K que es independiente del pH, la única parte que puede aportar nueva información al modelo una vez que el pH ya está incluido.

  • Verificación de ortogonalidad: los residuos MCO tienen siempre media exactamente cero y son ortogonales (no correlacionados) a los predictores del modelo. Confirmar que \(\bar{e}_Y \approx 0\), \(\bar{e}_{X_3} \approx 0\) y que \(\text{cor}(e_Y, \text{pH}) \approx 0\), \(\text{cor}(e_{X_3}, \text{pH}) \approx 0\) es sólo una verificación numérica de que la construcción es correcta.

  • Gráfico \(e_Y\) vs. \(e_{X_3}\): este gráfico es el central del apartado. La pendiente de la recta ajustada sobre esta nube de puntos mide el efecto de K sobre la biomasa una vez eliminada la influencia compartida con el pH, es decir, es el efecto parcial de K. La correlación entre ambos residuos es la correlación parcial de biomasa y K dado pH. Veremos en el apartado e) que esa pendiente es exactamente \(\hat\beta_3\) del modelo múltiple (Teorema de Frisch-Waugh-Lovell).

  • Dirección de la nube: la nube tiene pendiente negativa, lo que indica que los suelos con más potasio del que les correspondería por su pH tienden a tener menos biomasa de la que les correspondería por su pH. Dicho de otro modo, el efecto parcial de K sobre la biomasa es negativo: dado el mismo pH, un mayor contenido en K se asocia con menor biomasa.

  • Propósito global del apartado: este procedimiento de residualización sirve para aislar visualmente el efecto parcial de una variable, y es la demostración intuitiva de por qué el coeficiente de K en el modelo múltiple puede ser muy distinto de su correlación simple con la biomasa. No se trata simplemente de calcular residuos: se trata de descomponer la variación de cada variable en una parte explicada por el pH y una parte nueva, y luego relacionar esas partes nuevas entre sí.

(v) Visualización 3D interactiva de los residuos:

Puedes escanear el siguiente código QR con tu móvil para acceder a una visualización interactiva en 3D que hemos preparado. En ella se aprecia el plano de regresión múltiple y las líneas verticales negras que caen desde cada punto real hasta el plano (los residuos ortogonales).

# Si no tienes instalado el paquete qrcode, ejecuta en la consola: install.packages("qrcode")
library(qrcode)

# NOTA: Debes subir el archivo Vis3D.html (que acabo de crear en tu carpeta) a alguna web 
# gratuita como RPubs o GitHub Pages. Cuando tengas el enlace, ponlo aquí abajo:
url_visualizacion <- "https://rpubs.com/tu_usuario/visualizacion_3D" 

qr <- qr_code(url_visualizacion)
plot(qr)

2.5 d) Regresión de residuos sobre residuos

Se ajusta el modelo

\[e_{Y,i} = \alpha + \gamma\, e_{X_3,i} + u_i ,\]

donde \(e_Y\) son los residuos de biomasa ~ pH y \(e_{X_3}\) los de K ~ pH.

mod_res <- lm(eY ~ eX3)

# Tabla de coeficientes del modelo de residuos
coef_res <- as.data.frame(summary(mod_res)$coefficients)
names(coef_res) <- c("Estimación", "Error Estándar", "t", "p-valor")
kable(round(coef_res, 8), caption = "Coeficientes del modelo de residuos: eY ~ eX3") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
  column_spec(5, bold = TRUE, color = ifelse(coef_res$`p-valor` < 0.05, "green", "red"))
Coeficientes del modelo de residuos: eY ~ eX3
Estimación Error Estándar t p-valor
(Intercept) 0.0000000 57.6818902 0.000000 1.0000000
eX3 -0.5247918 0.1996662 -2.628345 0.0123111 
gamma_hat <- unname(coef(mod_res)["eX3"])

# Tabla de coeficientes con IC 95%
tabla_res <- cbind(
  Estimación = coef(mod_res),
  confint(mod_res),
  `p-valor`  = summary(mod_res)$coefficients[, 4]
)
kable(round(tabla_res, 8), caption = "Coeficientes e IC 95 %: regresión de residuos") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Coeficientes e IC 95 %: regresión de residuos
Estimación 2.5 % 97.5 % p-valor
(Intercept) 0.0000000 -116.770882 116.7708818 1.0000000
eX3 -0.5247918 -0.928995 -0.1205886 0.0123111

Interpretación:

  • γ̂ = -0.52: pendiente de la nube de puntos \(e_Y\) vs. \(e_{X_3}\). El cambio esperado en la biomasa por cada unidad adicional de K una vez descontado el efecto del pH en ambas variables es de -0.52 unidades.

  • Intercepto α̂ ≈ 0: es cero lógicamente (salvo error numérico) porque los dos conjuntos de residuos tienen media exactamente cero.

  • Significancia de \(\hat\gamma\):

    \[H_0:\gamma = 0 \quad \text{frente a} \quad H_1:\gamma \neq 0.\]

    El p-valor de \(\hat\gamma\) es 0.0123. Fijando \(\alpha = 0.05\), hay evidencia de que K aporta información sobre la biomasa más allá del pH.

  • Nota importante (inferencia): el error estándar y el p-valor de γ̂ que reporta esta regresión no son los del modelo múltiple.

2.6 e) Regresión múltiple y equivalencia con \(\hat\gamma\)

Se ajusta el modelo

\[\text{biomasa}_i = \beta_0 + \beta_2\,\text{pH}_i + \beta_3\,K_i + \varepsilon_i .\]

m_mult <- lm(biomasa ~ pH + K, data = suelo)

# Tabla de coeficientes
coef_m <- as.data.frame(summary(m_mult)$coefficients)
names(coef_m) <- c("Estimación", "Error Estándar", "t", "p-valor")
kable(round(coef_m, 6), caption = "Coeficientes del modelo múltiple: biomasa ~ pH + K") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
  column_spec(5, bold = TRUE, color = ifelse(coef_m$`p-valor` < 0.05, "green", "red"))
Coeficientes del modelo múltiple: biomasa ~ pH + K
Estimación Error Estándar t p-valor
(Intercept) -714.134638 268.973751 -2.655035 0.011626
pH 469.729258 48.791357 9.627305 0.000000
K -0.524792 0.202346 -2.593531 0.013530 
# Bondad de ajuste
gl_m <- data.frame(
  `R²`          = summary(m_mult)$r.squared,
  `R² ajustado` = summary(m_mult)$adj.r.squared,
  `Sigma (RSE)` = summary(m_mult)$sigma,
  check.names   = FALSE
)
kable(round(gl_m, 6), caption = "Bondad de ajuste: biomasa ~ pH + K") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Bondad de ajuste: biomasa ~ pH + K
R² R² ajustado Sigma (RSE)
0.721002 0.705921 369.7093 
beta3_hat <- unname(coef(m_mult)["K"])
gamma_hat <- unname(coef(mod_res)["eX3"])

# Comparación con todos los decimales
print(beta3_hat, digits = 15)
## [1] -0.524791779453377
print(gamma_hat, digits = 15)
## [1] -0.524791779453377
c(beta3 = beta3_hat,
  gamma = gamma_hat,
  diferencia = beta3_hat - gamma_hat)
##         beta3         gamma    diferencia 
## -5.247918e-01 -5.247918e-01  2.220446e-16
# Verificación numérica formal
all.equal(beta3_hat, gamma_hat)
## [1] TRUE
isTRUE(all.equal(beta3_hat, gamma_hat, tolerance = 1e-10))
## [1] TRUE
# Tabla comparativa
tabla_e <- data.frame(
  procedimiento = c("Modelo múltiple: beta3 (K)",
                    "Residuos sobre residuos: gamma"),
  estimacion    = c(beta3_hat, gamma_hat),
  error_est     = c(summary(m_mult)$coefficients["K", 2],
                    summary(mod_res)$coefficients["eX3", 2]),
  gl_residuales = c(df.residual(m_mult), df.residual(mod_res))
)
knitr::kable(tabla_e, digits = 8,
             caption = "Coincide la estimación; difieren error estándar y grados de libertad")
Coincide la estimación; difieren error estándar y grados de libertad
procedimiento estimacion error_est gl_residuales
Modelo múltiple: beta3 (K) -0.5247918 0.2023465 37
Residuos sobre residuos: gamma -0.5247918 0.1996662 38

Este es el resultado más significativo que podemos obtener.

Aplicando la fórmula con los valores del modelo múltiple (\(n = 40\), \(k = 2\), \(p = 3\)):

\[R^2 = \frac{SSR}{SST} = 1 - \frac{SSE}{SST} = \frac{13{,}069{,}455}{18{,}126{,}800} = 1 - \frac{5{,}057{,}345}{18{,}126{,}800} = 0.7210\]

El \(R^2\) ajustado (penaliza por usar \(k=2\) predictores en vez de 1):

\[R^2_{\text{adj}} = 1 - \frac{SSE/(n-p)}{SST/(n-1)} = 1 - \frac{5{,}057{,}345\,/\,(40-3)}{18{,}126{,}800\,/\,(40-1)} = 1 - \frac{136{,}685.0}{464{,}789.7} = 0.7059\]

El estadístico F global para el contraste \(H_0:\beta_2 = \beta_3 = 0\) frente a \(H_1:\exists\,j,\,\beta_j\neq 0\):

\[F = \frac{SSR/k}{SSE/(n-p)} = \frac{13{,}069{,}455\,/\,2}{5{,}057{,}345\,/\,(40-3)} = \frac{6{,}534{,}727.5}{136{,}685.0} = 47.81\]

Bajo \(H_0\) este estadístico sigue una \(F_{2,\,37}\), con p-valor \(\approx 0\).

Nótese que el \(SST = 18{,}126{,}800\) es el mismo que en el modelo simple del apartado b) porque la variable respuesta (biomasa) no cambia. El \(SSR\) aumenta de \(12{,}150{,}057\) a \(13{,}069{,}455\) al incorporar \(K\), confirmando que esta variable explica variabilidad adicional que el pH solo no recogía.

# Comparación con la regresión simple biomasa ~ K (apoya el apartado f)
m_simple_K <- lm(biomasa ~ K, data = suelo)
fwl_comp <- data.frame(
  Modelo             = c("Múltiple (biomasa ~ pH + K)", "Simple (biomasa ~ K)"),
  `Coef. K estimado` = c(beta3_hat, unname(coef(m_simple_K)["K"])),
  check.names = FALSE
)
fwl_comp[, sapply(fwl_comp, is.numeric)] <- round(fwl_comp[, sapply(fwl_comp, is.numeric)], 8)
kable(fwl_comp, caption = "Comparación coeficiente K: múltiple vs simple") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Comparación coeficiente K: múltiple vs simple
Modelo Coef. K estimado
Múltiple (biomasa ~ pH + K) -0.5247918
Simple (biomasa ~ K) -0.3450211

Interpretación:

  • β̂₃ = -0.52 y γ̂ = -0.52: ambas estimaciones coinciden numéricamente (la diferencia es del orden de 1e-15),
  • Nota Importante: Informandome en internet, he visto que esto verifica el teorema de Frisch-Waugh-Lovell: el coeficiente de \(K\) en la regresión múltiple es la pendiente de regresar los residuos de biomasa ~ pH sobre los residuos de K ~ pH.
  • Qué es lo que coincide y qué no: coincide la estimación puntual; el error estándar y los grados de libertad difieren (n − 3 en el modelo múltiple, n − 2 en la regresión de residuos), por lo que el p-valor de γ̂ del apartado d) no debe usarse para inferencia. La inferencia correcta sobre β₃ es la del modelo múltiple. Esto es importante ya que los interavlos de confianza cunado difieren el error estandar y los grados de libertad no son los mismos para dos modelos que dan la misma predicción puntual.
  • Coincidencia entre lenguajes: R y Python dan los mismos valor, salvo diferencias de redondeo, como lo estamos viendo en todos los ejerciicios aunque en algunas no lo escriba.

2.7 f) ¿Por qué \(\hat\beta_3\) no es el coeficiente de biomasa ~ K?

m_simple_K <- lm(biomasa ~ K, data = suelo)
m_mult     <- lm(biomasa ~ pH + K, data = suelo)
m_aux      <- lm(K ~ pH, data = suelo)          # regresión auxiliar de K sobre pH

b_simple   <- unname(coef(m_simple_K)["K"])     # efecto marginal de K
b_multiple <- unname(coef(m_mult)["K"])         # efecto parcial de K (controlando pH)
b_pH       <- unname(coef(m_mult)["pH"])        # efecto parcial de pH
d_aux      <- unname(coef(m_aux)["pH"])         # pendiente de K sobre pH

tabla_f <- data.frame(
  coeficiente = c("K en biomasa ~ K (simple)",
                  "K en biomasa ~ pH + K (múltiple)",
                  "pH en biomasa ~ pH + K (múltiple)",
                  "pH en K ~ pH (auxiliar)"),
  valor = c(b_simple, b_multiple, b_pH, d_aux)
)
knitr::kable(tabla_f, digits = 6, caption = "Ingredientes de la comparación")
Ingredientes de la comparación
coeficiente valor
K en biomasa ~ K (simple) -0.345021
K en biomasa ~ pH + K (múltiple) -0.524792
pH en biomasa ~ pH + K (múltiple) 469.729258
pH en K ~ pH (auxiliar) 22.251835 
# Identidad del sesgo por variable omitida (exacta en MCO):
#   beta_simple(K) = beta3 + beta2 * (pendiente de pH sobre K)
m_aux2 <- lm(pH ~ K, data = suelo)
d2     <- unname(coef(m_aux2)["K"])

sesgo_df <- data.frame(
  Cantidad             = c("β̂K (simple)", "β̂3 + sesgo", "Sesgo (β̂2 · δ̂)", "Diferencia"),
  Valor                = c(b_simple,
                           b_multiple + b_pH * d2,
                           b_pH * d2,
                           b_simple - (b_multiple + b_pH * d2))
)
sesgo_df[, sapply(sesgo_df, is.numeric)] <- round(sesgo_df[, sapply(sesgo_df, is.numeric)], 8)
kable(sesgo_df, caption = "Verificación numérica de la identidad del sesgo por variable omitida") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Verificación numérica de la identidad del sesgo por variable omitida
Cantidad Valor
β̂K (simple)
-0.345021
 
# Cuánto se relacionan los predictores (mecanismo del sesgo)
cor_ph_k <- data.frame(`cor(pH, K)` = cor(suelo$pH, suelo$K), check.names = FALSE)
kable(round(cor_ph_k, 6), caption = "Correlación entre pH y K") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Correlación entre pH y K
cor(pH, K)
0.092282

\(\beta_3\) es el efecto parcial de \(K\) sobre la biomasa: el cambio esperado en la biomasa por cada unidad adicional de potasio manteniendo constante el pH. Por otro lado el coeficiente de la regresión simple biomasa ~ K es un efecto marginal: compara suelos con distinto potasio sin fijar el pH, de modo que también recoge parte del efecto del pH siempre que ambas variables estén relacionadas.

La relacion entre ambos coeficientes es

\[\hat\beta_K^{\,simple} \;=\; \hat\beta_3 \;+\; \hat\beta_2\,\hat\delta ,\]

donde \(\hat\delta\) es la pendiente de la regresión de \(\text{pH}\) sobre \(K\). El segundo término es el sesgo por variable omitida: mide cuánto del efecto del pH se añade adicionalmente en el coeficiente de \(K\) cuando el pH se deja fuera del modelo. Por tanto, los dos coeficientes solo coinciden si \(\hat\beta_2 = 0\), es decir, si el pH no afecta a la biomasa, o si \(\hat\delta = 0\), que significa que pH y \(K\) no están correlacionados (son ortogonales). En estos datos, cor(pH, K) = 0.092 y \(\hat\beta_2\) = 469.729, así que el sesgo vale 0.180 y el coeficiente simple (\(\hat\beta_K^{simple}\) = -0.345) difiere del parcial (\(\hat\beta_3\) = -0.525). La identidad se verifica numéricamente: \(-0.525 + 0.180 = -0.345\), con una diferencia de 0 entre ambos lados.

Es verdad que la correlacion entre pH y \(K\) es baja (0.092) pero aun asi el sesgo no es despreciable porque el pH tiene un efecto parcial muy grande sobre la biomasa (\(\hat\beta_2\) = 469.729): el sesgo depende del producto \(\hat\beta_2 \cdot \hat\delta\), no solo de la correlación. Adicionalmente a esto, como el sesgo es positivo, la regresión simple subestima en valor absoluto el efecto negativo de \(K\): pasa de -0.525 (efecto parcial, controlando por pH) a -0.345 (efecto marginal). El signo se mantiene, pero la magnitud del efecto de \(K\) es un 34 % menor en la regresión simple.

Esto se conecta con lo visto en el apartado c) ya que ahi vemos que \(\hat\beta_3\) es la pendiente entre las partes de la biomasa y de \(K\) no explicadas por el pH (teorema de Frisch-Waugh-Lovell), mientras que la regresión simple usa las variables completas.

3 Ejercicio 3. Selección de variables y esperanza de vida

3.1 a) EDA

esp <- read.table("esperanza.txt", header = TRUE)

cat("Dimensiones:", dim(esp), "\n")
## Dimensiones: 45 8
kable(head(esp), caption = "Primeras filas del dataset (esperanza de vida)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Primeras filas del dataset (esperanza de vida)
habitantes ingresos analfabetismo esp_vida asesinatos universitarios heladas area
25715 3804 2.1 69.05 15.1 41.3 66 50708
22466 6495 1.5 69.31 11.3 66.7 199 566432
24213 3558 1.9 70.66 10.1 39.9 114 51945
43302 5294 1.1 71.71 10.3 62.6 70 156361
24646 5064 0.7 72.06 6.8 63.9 217 103766
22686 4989 0.9 70.06 6.2 54.6 156 1982 
# Nombre de la variable respuesta: la columna que NO es un predictor del enunciado
predictores <- c("habitantes", "analfabetismo", "ingresos", "asesinatos",
                 "universitarios", "heladas", "area")
resp <- setdiff(names(esp), predictores)
resp        # debe ser un único nombre (la esperanza de vida)
## [1] "esp_vida"
# --- Datos faltantes ---
kable(t(colSums(is.na(esp))), caption = "Valores faltantes por variable") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Valores faltantes por variable
habitantes ingresos analfabetismo esp_vida asesinatos universitarios heladas area
0 0 0 0 0 0 0 0 
# --- Descriptivos ---
desc <- data.frame(
  media = sapply(esp, mean),
  sd    = sapply(esp, sd),
  cv    = sapply(esp, function(x) sd(x) / mean(x)),
  min   = sapply(esp, min),
  max   = sapply(esp, max)
)
kable(round(desc, 3), caption = "Estadísticos descriptivos: Esperanza de Vida") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Estadísticos descriptivos: Esperanza de Vida
media sd cv min max
habitantes 26431.600 4632.280 0.175 22466.00 43302.0
ingresos 4596.178 619.595 0.135 3278.00 6495.0
analfabetismo 1.178 0.626 0.531 0.50 2.8
esp_vida 70.820 1.356 0.019 67.96 73.6
asesinatos 7.607 3.656 0.481 1.40 15.1
universitarios 52.809 8.443 0.160 37.80 67.3
heladas 177.244 54.265 0.306 56.00 268.0
area 72315.911 89003.029 1.231 1049.00 566432.0 
library(ggcorrplot)

# --- Matriz de correlaciones (Gráfico) ---
ggcorrplot(cor(esp),
           hc.order = FALSE,
           type = "lower",
           lab = TRUE,
           lab_size = 3.5,
           method = "square",
           colors = c("#2166ac", "#f7f7f7", "#b2182b"),
           title = "Matriz de Correlación Completa: Esperanza de Vida y Predictores",
           ggtheme = theme_minimal())

# Correlación de cada predictor con la respuesta (graficado)
cor_resp <- sort(cor(esp[, predictores], esp[[resp]])[, 1], decreasing = TRUE)
df_cor_resp <- data.frame(
  predictor = factor(names(cor_resp), levels = names(cor_resp)),
  correlacion = as.numeric(cor_resp)
)

ggplot(df_cor_resp, aes(x = reorder(predictor, correlacion), y = correlacion, fill = correlacion)) +
  geom_col(width = 0.6) +
  geom_text(aes(label = round(correlacion, 3)), 
            hjust = ifelse(df_cor_resp$correlacion >= 0, -0.25, 1.25), size = 3.8, fontface = "bold") +
  coord_flip() +
  scale_fill_gradient2(low = "#2166ac", mid = "#f7f7f7", high = "#b2182b", midpoint = 0, limits = c(-1, 1)) +
  scale_y_continuous(limits = c(-1, 1)) +
  labs(title = "Correlación de cada Predictor con Esperanza de Vida",
       x = "Predictor", y = "Coeficiente de Correlación (r)") +
  theme_minimal() +
  theme(legend.position = "none")

library(ggplot2)
library(patchwork)

graf <- lapply(predictores, function(v) {
  ggplot(esp, aes(.data[[v]], .data[[resp]])) +
    geom_point(color = "steelblue", alpha = 0.7) +
    geom_smooth(method = "lm", color = "red") +
    labs(title = paste(resp, "vs.", v), x = v, y = resp)
})

wrap_plots(graf, ncol = 3)

# Relaciones entre los propios predictores
pairs(esp[, predictores], pch = 19, col = "steelblue", cex = 0.6)

# Correlaciones entre predictores (busca las de mayor valor absoluto)
cp <- cor(esp[, predictores])
cp[upper.tri(cp, diag = TRUE)] <- NA
pares_cor <- as.data.frame(as.table(cp))
pares_cor <- na.omit(pares_cor)
names(pares_cor) <- c("Variable 1", "Variable 2", "Correlación")
top5 <- head(pares_cor[order(-abs(pares_cor$Correlación)), ], 5)
top5$Correlación <- round(top5$Correlación, 5)
kable(top5,
      caption = "Top 5 pares de predictores con mayor correlación (en valor absoluto)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Top 5 pares de predictores con mayor correlación (en valor absoluto)
Variable 1 Variable 2 Correlación
11 asesinatos analfabetismo 0.70845
12 universitarios analfabetismo -0.67952
13 heladas analfabetismo -0.65479
19 universitarios ingresos 0.64165
27 heladas asesinatos -0.54625

Interpretación :

  • Estructura: 50 observaciones (los estados de EE. UU.) y 8 variables: la esperanza de vida y 7 predictores, todas numéricas y sin datos faltantes como lo podemos ver con colSum(.isna).
  • Predictores más asociados con la esperanza de vida: Aqui se ve claramente en la matriz que la variable con mayor r absoluto y por ende mas asociado es “asesinatos” con una r absoluta que casi tiene un valor de 0.8.
  • Predictores con asociación débil: En este apartado tenemos a habitantes y area, con una r absoluta menor que 0.1 respecto a la esperanza de vida. Estos dos son candidatos a quedar fuera de los modelos de selección.
  • Relaciones entre predictores: Para ver si hay multicolinealidad, vemos por parejas de predictores las correlaciones entre ellas. Viendo las que mayor r absoluto tienen entre ellas ya que hay muchas variables explicativas y no tendría sentido analizar cada una, veos que, por un lado, con signo positivo es decir, si en una vbariable se obtiene una valoracion alta en la otra también, tenemos a asesinatos y analfabetismo y universitarios y ingresos, las dos con una r mayor a 0.64. Estas asociaciones tienen sentido real ya que el estatus social y económico alto de la gente de un grupo se entiende como. POr otro lado, con signo negativo y significando que uhna valoración alta en un factor se asocia con un valor bajo en la otra, tenemos a universitarios y analfabetismo por un lado, que es una asociación inversa trivial, y por otro heladas y analfabetismo y heladas. Me parece curiosa está asociación. Creo que obviamente no tendrá una causalidad la una con la otra pero es curioso observarlo. Este analisis es importanete a la hora de hacer el analisis multivariado ya que aunque individualmente sean muy interesantes a la hora de crear un modelo, si dos predictores aportan información parecida, al incluir uno el otro puede dejar de ser significativo, y el orden de entrada en la selección dependerá de esa estructura.

3.2 b) Selección progresiva (forward) con add1()

Se parte del modelo nulo y, en cada etapa, se calcula con add1() el test F parcial de cada variable candidata. Ingresa la de menor p-valor, siempre que sea < 0.01; el proceso se detiene cuando ninguna de ellas obtiene un p–valor menor que 0.01.

alpha <- 0.01
# --- Etapa 0: modelo nulo y primer add1() ---
mod_actual <- lm(reformulate("1", resp), data = esp)
scope_f    <- reformulate(predictores)          # ~ todos los predictores

add1(mod_actual, scope = scope_f, test = "F")
## Single term additions
## 
## Model:
## esp_vida ~ 1
##                Df Sum of Sq    RSS     AIC F value    Pr(>F)    
## <none>                      80.852  28.368                      
## habitantes      1     0.365 80.487  30.164  0.1952   0.66082    
## analfabetismo   1    28.769 52.083  10.578 23.7521 1.530e-05 ***
## ingresos        1     7.820 73.032  25.791  4.6044   0.03758 *  
## asesinatos      1    50.882 29.970 -14.290 73.0022 8.152e-11 ***
## universitarios  1    28.097 52.755  11.155 22.9018 2.036e-05 ***
## heladas         1     4.823 76.029  27.601  2.7277   0.10591    
## area            1     0.775 80.077  29.935  0.4160   0.52239    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Para no repetir el proceso a mano, el siguiente chunk lo automatiza y muestra el add1() de cada etapa:

mod_actual <- lm(reformulate("1", resp), data = esp)
incluidas  <- character(0)
historial  <- data.frame(paso = integer(0), variable = character(0),
                         F = numeric(0), p_valor = numeric(0))

repeat {
  candidatas <- setdiff(predictores, incluidas)
  if (length(candidatas) == 0) break

  tab <- add1(mod_actual, scope = reformulate(predictores), test = "F")
  tab <- tab[candidatas, , drop = FALSE]              # solo las aún no incluidas

  cat("\n--- Paso", length(incluidas) + 1, ": modelo actual =",
      deparse(formula(mod_actual)), "---\n")
  print(tab[, c("Df", "Sum of Sq", "RSS", "AIC", "F value", "Pr(>F)")])

  mejor <- rownames(tab)[which.min(tab[["Pr(>F)"]])]
  p_min <- min(tab[["Pr(>F)"]])
  f_min <- tab[mejor, "F value"]

  if (p_min < alpha) {
    incluidas  <- c(incluidas, mejor)
    historial  <- rbind(historial,
                        data.frame(paso = length(incluidas), variable = mejor,
                                   F = f_min, p_valor = p_min))
    mod_actual <- update(mod_actual, reformulate(incluidas, resp))
    cat("=> Ingresa:", mejor, " (p =", format(p_min, digits = 4), ")\n")
  } else {
    cat("=> Ninguna variable cumple p <", alpha, ". Fin de la selección.\n")
    break
  }
}
## 
## --- Paso 1 : modelo actual = esp_vida ~ 1 ---
##                Df Sum of Sq    RSS     AIC F value    Pr(>F)    
## habitantes      1     0.365 80.487  30.164  0.1952   0.66082    
## analfabetismo   1    28.769 52.083  10.578 23.7521 1.530e-05 ***
## ingresos        1     7.820 73.032  25.791  4.6044   0.03758 *  
## asesinatos      1    50.882 29.970 -14.290 73.0022 8.152e-11 ***
## universitarios  1    28.097 52.755  11.155 22.9018 2.036e-05 ***
## heladas         1     4.823 76.029  27.601  2.7277   0.10591    
## area            1     0.775 80.077  29.935  0.4160   0.52239    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## => Ingresa: asesinatos  (p = 8.152e-11 )
## 
## --- Paso 2 : modelo actual = esp_vida ~ asesinatos ---
##                Df Sum of Sq    RSS     AIC F value  Pr(>F)  
## habitantes      1    3.2826 26.688 -17.510  5.1659 0.02821 *
## analfabetismo   1    0.1933 29.777 -12.581  0.2726 0.60435  
## ingresos        1    1.0796 28.891 -13.941  1.5694 0.21722  
## universitarios  1    4.2458 25.725 -19.165  6.9320 0.01180 *
## heladas         1    4.1207 25.850 -18.946  6.6953 0.01322 *
## area            1    0.4782 29.492 -13.014  0.6811 0.41388  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## => Ninguna variable cumple p < 0.01 . Fin de la selección.
mod_forward <- mod_actual
# --- Orden de ingreso y p-valores ---
kable(historial, digits = 6, row.names = FALSE,
      caption = "Selección progresiva (forward): orden de ingreso (alfa = 0.01)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Selección progresiva (forward): orden de ingreso (alfa = 0.01)
paso variable F p_valor
1 asesinatos 73.00217 0 
# --- Modelo final ---
cat("Fórmula del modelo final (forward):", deparse(formula(mod_forward)), "\n")
## Fórmula del modelo final (forward): esp_vida ~ asesinatos
coef_fw <- as.data.frame(summary(mod_forward)$coefficients)
names(coef_fw) <- c("Estimación", "Error Estándar", "t", "p-valor")
kable(round(coef_fw, 6), caption = "Coeficientes del modelo final (forward)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
  column_spec(5, bold = TRUE, color = ifelse(coef_fw$`p-valor` < 0.05, "green", "red"))
Coeficientes del modelo final (forward)
Estimación Error Estándar t p-valor
(Intercept) 73.057001 0.289915 251.994794 0
asesinatos -0.294114 0.034423 -8.544131 0 
gl_fw <- data.frame(`R²` = summary(mod_forward)$r.squared,
                    `R² ajustado` = summary(mod_forward)$adj.r.squared,
                    `Sigma` = summary(mod_forward)$sigma, check.names = FALSE)
kable(round(gl_fw, 6), caption = "Bondad de ajuste: modelo forward") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Bondad de ajuste: modelo forward
R² R² ajustado Sigma
0.629317 0.620697 0.834858

Interpretación:

  • Orden de ingreso: la única variable que entra al modelo es asesinatos, con un p-valor de

    8.152e-11

    Es la variable de mayor asociación marginal con la esperanza de vida (coincide con la de mayor |r| del EDA).

  • Detención: el proceso para cuando ninguna variable restante tiene p < 0.01; las variables excluidas son todas las demás.

  • Modelo final (forward): esp_vida~ asesinatos, con R² ajustado = 0.6207

  • Nota importante: el p-valor de cada variable depende de qué otras ya estén en el modelo; por eso un predictor correlacionado con otro ya incluido puede perder significancia. Es por eso que los p-valpores de las variables que al principio eran casi 0 suben a mas de 0.1

  • Nota importante 2: Si hubieramos sido mas plausibles con alpha y tomado un nivel de significancia mayor, con add1() añadiriamos al modelo la variable “universitarios”. Esto puede ser interesnate para considerar un modelo que quizá obtenga mejores resultados predictivos por ejemplo.

3.3 c) Selección regresiva (backward) con drop1()

Se parte del modelo completo y, en cada etapa, drop1() calcula el test F parcial de cada variable incluida (efecto de quitarla, dado que las demás permanecen). Se elimina la de mayor p-valor, siempre que sea > 0.01; el proceso se detiene cuando todas son significativas.

# --- Etapa 0: modelo completo y primer drop1() ---
mod_completo <- lm(reformulate(predictores, resp), data = esp)
summary(mod_completo)
## 
## Call:
## lm(formula = reformulate(predictores, resp), data = esp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.13551 -0.45779  0.00357  0.54006  1.37940 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     7.145e+01  2.030e+00  35.188  < 2e-16 ***
## habitantes      5.222e-05  2.738e-05   1.907   0.0643 .  
## analfabetismo  -3.699e-02  3.613e-01  -0.102   0.9190    
## ingresos       -2.527e-04  2.537e-04  -0.996   0.3257    
## asesinatos     -3.198e-01  4.689e-02  -6.819 4.94e-08 ***
## universitarios  5.351e-02  2.249e-02   2.380   0.0226 *  
## heladas        -6.931e-03  2.907e-03  -2.384   0.0224 *  
## area            4.635e-07  1.626e-06   0.285   0.7772    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7055 on 37 degrees of freedom
## Multiple R-squared:  0.7722, Adjusted R-squared:  0.7292 
## F-statistic: 17.92 on 7 and 37 DF,  p-value: 3.915e-10
drop1(mod_completo, test = "F")
## Single term deletions
## 
## Model:
## esp_vida ~ habitantes + analfabetismo + ingresos + asesinatos + 
##     universitarios + heladas + area
##                Df Sum of Sq    RSS     AIC F value   Pr(>F)    
## <none>                      18.414 -24.209                     
## habitantes      1    1.8099 20.224 -21.991  3.6366  0.06431 .  
## analfabetismo   1    0.0052 18.419 -26.197  0.0105  0.91901    
## ingresos        1    0.4937 18.908 -25.019  0.9919  0.32574    
## asesinatos      1   23.1439 41.558  10.419 46.5037 4.94e-08 ***
## universitarios  1    2.8181 21.232 -19.802  5.6624  0.02260 *  
## heladas         1    2.8283 21.242 -19.780  5.6831  0.02237 *  
## area            1    0.0404 18.455 -26.111  0.0813  0.77719    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El siguiente chunk automatiza el proceso y muestra el drop1() de cada etapa:

mod_actual <- lm(reformulate(predictores, resp), data = esp)
incluidas  <- predictores
historial  <- data.frame(paso = integer(0), variable_eliminada = character(0),
                         F = numeric(0), p_valor = numeric(0))

repeat {
  if (length(incluidas) == 0) break

  tab <- drop1(mod_actual, test = "F")
  tab <- tab[incluidas, , drop = FALSE]              # quita la fila <none>

  cat("\n--- Paso", nrow(historial) + 1, ": modelo actual =",
      deparse(formula(mod_actual)), "---\n")
  print(tab[, c("Df", "Sum of Sq", "RSS", "AIC", "F value", "Pr(>F)")])

  peor  <- rownames(tab)[which.max(tab[["Pr(>F)"]])]
  p_max <- max(tab[["Pr(>F)"]])
  f_max <- tab[peor, "F value"]

  if (p_max > alpha) {
    historial  <- rbind(historial,
                        data.frame(paso = nrow(historial) + 1,
                                   variable_eliminada = peor,
                                   F = f_max, p_valor = p_max))
    incluidas  <- setdiff(incluidas, peor)
    mod_actual <- if (length(incluidas) > 0) {
      lm(reformulate(incluidas, resp), data = esp)
    } else {
      lm(reformulate("1", resp), data = esp)
    }
    cat("=> Se elimina:", peor, " (p =", format(p_max, digits = 4), ")\n")
  } else {
    cat("=> Todas las variables son significativas (p <=", alpha,
        "). Fin de la selección.\n")
    break
  }
}
## 
## --- Paso 1 : modelo actual = esp_vida ~ habitantes + analfabetismo + ingresos + asesinatos +      universitarios + heladas + area ---
##                Df Sum of Sq    RSS     AIC F value   Pr(>F)    
## habitantes      1    1.8099 20.224 -21.991  3.6366  0.06431 .  
## analfabetismo   1    0.0052 18.419 -26.197  0.0105  0.91901    
## ingresos        1    0.4937 18.908 -25.019  0.9919  0.32574    
## asesinatos      1   23.1439 41.558  10.419 46.5037 4.94e-08 ***
## universitarios  1    2.8181 21.232 -19.802  5.6624  0.02260 *  
## heladas         1    2.8283 21.242 -19.780  5.6831  0.02237 *  
## area            1    0.0404 18.455 -26.111  0.0813  0.77719    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## => Se elimina: analfabetismo  (p = 0.919 )
## 
## --- Paso 2 : modelo actual = esp_vida ~ habitantes + ingresos + asesinatos + universitarios +      heladas + area ---
##                Df Sum of Sq    RSS     AIC F value    Pr(>F)    
## habitantes      1    1.9255 20.345 -23.723  3.9724  0.053469 .  
## ingresos        1    0.4931 18.912 -27.008  1.0174  0.319518    
## asesinatos      1   26.1458 44.565  11.563 53.9400 8.452e-09 ***
## universitarios  1    3.9038 22.323 -19.547  8.0538  0.007247 ** 
## heladas         1    3.8665 22.286 -19.622  7.9769  0.007505 ** 
## area            1    0.0352 18.455 -28.111  0.0727  0.788930    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## => Se elimina: area  (p = 0.7889 )
## 
## --- Paso 3 : modelo actual = esp_vida ~ habitantes + ingresos + asesinatos + universitarios +      heladas ---
##                Df Sum of Sq    RSS     AIC F value    Pr(>F)    
## habitantes      1     1.896 20.350 -25.711  4.0058  0.052335 .  
## ingresos        1     0.458 18.913 -29.008  0.9679  0.331277    
## asesinatos      1    33.535 51.990  16.497 70.8700 2.641e-10 ***
## universitarios  1     4.583 23.038 -20.128  9.6860  0.003468 ** 
## heladas         1     3.832 22.286 -21.621  8.0978  0.007028 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## => Se elimina: ingresos  (p = 0.3313 )
## 
## --- Paso 4 : modelo actual = esp_vida ~ habitantes + asesinatos + universitarios + heladas ---
##                Df Sum of Sq    RSS     AIC F value    Pr(>F)    
## habitantes      1     1.489 20.402 -27.597  3.1495  0.083563 .  
## asesinatos      1    33.746 52.658  15.072 71.3718 1.972e-10 ***
## universitarios  1     4.764 23.676 -20.898 10.0753  0.002889 ** 
## heladas         1     4.127 23.040 -22.125  8.7293  0.005226 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## => Se elimina: habitantes  (p = 0.08356 )
## 
## --- Paso 5 : modelo actual = esp_vida ~ asesinatos + universitarios + heladas ---
##                Df Sum of Sq    RSS     AIC F value    Pr(>F)    
## asesinatos      1    32.259 52.660  13.074  64.828 5.635e-10 ***
## universitarios  1     5.448 25.850 -18.946  10.949  0.001958 ** 
## heladas         1     5.323 25.725 -19.165  10.697  0.002179 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## => Todas las variables son significativas (p <= 0.01 ). Fin de la selección.
mod_backward <- mod_actual
# --- Orden de eliminación y p-valores ---
kable(historial, digits = 6, row.names = FALSE,
      caption = "Selección regresiva (backward): orden de eliminación (alfa = 0.01)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Selección regresiva (backward): orden de eliminación (alfa = 0.01)
paso variable_eliminada F p_valor
1 analfabetismo 0.010481 0.919010
2 area 0.072682 0.788930
3 ingresos 0.967869 0.331277
4 habitantes 3.149462 0.083563 
# --- Modelo final ---
cat("Fórmula del modelo final (backward):", deparse(formula(mod_backward)), "\n")
## Fórmula del modelo final (backward): esp_vida ~ asesinatos + universitarios + heladas
coef_bw <- as.data.frame(summary(mod_backward)$coefficients)
names(coef_bw) <- c("Estimación", "Error Estándar", "t", "p-valor")
kable(round(coef_bw, 6), caption = "Coeficientes del modelo final (backward)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
  column_spec(5, bold = TRUE, color = ifelse(coef_bw$`p-valor` < 0.05, "green", "red"))
Coeficientes del modelo final (backward)
Estimación Error Estándar t p-valor
(Intercept) 71.935735 1.026193 70.099582 0.000000
asesinatos -0.301891 0.037495 -8.051591 0.000000
universitarios 0.048241 0.014579 3.308882 0.001958
heladas -0.007713 0.002358 -3.270686 0.002179 
gl_bw <- data.frame(`R²` = summary(mod_backward)$r.squared,
                    `R² ajustado` = summary(mod_backward)$adj.r.squared,
                    `Sigma` = summary(mod_backward)$sigma, check.names = FALSE)
kable(round(gl_bw, 6), caption = "Bondad de ajuste: modelo backward") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Bondad de ajuste: modelo backward
R² R² ajustado Sigma
0.747667 0.729204 0.705409

Interpretacion:

  • Orden de eliminación: primero sale analfabetirsmo, luego area, luego ingresos y po rultimo habitantes.. Cada variable se elimina cuando, dadas las demás, su aporte no es significativo al 1 %.

    Vemos que el p-valor de analfbaetismo es altisimo. Esto es sin dudad por la multicolinealidad que hemos visto que tienen analfabetismo y universitarios. Ya que estan fuertemente e inversamente relacionados. Luego si el modelo primero contuviese analfabetismno y luego universitarios, yo creo que se qeudaría fuera universitarios y no al reves.

  • Detención: el proceso para cuando todas las variables restantes tienen p ≤ 0.01; las que quedan son asesinatos universitarios y heladas.

  • Modelo final (backward): resp ~ asesinatos + universitarios + heladas, con R² ajustado = 0.7292

  • Comparación con forward: El modelo NO coincide con el apartado b (del add1). Esto se debe a que backward evalúa cada variable en presencia de todas las demás, y forward solo en presencia de las ya incluidas; con predictores correlacionados (como ya hemos visto anteriormente que si sucede) los caminos pueden divergir como es el caso de estas variables y estos modelos.

3.4 d) Selección stepwise bidireccional

Se combinan las dos direcciones: tras cada ingreso (add1()), se revisa si las variables ya incluidas siguen siendo significativas (drop1()), y se elimina la que deje de cumplir el criterio del 1 %.

ajustar <- function(vars) {
  if (length(vars) == 0) lm(reformulate("1", resp), data = esp)
  else                   lm(reformulate(vars, resp), data = esp)
}

incluidas <- character(0)
mod_actual <- ajustar(incluidas)
trayectoria <- data.frame(paso = integer(0), accion = character(0),
                          variable = character(0), F = numeric(0),
                          p_valor = numeric(0), modelo = character(0),
                          stringsAsFactors = FALSE)
registrar <- function(accion, variable, F, p) {
  trayectoria <<- rbind(trayectoria,
    data.frame(paso = nrow(trayectoria) + 1, accion = accion, variable = variable,
               F = F, p_valor = p,
               modelo = paste(incluidas, collapse = " + "),
               stringsAsFactors = FALSE))
}

max_iter <- 4 * length(predictores)     # tope de seguridad contra ciclos
iter <- 0
repeat {
  iter <- iter + 1
  if (iter > max_iter) { cat("Se alcanzó el tope de iteraciones.\n"); break }

  # ---------- 1) PASO DE INGRESO (add1) ----------
  candidatas <- setdiff(predictores, incluidas)
  entro <- FALSE
  if (length(candidatas) > 0) {
    tab_a <- add1(mod_actual, scope = reformulate(predictores), test = "F")
    tab_a <- tab_a[candidatas, , drop = FALSE]
    cat("\n=== add1: modelo actual =", deparse(formula(mod_actual)), "===\n")
    print(tab_a[, c("Df", "Sum of Sq", "RSS", "AIC", "F value", "Pr(>F)")])

    mejor <- rownames(tab_a)[which.min(tab_a[["Pr(>F)"]])]
    p_min <- min(tab_a[["Pr(>F)"]])
    if (p_min < alpha) {
      incluidas  <- c(incluidas, mejor)
      mod_actual <- ajustar(incluidas)
      registrar("ingresa", mejor, tab_a[mejor, "F value"], p_min)
      cat("=> Ingresa:", mejor, "(p =", format(p_min, digits = 4), ")\n")
      entro <- TRUE
    }
  }

  # ---------- 2) PASO DE REVISIÓN (drop1) ----------
  salio <- FALSE
  repeat {
    if (length(incluidas) == 0) break
    tab_d <- drop1(mod_actual, test = "F")
    tab_d <- tab_d[incluidas, , drop = FALSE]
    cat("\n=== drop1: modelo actual =", deparse(formula(mod_actual)), "===\n")
    print(tab_d[, c("Df", "Sum of Sq", "RSS", "AIC", "F value", "Pr(>F)")])

    peor  <- rownames(tab_d)[which.max(tab_d[["Pr(>F)"]])]
    p_max <- max(tab_d[["Pr(>F)"]])
    if (p_max > alpha) {
      incluidas  <- setdiff(incluidas, peor)
      mod_actual <- ajustar(incluidas)
      registrar("sale", peor, tab_d[peor, "F value"], p_max)
      cat("=> Sale:", peor, "(p =", format(p_max, digits = 4), ")\n")
      salio <- TRUE
    } else break
  }

  # ---------- 3) CRITERIO DE PARADA ----------
  if (!entro && !salio) {
    cat("\n=> No es posible agregar ni eliminar variables. Fin de la selección.\n")
    break
  }
}
## 
## === add1: modelo actual = esp_vida ~ 1 ===
##                Df Sum of Sq    RSS     AIC F value    Pr(>F)    
## habitantes      1     0.365 80.487  30.164  0.1952   0.66082    
## analfabetismo   1    28.769 52.083  10.578 23.7521 1.530e-05 ***
## ingresos        1     7.820 73.032  25.791  4.6044   0.03758 *  
## asesinatos      1    50.882 29.970 -14.290 73.0022 8.152e-11 ***
## universitarios  1    28.097 52.755  11.155 22.9018 2.036e-05 ***
## heladas         1     4.823 76.029  27.601  2.7277   0.10591    
## area            1     0.775 80.077  29.935  0.4160   0.52239    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## => Ingresa: asesinatos (p = 8.152e-11 )
## 
## === drop1: modelo actual = esp_vida ~ asesinatos ===
##            Df Sum of Sq    RSS    AIC F value    Pr(>F)    
## asesinatos  1    50.882 80.852 28.368  73.002 8.152e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## === add1: modelo actual = esp_vida ~ asesinatos ===
##                Df Sum of Sq    RSS     AIC F value  Pr(>F)  
## habitantes      1    3.2826 26.688 -17.510  5.1659 0.02821 *
## analfabetismo   1    0.1933 29.777 -12.581  0.2726 0.60435  
## ingresos        1    1.0796 28.891 -13.941  1.5694 0.21722  
## universitarios  1    4.2458 25.725 -19.165  6.9320 0.01180 *
## heladas         1    4.1207 25.850 -18.946  6.6953 0.01322 *
## area            1    0.4782 29.492 -13.014  0.6811 0.41388  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## === drop1: modelo actual = esp_vida ~ asesinatos ===
##            Df Sum of Sq    RSS    AIC F value    Pr(>F)    
## asesinatos  1    50.882 80.852 28.368  73.002 8.152e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## => No es posible agregar ni eliminar variables. Fin de la selección.
mod_stepwise <- mod_actual
# --- Trayectoria completa ---
kable(trayectoria, digits = 6, row.names = FALSE,
      caption = "Stepwise bidireccional: trayectoria completa (alfa = 0.01)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Stepwise bidireccional: trayectoria completa (alfa = 0.01)
paso accion variable F p_valor modelo
1 ingresa asesinatos 73.00217 0 asesinatos 
# --- Modelo final ---
cat("Fórmula del modelo final (stepwise):", deparse(formula(mod_stepwise)), "\n")
## Fórmula del modelo final (stepwise): esp_vida ~ asesinatos
coef_sw <- as.data.frame(summary(mod_stepwise)$coefficients)
names(coef_sw) <- c("Estimación", "Error Estándar", "t", "p-valor")
kable(round(coef_sw, 6), caption = "Coeficientes del modelo final (stepwise)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
  column_spec(5, bold = TRUE, color = ifelse(coef_sw$`p-valor` < 0.05, "green", "red"))
Coeficientes del modelo final (stepwise)
Estimación Error Estándar t p-valor
(Intercept) 73.057001 0.289915 251.994794 0
asesinatos -0.294114 0.034423 -8.544131 0 
gl_sw <- data.frame(`R²` = summary(mod_stepwise)$r.squared,
                    `R² ajustado` = summary(mod_stepwise)$adj.r.squared,
                    `Sigma` = summary(mod_stepwise)$sigma, check.names = FALSE)
kable(round(gl_sw, 6), caption = "Bondad de ajuste: modelo stepwise") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Bondad de ajuste: modelo stepwise
R² R² ajustado Sigma
0.629317 0.620697 0.834858

Interpretacion:

  • Trayectoria: el proceso empieza con el modelo nulo; entra asesinatos con el mismo p-valor del forward selection (trivial) y en el siguiente paso no se pude añadir ninguna variable al modelo (p-vallores>0.01) y tampoco se pude eliminar ninguno, es decir asesinatos obviamente. Este modelo coincide con el forward. Al tener unicamente una variable en el forward, podiamos haber previsto que iba a ocurrir esto.
  • Modelo final (stepwise): resp ~ asesinatos, con R² ajustado = 0.6207 (mismo que forward obvio)
  • Comparación con forward y backward: Con el stepwise direccional, vemos que se obtiene el mismo modelo que con el forward y distinto al backward. Tendriamos que hacer distintas pruebsa para ver cual de los dos modelos es mejor en cada situación. V9iendo solo el R^2, podemos decir que el modelo del backward explica un 10% más de la variabilidad que el otro modelo. Veremos otros criterios en el siguiente apartado.

3.5 e) Comparación de modelos

Me voy a referir al modelo de solo asesinatos como el modelo forward y al de heladas, asesinatos y universiatrios como el modelo backward en este apartado.

predictores <- c("habitantes", "analfabetismo", "ingresos", "asesinatos",
                 "universitarios", "heladas", "area")
resp <- setdiff(names(esp), predictores)
stopifnot(length(resp) == 1)

# Modelo completo: su sigma^2 se usa como referencia para el Cp de Mallows
mod_completo <- lm(reformulate(predictores, resp), data = esp)
s2_completo  <- summary(mod_completo)$sigma^2
n <- nrow(esp)
# Indicadores de un modelo (p = número de parámetros, incluido el intercepto)
indicadores <- function(mod) {
  p   <- length(coef(mod))
  sse <- sum(resid(mod)^2)
  c(p     = p,
    R2_aj = summary(mod)$adj.r.squared,
    Cp    = sse / s2_completo - (n - 2 * p),
    AIC   = AIC(mod),
    BIC   = BIC(mod))
}

tabla_e <- rbind(Forward  = indicadores(mod_forward),
                 Backward = indicadores(mod_backward),
                 Stepwise = indicadores(mod_stepwise))
tabla_e <- as.data.frame(tabla_e)
kable(round(tabla_e, 4),
      caption = "Comparación de modelos: R²aj ↑ (max), Cp ↓, AIC ↓, BIC ↓ (min)") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
  column_spec(2, bold = TRUE)
Comparación de modelos: R²aj ↑ (max), Cp ↓, AIC ↓, BIC ↓ (min)
p R2_aj Cp AIC BIC
Forward 2 0.6207 19.2205 115.4142 120.8342
Backward 4 0.7292 3.9936 102.1074 111.1407
Stepwise 2 0.6207 19.2205 115.4142 120.8342 
# Variables de cada modelo
vars_df <- data.frame(
  Modelo    = c("Forward", "Backward", "Stepwise"),
  Variables = c(
    paste(attr(terms(mod_forward),  "term.labels"), collapse = ", "),
    paste(attr(terms(mod_backward), "term.labels"), collapse = ", "),
    paste(attr(terms(mod_stepwise), "term.labels"), collapse = ", ")
  )
)
kable(vars_df, caption = "Variables incluidas en cada modelo final") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
Variables incluidas en cada modelo final
Modelo Variables
Forward asesinatos
Backward asesinatos, universitarios, heladas
Stepwise asesinatos 
# Convención de add1/drop1/step: extractAIC() omite constantes
# extractAIC = n*log(SSE/n) + 2p ; AIC() = n*log(SSE/n) + n*(1+log(2*pi)) + 2(p+1)
tabla_conv <- t(sapply(list(Forward = mod_forward, Backward = mod_backward,
                            Stepwise = mod_stepwise),
                       function(m) c(`AIC()` = AIC(m),
                                     `extractAIC()` = extractAIC(m)[2])))
kable(round(tabla_conv, 4),
      caption = "AIC() frente a extractAIC(): difieren por una constante") |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
AIC() frente a extractAIC(): difieren por una constante
AIC() extractAIC()
Forward 115.4142 -14.2902
Backward 102.1074 -27.5971
Stepwise 115.4142 -14.2902

Selección del modelo (completa con tus resultados):

  • R² ajustado (se prefiere grande): mide la proporción de variabilidad explicada penalizando por el número de predictores; a diferencia del R², no aumenta automáticamente al añadir variables. El modelo con mayor R²aj es como ya he indicado anyteriormente el del Backward. Con una explicaicon de la variabilidad 10% mayor (72%)
  • Cp de Mallows (se prefiere pequeño y cercano a p): compara el error de cada submodelo con el del modelo completo; un Cp ≈ p indica poco sesgo por variables omitidas, y un Cp muy superior a p señala que faltan variables relevantes. El menor Cp lo tiene también el backward (4 frente a 19.22 del forward). Tiene sentido ya que se omiten menos varibles
  • AIC (se prefiere pequeño): equilibra ajuste y complejidad con penalización 2 por parámetro. El menor AIC lo tiene también el Backward. Aunque se utilicen más variables (y se penalice más), es significativamente mejor que el forward, por eso da un mejor resultado en el AIC también.
  • BIC (se prefiere pequeño): penaliza más la complejidad (log(n) por parámetro, con n = 50 ⇒ ≈ 3.9 por parámetro), por lo que tiende a modelos más parsimoniosos. El menor BIC lo tiene de nuevo el backward, con un 111 frente al 120 del forward. Misma exolicación que para AIC.
  • Modelo elegido: Como los 4 criterios coinciden y diocen que el backward da el mejor modelo (asesinatos + universitarios + heladas), nos quedamos con ese modelo. Si hubiera discrepancias, se prioriza la parsimonia y la interpretación, y se explica la discrepancia y porque puede haber sucedido eso.