Universidad Pontificia Bolivariana

Regresión lineal, semestre 2026-20

Modelo de regresión lineal múltiple para la resistencia a la compresión del concreto

Integrantes

Katerine Correa Lugo

Tomás Wrueta Lugo

Nicolás Sánchez Espitia

Alejandro Castillo Villadiego

Profesor

Delio Salgado Chamorro

Fecha de entrega

10 de noviembre de 2026

0. Introducción

0.1 Problema analizado

La resistencia a la compresión es la propiedad que más se utiliza para evaluar la calidad del concreto, porque de ella depende que una estructura soporte las cargas para las que fue diseñada. Su valor no se conoce en el momento de preparar la mezcla: se obtiene rompiendo probetas en el laboratorio después de un tiempo de curado, lo que hace lento y costoso el control de calidad.

Este trabajo ajusta un modelo de regresión lineal múltiple para explicar la resistencia a la compresión (csMPa, en MPa) a partir de la composición de la mezcla y de la edad de la probeta. Un modelo adecuado permitiría estimar la resistencia esperada de una mezcla antes de ensayarla y conocer qué componentes la aumentan y cuáles la disminuyen.

La regresión lineal múltiple es apropiada porque la variable de interés es continua, se dispone de varias variables explicativas cuantitativas y se quiere medir el aporte de cada una manteniendo constantes las demás. Los datos corresponden al conjunto Concrete Compressive Strength de Yeh (1998), disponible en Kaggle (ver referencias), en el archivo Concrete_Data_Yeh.csv. Cada fila es una mezcla de concreto ensayada a una edad determinada.

0.2 Variables del problema

La variable respuesta es la resistencia a la compresión y las ocho variables regresoras son los componentes de la mezcla y la edad.

Tabla 1. Variables del estudio
Variable Descripción Unidad Rol
csMPa Resistencia a la compresión MPa Respuesta (\(y\))
cement Cemento kg/m³ Regresora \(x_1\)
slag Escoria de alto horno kg/m³ Regresora \(x_2\)
flyash Ceniza volante kg/m³ Regresora \(x_3\)
water Agua kg/m³ Regresora \(x_4\)
superplasticizer Superplastificante kg/m³ Regresora \(x_5\)
coarseaggregate Agregado grueso kg/m³ Regresora \(x_6\)
fineaggregate Agregado fino kg/m³ Regresora \(x_7\)
age Edad de la muestra días Regresora \(x_8\)

Todas las variables son cuantitativas, por lo que no hay variables cualitativas que requieran tratamiento mediante variables indicadoras.

0.3 Análisis descriptivo

Primero se cargan los datos y se definen las variables. Después se revisan los estadísticos básicos, los valores faltantes y las correlaciones.

library(nortest)
library(lmtest)
library(car)
library(sandwich)

Concreto <- read.csv("Concrete_Data_Yeh.csv")

csMPa            <- Concreto$csMPa
cement           <- Concreto$cement
slag             <- Concreto$slag
flyash           <- Concreto$flyash
water            <- Concreto$water
superplasticizer <- Concreto$superplasticizer
coarseaggregate  <- Concreto$coarseaggregate
fineaggregate    <- Concreto$fineaggregate
age              <- Concreto$age

n <- length(csMPa)
k <- 8
summary(Concreto)
##      cement           slag           flyash           water      
##  Min.   :102.0   Min.   :  0.0   Min.   :  0.00   Min.   :121.8  
##  1st Qu.:192.4   1st Qu.:  0.0   1st Qu.:  0.00   1st Qu.:164.9  
##  Median :272.9   Median : 22.0   Median :  0.00   Median :185.0  
##  Mean   :281.2   Mean   : 73.9   Mean   : 54.19   Mean   :181.6  
##  3rd Qu.:350.0   3rd Qu.:142.9   3rd Qu.:118.30   3rd Qu.:192.0  
##  Max.   :540.0   Max.   :359.4   Max.   :200.10   Max.   :247.0  
##  superplasticizer coarseaggregate  fineaggregate        age        
##  Min.   : 0.000   Min.   : 801.0   Min.   :594.0   Min.   :  1.00  
##  1st Qu.: 0.000   1st Qu.: 932.0   1st Qu.:731.0   1st Qu.:  7.00  
##  Median : 6.400   Median : 968.0   Median :779.5   Median : 28.00  
##  Mean   : 6.205   Mean   : 972.9   Mean   :773.6   Mean   : 45.66  
##  3rd Qu.:10.200   3rd Qu.:1029.4   3rd Qu.:824.0   3rd Qu.: 56.00  
##  Max.   :32.200   Max.   :1145.0   Max.   :992.6   Max.   :365.00  
##      csMPa      
##  Min.   : 2.33  
##  1st Qu.:23.71  
##  Median :34.45  
##  Mean   :35.82  
##  3rd Qu.:46.13  
##  Max.   :82.60
colSums(is.na(Concreto))
##           cement             slag           flyash            water 
##                0                0                0                0 
## superplasticizer  coarseaggregate    fineaggregate              age 
##                0                0                0                0 
##            csMPa 
##                0
round(cor(Concreto), 3)
##                  cement   slag flyash  water superplasticizer coarseaggregate
## cement            1.000 -0.275 -0.397 -0.082            0.092          -0.109
## slag             -0.275  1.000 -0.324  0.107            0.043          -0.284
## flyash           -0.397 -0.324  1.000 -0.257            0.378          -0.010
## water            -0.082  0.107 -0.257  1.000           -0.658          -0.182
## superplasticizer  0.092  0.043  0.378 -0.658            1.000          -0.266
## coarseaggregate  -0.109 -0.284 -0.010 -0.182           -0.266           1.000
## fineaggregate    -0.223 -0.282  0.079 -0.451            0.223          -0.178
## age               0.082 -0.044 -0.154  0.278           -0.193          -0.003
## csMPa             0.498  0.135 -0.106 -0.290            0.366          -0.165
##                  fineaggregate    age  csMPa
## cement                  -0.223  0.082  0.498
## slag                    -0.282 -0.044  0.135
## flyash                   0.079 -0.154 -0.106
## water                   -0.451  0.278 -0.290
## superplasticizer         0.223 -0.193  0.366
## coarseaggregate         -0.178 -0.003 -0.165
## fineaggregate            1.000 -0.156 -0.167
## age                     -0.156  1.000  0.329
## csMPa                   -0.167  0.329  1.000
pairs(Concreto, pch = 20, cex = 0.4)

Principales comportamientos. La base tiene 1030 observaciones y 9 variables, con 0 valores faltantes. La resistencia tiene una media de 35.82 MPa y una desviación estándar de 16.71 MPa, con valores entre 2.33 y 82.6 MPa, de modo que hay mezclas de resistencia muy distinta. Las variables slag, flyash y superplasticizer tienen muchos valores iguales a cero (45.7 %, 55 % y 36.8 % de las mezclas, respectivamente). Esos ceros no son datos faltantes: indican que la mezcla no llevó ese componente. La edad tiene una distribución muy asimétrica: su media es 45.7 días, su mediana 28 y su máximo 365, porque la mayoría de los ensayos se hicieron a edades tempranas y solo unos pocos a edades largas. En la matriz de dispersión se observa además que age toma pocos valores distintos, ya que los ensayos se programan a edades estándar.

Entre las regresoras, la correlación lineal más alta con la resistencia es la de cement (0.498) y la más baja es la de water (-0.29). Estas correlaciones son solo una primera aproximación, porque no controlan el efecto de las demás variables; ese control lo hace el modelo de regresión.

1. Modelo y estimación de parámetros

1.1 Explicación del procedimiento

El modelo de regresión lineal múltiple expresa la resistencia como una combinación lineal de las ocho regresoras más un error aleatorio. Cada coeficiente \(\beta_j\) mide el cambio promedio de la resistencia cuando la variable \(x_j\) aumenta una unidad y las demás permanecen constantes. Los parámetros se estiman por mínimos cuadrados ordinarios (MCO), que escogen los valores que hacen mínima la suma de los cuadrados de los residuales.

1.2 Fórmulas

El modelo general, para \(i = 1, \dots, n\), es

\[csMPa_i = \beta_0 + \beta_1\, cement_i + \beta_2\, slag_i + \beta_3\, flyash_i + \beta_4\, water_i + \beta_5\, superplasticizer_i + \beta_6\, coarseaggregate_i + \beta_7\, fineaggregate_i + \beta_8\, age_i + \varepsilon_i, \qquad \varepsilon_i \sim N(0, \sigma^2)\]

En forma matricial, \(\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\), y el estimador de MCO es

\[\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y}.\]

1.3 Código y resultados

modelo <- lm(csMPa ~ cement + slag + flyash + water + superplasticizer +
               coarseaggregate + fineaggregate + age)
summary(modelo)
## 
## Call:
## lm(formula = csMPa ~ cement + slag + flyash + water + superplasticizer + 
##     coarseaggregate + fineaggregate + age)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.654  -6.302   0.703   6.569  34.450 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -23.331214  26.585504  -0.878 0.380372    
## cement             0.119804   0.008489  14.113  < 2e-16 ***
## slag               0.103866   0.010136  10.247  < 2e-16 ***
## flyash             0.087934   0.012583   6.988 5.02e-12 ***
## water             -0.149918   0.040177  -3.731 0.000201 ***
## superplasticizer   0.292225   0.093424   3.128 0.001810 ** 
## coarseaggregate    0.018086   0.009392   1.926 0.054425 .  
## fineaggregate      0.020190   0.010702   1.887 0.059491 .  
## age                0.114222   0.005427  21.046  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.4 on 1021 degrees of freedom
## Multiple R-squared:  0.6155, Adjusted R-squared:  0.6125 
## F-statistic: 204.3 on 8 and 1021 DF,  p-value: < 2.2e-16

La ecuación estimada es

\[\widehat{csMPa} = -23.3312 +0.1198\,cement +0.1039\,slag +0.0879\,flyash -0.1499\,water +0.2922\,superplasticizer +0.0181\,coarseaggregate +0.0202\,fineaggregate +0.1142\,age\]

1.4 Interpretación

Los coeficientes de cement, slag, flyash, superplasticizer, coarseaggregate, fineaggregate, age son positivos, de modo que aumentar esas cantidades (o la edad) eleva la resistencia estimada. Los coeficientes de water son negativos, de modo que aumentar esas cantidades la reduce. En todos los casos el efecto corresponde a un aumento de 1 kg/m³ en los componentes o de 1 día en la edad, con las demás variables constantes.

El intercepto estimado es -23.33 MPa y corresponde a una mezcla con todas las cantidades en cero. Esa situación no existe en los datos (por ejemplo, el agua mínima observada es de 121.8 kg/m³), por lo que el intercepto sirve para ajustar el modelo y no tiene interpretación práctica. Además, los componentes de una mezcla no varían de forma independiente en la realidad, de modo que la expresión “manteniendo constantes las demás” debe leerse con cuidado.

El modelo tiene un \(R^2\) de 0.6155, un \(R^2\) ajustado de 0.6125 y un error estándar residual de 10.4 MPa.

1.5 Conclusión del apartado

Se estimó el modelo completo con las ocho regresoras, que explica el 61.6 % de la variabilidad de la resistencia del concreto. Esta ecuación es el modelo base sobre el cual se desarrollan las pruebas y los diagnósticos de los apartados siguientes.

2. Inferencia sobre el modelo completo

2.1 Explicación del procedimiento

La prueba F global responde a una primera pregunta: ¿el modelo, en conjunto, explica la resistencia mejor que un modelo sin regresoras, que solo usaría la media? Para resolverla se separa la variabilidad total de la resistencia en dos partes: la que explica la regresión y la que queda en los residuales. Esa información se resume en la tabla de análisis de varianza (ANOVA). Se trabaja con un nivel de significancia \(\alpha = 0.05\).

2.2 Hipótesis

\[H_0: \beta_1 = \beta_2 = \cdots = \beta_8 = 0 \qquad \text{vs.} \qquad H_1: \text{al menos un } \beta_j \neq 0, \; j = 1, \dots, 8\]

2.3 Fórmulas

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

\[MSR = \frac{SSR}{k}, \qquad MSE = \frac{SSE}{n - k - 1}, \qquad F_0 = \frac{MSR}{MSE}\]

Se rechaza \(H_0\) si el estadístico de prueba supera al estadístico teórico, es decir, si \(F_0 > F_{\alpha;\,k,\,n-k-1}\).

2.4 Código y resultados

SST <- sum((csMPa - mean(csMPa))^2)
SSE <- sum(modelo$residuals^2)
SSR <- SST - SSE

gl_reg <- k
gl_res <- n - k - 1
gl_tot <- n - 1

CMR <- SSR / gl_reg
CME <- SSE / gl_res

F0    <- CMR / CME
F_teo <- qf(0.95, gl_reg, gl_res)
p_F   <- pf(F0, gl_reg, gl_res, lower.tail = FALSE)

tabla_anova <- data.frame(
  Fuente = c("Regresión", "Residual", "Total"),
  Suma_cuadrados = c(SSR, SSE, SST),
  gl = c(gl_reg, gl_res, gl_tot),
  Cuadrado_medio = c(CMR, CME, NA),
  F0 = c(F0, NA, NA))
knitr::kable(tabla_anova, digits = 4, na = "",
             caption = "Tabla 2. Análisis de varianza del modelo completo")
Tabla 2. Análisis de varianza del modelo completo
Fuente Suma_cuadrados gl Cuadrado_medio F0
Regresión 176762.0 8 22095.2542 204.3167
Residual 110413.2 1021 108.1422 NA
Total 287175.2 1029 NA NA 
c(F_prueba = F0, F_teorico = F_teo, p_valor = p_F)
##      F_prueba     F_teorico       p_valor 
##  2.043167e+02  1.947456e+00 6.285978e-206

2.5 Comparación con el criterio de decisión

El estadístico de prueba es \(F_0 =\) 204.32 y el estadístico teórico es \(F_{0.05;\,8,\,1021} =\) 1.9475. Como el estadístico de prueba F0 es mayor que el estadístico teórico y el p-valor es < 1e-04, menor que 0.05, se rechaza H0.

2.6 Interpretación y conclusión

Con un nivel de significancia del 5 %, hay evidencia de que al menos una de las variables, ya sea un componente de la mezcla o la edad, está relacionada con la resistencia del concreto. Esto significa que el modelo completo explica la resistencia mejor que su promedio.

3. Inferencias para cada uno de los coeficientes

3.1 Explicación del procedimiento

Una vez establecido que el modelo es significativo en conjunto, interesa saber qué variables aportan de manera individual. La prueba t evalúa, para cada coeficiente, si su valor es distinto de cero cuando las demás variables ya están incluidas en el modelo. Se incluye también el intercepto y se usa \(\alpha = 0.05\).

3.2 Hipótesis

\[H_0: \beta_j = 0 \qquad \text{vs.} \qquad H_1: \beta_j \neq 0, \qquad j = 0, 1, \dots, 8\]

3.3 Fórmulas

\[\text{ee}(\hat{\beta}_j) = \sqrt{MSE\,[(\mathbf{X}^\top\mathbf{X})^{-1}]_{jj}}, \qquad t_0 = \frac{\hat{\beta}_j}{\text{ee}(\hat{\beta}_j)}\]

Se rechaza \(H_0\) si el estadístico de prueba supera al estadístico teórico, es decir, si \(|t_0| > t_{\alpha/2;\,n-k-1}\).

3.4 Código y resultados

tab_coef <- summary(modelo)$coefficients
gl       <- n - k - 1
t_crit   <- qt(0.975, gl)

inf_t <- data.frame(
  Coeficiente = rownames(tab_coef),
  Estimado    = tab_coef[, "Estimate"],
  Error_est   = tab_coef[, "Std. Error"],
  t0          = tab_coef[, "t value"],
  t_teorico   = t_crit,
  p_valor     = format.pval(tab_coef[, "Pr(>|t|)"], digits = 4, eps = 0.0001),
  Decision    = ifelse(abs(tab_coef[, "t value"]) > t_crit,
                       "Rechaza H0", "No rechaza H0"),
  row.names   = NULL)
knitr::kable(inf_t, digits = 4,
             caption = "Tabla 3. Pruebas t individuales sobre los coeficientes")
Tabla 3. Pruebas t individuales sobre los coeficientes
Coeficiente Estimado Error_est t0 t_teorico p_valor Decision
(Intercept) -23.3312 26.5855 -0.8776 1.9623 0.3803719 No rechaza H0
cement 0.1198 0.0085 14.1129 1.9623 < 1e-04 Rechaza H0
slag 0.1039 0.0101 10.2474 1.9623 < 1e-04 Rechaza H0
flyash 0.0879 0.0126 6.9882 1.9623 < 1e-04 Rechaza H0
water -0.1499 0.0402 -3.7314 1.9623 0.0002009 Rechaza H0
superplasticizer 0.2922 0.0934 3.1279 1.9623 0.0018102 Rechaza H0
coarseaggregate 0.0181 0.0094 1.9257 1.9623 0.0544248 No rechaza H0
fineaggregate 0.0202 0.0107 1.8867 1.9623 0.0594908 No rechaza H0
age 0.1142 0.0054 21.0464 1.9623 < 1e-04 Rechaza H0

3.5 Comparación con el criterio de decisión

El estadístico teórico es \(t_{0.025;\,1021} =\) 1.9623. Se rechaza \(H_0\) cuando el estadístico de prueba \(|t_0|\) supera al estadístico teórico, lo que equivale a un p-valor menor que 0.05. Con este criterio se rechaza \(H_0\) para 6 de los 8 coeficientes de las variables regresoras.

3.6 Interpretación y conclusión

Con las demás variables en el modelo, tienen un efecto significativo y positivo sobre la resistencia cement, slag, flyash, superplasticizer, age, de modo que un aumento en estas cantidades incrementa la resistencia media del concreto. Tienen un efecto significativo y negativo water, de modo que un aumento en estas cantidades la reduce. La magnitud de cada efecto, en MPa por unidad, es la de la columna Estimado. Las variables coarseaggregate, fineaggregate no resultan significativas. Esto no quiere decir que no influyan en la resistencia, sino que, con las demás variables en el modelo, su aporte adicional no se distingue de cero.

También conviene distinguir entre significancia estadística e importancia práctica. Con 1030 observaciones, incluso efectos pequeños pueden resultar significativos, por lo que el tamaño del coeficiente debe leerse junto con el p-valor y con el rango real de cada variable.

En conclusión, con un nivel de significancia del 5 %, 6 de las 8 variables regresoras son significativas en el modelo completo. Las que no lo son se vuelven a evaluar en la selección de variables del apartado 8.

4. Verificación del supuesto de normalidad

4.1 Explicación del supuesto

Las pruebas t y F de los apartados anteriores se construyen suponiendo que los errores del modelo siguen una distribución normal, \(\varepsilon_i \sim N(0, \sigma^2)\). Como los errores no se observan, el supuesto se verifica con los residuales, \(e_i = y_i - \hat{y}_i\). Primero se examina la forma de su distribución con gráficos y después se contrasta con tres pruebas formales: Lilliefors, Anderson-Darling y Shapiro-Wilk.

4.2 Análisis gráfico

residuales <- modelo$residuals
qqnorm(residuales, main = "Gráfico Q-Q de los residuales")
qqline(residuales, col = "red")

hist(residuales, breaks = 30, col = "green",
     main = "Histograma de residuales", xlab = "Residuales")

En el gráfico Q-Q se comparan los cuantiles de los residuales con los de una normal: si el supuesto se cumple, los puntos se ubican sobre la línea roja. Las desviaciones respecto a la línea, sobre todo en los extremos, indican que las colas de la distribución se apartan de la normal. En el histograma se espera una forma simétrica y de campana.

4.3 Pruebas formales

Hipótesis.

\[H_0: \varepsilon_i \sim \text{Normal} \qquad \text{vs.} \qquad H_1: \varepsilon_i \not\sim \text{Normal}\]

con \(\alpha = 0.05\).

Estadísticos. Con los residuales ordenados \(e_{(1)} \le \cdots \le e_{(n)}\) y \(z_{(i)} = (e_{(i)} - \bar{e})/s_e\), las tres pruebas usan los siguientes estadísticos.

Lilliefors (Kolmogorov-Smirnov con media y varianza estimadas):

\[D = \max_i \left\{ \frac{i}{n} - \Phi(z_{(i)}),\; \Phi(z_{(i)}) - \frac{i-1}{n} \right\}\]

Con \(n > 30\) y \(\alpha = 0.05\), el estadístico teórico aproximado es \(D_{teo} = 0.886/\sqrt{n}\) (Lilliefors, 1967), y se rechaza \(H_0\) si el estadístico de prueba supera al teórico, es decir, si \(D > D_{teo}\).

Anderson-Darling (Anderson y Darling, 1952):

\[A^2 = -n - \frac{1}{n}\sum_{i=1}^{n}(2i-1)\left[\ln \Phi(z_{(i)}) + \ln\left(1 - \Phi(z_{(n+1-i)})\right)\right]\]

Shapiro-Wilk (Shapiro y Wilk, 1965):

\[W = \frac{\left(\sum_{i=1}^{n} a_i\, e_{(i)}\right)^2}{\sum_{i=1}^{n}(e_i - \bar{e})^2}\]

Para Anderson-Darling y Shapiro-Wilk no hay un único estadístico teórico de tabla, por lo que la decisión se toma con el p-valor: se rechaza \(H_0\) si es menor que \(\alpha\). Las tres pruebas se aplican con el paquete nortest (Gross y Ligges, 2015).

Código y resultados.

lt <- lillie.test(residuales)
lt
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuales
## D = 0.033815, p-value = 0.007526
ad <- ad.test(residuales)
ad
## 
##  Anderson-Darling normality test
## 
## data:  residuales
## A = 1.2398, p-value = 0.003149
sw <- shapiro.test(residuales)
sw
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.99532, p-value = 0.002993

El estadístico de prueba de Lilliefors y su estadístico teórico se calculan de la siguiente manera:

ord   <- sort(residuales)
z     <- (ord - mean(residuales)) / sd(residuales)
Fz    <- pnorm(z)
i     <- 1:n
D     <- max(pmax(i / n - Fz, Fz - (i - 1) / n))
Dcrit <- 0.886 / sqrt(n)

data.frame(D_manual = D, D_lillie_test = unname(lt$statistic),
           D_teorico = Dcrit, p_valor = lt$p.value)
##     D_manual D_lillie_test  D_teorico     p_valor
## 1 0.03381455    0.03381455 0.02760674 0.007525946

4.4 Comparación con el criterio de decisión

Para Lilliefors, el estadístico de prueba es \(D =\) 0.0338 y el estadístico teórico es 0.0276. Como el estadístico de prueba D es mayor que el estadístico teórico, se rechaza H0. Los p-valores de las pruebas de Lilliefors, Anderson-Darling y Shapiro-Wilk son 0.007526, 0.003149 y 0.002993, respectivamente. Los tres son menores que 0.05, por lo que se rechaza H0 en las tres pruebas.

4.5 Interpretación, conclusión y recomendación

Las pruebas indican que los residuales no siguen una distribución normal; es decir, los errores que comete el modelo al predecir la resistencia del concreto no se distribuyen de forma normal. Por esta razón, los p-valores y la significancia de los apartados 2 y 3 deben tomarse con cautela. Como el supuesto no se cumple, se recomienda transformar alguna variable, en particular la edad, que es muy asimétrica; revisar la presencia de observaciones atípicas o influyentes; y emplear errores estándar robustos para la inferencia. Estas tres acciones se aplican en el apartado 7, con una transformación logarítmica de la edad, la distancia de Cook y los errores estándar HC3.

5. Verificación del supuesto de varianza constante

5.1 Explicación del supuesto

El modelo supone que los errores tienen la misma varianza en todas las observaciones, \(\text{Var}(\varepsilon_i) = \sigma^2\), lo que se conoce como homocedasticidad. Si la varianza cambia, las estimaciones de los coeficientes siguen siendo insesgadas, pero los errores estándar y las pruebas dejan de ser confiables. El supuesto se revisa con un gráfico de los residuales contra los valores ajustados y con la prueba de Breusch-Pagan (Breusch y Pagan, 1979), en su versión studentizada y en la original.

5.2 Análisis gráfico

ajustados <- modelo$fitted.values
plot(ajustados, residuales, pch = 20, main = "Residuales vs valores ajustados",
     xlab = "Valores ajustados", ylab = "Residuales")
abline(h = 0, col = "red")

Si la varianza es constante, los puntos forman una banda de ancho parejo alrededor de cero, sin forma de embudo ni de curva.

5.3 Prueba formal: Breusch-Pagan

Hipótesis.

\[H_0: \text{Var}(\varepsilon_i) = \sigma^2 \text{ para todo } i \qquad \text{vs.} \qquad H_1: \text{la varianza depende de las regresoras}\]

con \(\alpha = 0.05\).

Estadístico. Se ajusta una regresión auxiliar de \(e_i^2\) sobre las regresoras y se calcula

\[LM = n\,R^2_{aux} \;\overset{H_0}{\sim}\; \chi^2_{k}\]

La versión original supone errores normales; la studentizada (Koenker, 1981) no necesita ese supuesto. Se rechaza \(H_0\) si el estadístico de prueba supera al estadístico teórico, es decir, si \(LM > \chi^2_{\alpha;\,k}\).

Código y resultados.

bp_stud <- bptest(modelo)                   # Breusch-Pagan studentizado
bp_stud
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 137.2, df = 8, p-value < 2.2e-16
bp_orig <- bptest(modelo, studentize = F)   # Breusch-Pagan original
bp_orig
## 
##  Breusch-Pagan test
## 
## data:  modelo
## BP = 140.27, df = 8, p-value < 2.2e-16
qchisq(0.05, k, lower.tail = F)             # estadístico teórico
## [1] 15.50731

5.4 Comparación con el criterio de decisión

El estadístico teórico es \(\chi^2_{0.05;\,8} =\) 15.5073. La prueba studentizada da \(LM =\) 137.1963 (p-valor < 1e-04) y la original da \(LM =\) 140.2684 (p-valor < 1e-04). En ambas versiones el estadístico de prueba supera al estadístico teórico y el p-valor es menor que 0.05, por lo que se rechaza H0.

5.5 Interpretación, conclusión y recomendación

Hay evidencia de que la varianza de los errores no es constante: la dispersión de los errores de predicción cambia según las características de la mezcla, de modo que el modelo no predice con la misma precisión para todas las resistencias. En consecuencia, los errores estándar de MCO no son confiables y las pruebas t y F pueden conducir a conclusiones engañosas. Se recomienda transformar la variable respuesta o las regresoras para estabilizar la varianza, usar errores estándar robustos a heterocedasticidad (HC3) para la inferencia y revisar si la forma funcional del modelo está bien especificada, porque una relación no lineal mal modelada también produce este patrón. En el apartado 7 se aplican la transformación de la edad y los errores estándar HC3.

6. Verificación del supuesto de independencia

6.1 Explicación del supuesto

El modelo supone que los errores no están correlacionados entre sí. Este supuesto suele revisarse cuando las observaciones tienen un orden natural, por ejemplo en el tiempo. En esta base, el orden de las filas corresponde a cómo está organizado el conjunto de mezclas y no a una secuencia temporal, por lo que el resultado debe leerse con cautela. La verificación se hace con un gráfico de los residuales según el orden de las observaciones y con la prueba de Durbin-Watson (Durbin y Watson, 1950).

6.2 Análisis gráfico

orden <- 1:n
plot(orden, residuales, pch = 20, main = "Residuales según el orden de las observaciones",
     xlab = "Orden de las observaciones", ylab = "Residuales")
abline(h = 0, col = "red")

Si los errores son independientes, los residuales no muestran patrón alguno: ni rachas largas del mismo signo ni tendencias.

6.3 Prueba formal: Durbin-Watson

Hipótesis.

\[H_0: \rho = 0 \qquad \text{vs.} \qquad H_1: \rho \neq 0\]

donde \(\rho\) es la autocorrelación de primer orden de los errores, \(\varepsilon_t = \rho\,\varepsilon_{t-1} + u_t\), con \(\alpha = 0.05\).

Estadístico.

\[d = \frac{\sum_{t=2}^{n}(e_t - e_{t-1})^2}{\sum_{t=1}^{n} e_t^2} \approx 2\,(1 - \hat{\rho})\]

Un valor de \(d\) cercano a 2 indica ausencia de autocorrelación. Las tablas de Durbin-Watson no cubren \(n = 1030\), así que el estadístico teórico se obtiene con la aproximación normal: \(d\) tiene media 2 y varianza aproximada \(4/n\), por lo que se rechaza \(H_0\) si \(d\) queda fuera del intervalo

\[2 \pm z_{\alpha/2}\,\frac{2}{\sqrt{n}}\]

Código y resultados.

dw <- durbinWatsonTest(modelo, alternative = "two.sided")
dw
##  lag Autocorrelation D-W Statistic p-value
##    1       0.3559936      1.281638       0
##  Alternative hypothesis: rho != 0
# Las tablas de Durbin-Watson no cubren n = 1030, por eso se usa el p-valor
# La base esta ordenada por mezclas (no por tiempo): interpretar con cautela

El estadístico y los límites de no rechazo se calculan así:

d      <- sum(diff(residuales)^2) / sum(residuales^2)
r1     <- sum(residuales[-1] * residuales[-n]) / sum(residuales^2)
z_crit <- qnorm(0.975)
d_inf  <- 2 - z_crit * 2 / sqrt(n)
d_sup  <- 2 + z_crit * 2 / sqrt(n)

data.frame(d = d, r1 = r1, limite_inferior = d_inf, limite_superior = d_sup)
##          d        r1 limite_inferior limite_superior
## 1 1.281638 0.3559936         1.87786         2.12214

6.4 Comparación con el criterio de decisión

El estadístico de prueba es \(d =\) 1.2816, con una autocorrelación estimada de 0.356, y la región de no rechazo es el intervalo (1.8779; 2.1221). Como d cae fuera del intervalo, se rechaza H0. El p-valor de durbinWatsonTest es < 1e-04.

6.5 Interpretación, conclusión y recomendación

Hay evidencia de autocorrelación entre los residuales. Sin embargo, como las filas no están ordenadas en el tiempo, este resultado puede reflejar que la base está agrupada por tipos de mezcla, con mezclas parecidas una junto a otra, y no una dependencia temporal real. En conclusión, el supuesto de independencia no se cumple según la prueba de Durbin-Watson. Si el orden de los ensayos tuviera un significado real, como la fecha o el lote de producción, se recomienda repetir la prueba con ese orden y, de confirmarse la dependencia, ajustar el modelo teniéndola en cuenta. Si no lo tiene, se recomienda interpretar la prueba con cautela y usar errores estándar robustos, como se hace en el apartado 7.

7. Diagnóstico de la multicolinealidad

7.1 Explicación del procedimiento

Hay multicolinealidad cuando algunas regresoras están fuertemente relacionadas entre sí. En ese caso los coeficientes se vuelven inestables, sus errores estándar aumentan y resulta difícil separar el efecto de cada variable. En una mezcla de concreto es un problema posible, porque los componentes comparten el mismo volumen y se ajustan unos en función de otros. Para diagnosticarla se calcula el factor de inflación de la varianza (VIF) de cada regresora, con la función vif del paquete car (Fox y Weisberg, 2019).

7.2 Fórmula

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

donde \(R_j^2\) es el coeficiente de determinación de la regresión de \(x_j\) sobre las demás regresoras. Como criterio, un \(VIF < 5\) indica que no hay problema, un \(VIF\) entre 5 y 10 indica multicolinealidad moderada y un \(VIF > 10\) indica multicolinealidad severa.

7.3 Código y resultados

vif(modelo)
##           cement             slag           flyash            water 
##         7.488944         7.276963         6.170634         7.003957 
## superplasticizer  coarseaggregate    fineaggregate              age 
##         2.963776         5.074617         7.005081         1.118367

7.4 Comparación con el criterio de decisión

El mayor VIF es el de cement (7.489). Ninguna variable supera el valor de 10, pero 6 variables tienen un VIF de 5 o más (cement, slag, flyash, water, coarseaggregate, fineaggregate), lo que indica multicolinealidad moderada.

7.5 Interpretación, conclusión y recomendación

Algunos componentes de la mezcla están relacionados entre sí, de modo que sus coeficientes pueden ser inestables y es más difícil separar el efecto individual de cada ingrediente sobre la resistencia. Frente a esto se recomienda eliminar una de las variables muy correlacionadas o combinarlas, centrar las variables cuando hay términos relacionados entre sí y revisar la selección de variables del apartado 8, que puede reducir el problema.

7.6 Observaciones influyentes

Una observación influyente es una mezcla que, al retirarla, cambia de forma apreciable los coeficientes estimados. La distancia de Cook (Cook, 1977) mide esa influencia, y suele señalarse como influyente toda observación con \(D_i > 4/n\).

\[D_i = \frac{r_i^2}{k + 1}\cdot\frac{h_{ii}}{1 - h_{ii}}\]

cook <- cooks.distance(modelo)
plot(cook, type = "h", main = "Distancia de Cook",
     ylab = "Distancia de Cook", xlab = "Observación")
abline(h = 4/n, col = "red")

sum(cook > 4/n)
## [1] 79
sum(abs(rstandard(modelo)) > 3)
## [1] 2

Hay 79 observaciones con distancia de Cook mayor que \(4/n =\) 0.00388 y 2 con residual estandarizado de valor absoluto mayor que 3.

7.7 Modelo mejorado: logaritmo de la edad

La edad influye sobre la resistencia de forma no lineal: el concreto gana resistencia con rapidez en los primeros días y después más despacio. Por eso se prueba un modelo con el logaritmo de la edad y su cuadrado. En modelo2 el término cuadrático se calcula sin centrar. En modelo3 el logaritmo de la edad se centra en su media antes de elevarlo al cuadrado, lo que reduce la multicolinealidad entre ambos términos.

\[logage_c = \ln(age) - \overline{\ln(age)}, \qquad logage\_c2 = logage_c^{\,2}\]

logage   <- log(age)
logage_c <- logage - mean(logage)
logage_c2 <- logage_c^2

modelo2 <- lm(csMPa ~ cement + slag + flyash + water + superplasticizer +
                coarseaggregate + fineaggregate + logage + I(logage^2))
modelo3 <- lm(csMPa ~ cement + slag + flyash + water + superplasticizer +
                coarseaggregate + fineaggregate + logage_c + logage_c2)
vif(modelo2)
##           cement             slag           flyash            water 
##         7.582262         7.281209         6.141904         7.042623 
## superplasticizer  coarseaggregate    fineaggregate           logage 
##         2.984403         5.112441         7.053790        18.895152 
##      I(logage^2) 
##        19.092992
vif(modelo3)
##           cement             slag           flyash            water 
##         7.582262         7.281209         6.141904         7.042623 
## superplasticizer  coarseaggregate    fineaggregate         logage_c 
##         2.984403         5.112441         7.053790         1.072679 
##        logage_c2 
##         1.108075

Los supuestos de los apartados 4, 5 y 6 se verifican de nuevo sobre modelo3 con la función supuestos.

# Funcion que repite los supuestos para un modelo
supuestos <- function(m, nombre) {
  r <- m$residuals
  qqnorm(r, main = paste("Gráfico Q-Q", nombre)); qqline(r, col = "red")
  plot(m$fitted.values, r, pch = 20, xlab = "Valores ajustados",
       ylab = "Residuales", main = paste("Residuales vs. valores ajustados -", nombre))
  abline(h = 0, col = "red")
  print(lillie.test(r)); print(ad.test(r)); print(shapiro.test(r))
  print(bptest(m)); print(bptest(m, studentize = F))
  print(durbinWatsonTest(m, alternative = "two.sided"))
}

res2 <- modelo3$residuals
hist(res2, breaks = 30, col = "green", main = "Residuales del modelo mejorado")

supuestos(modelo3, "Modelo mejorado")

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  r
## D = 0.038886, p-value = 0.0008948
## 
## 
##  Anderson-Darling normality test
## 
## data:  r
## A = 2.128, p-value = 2.082e-05
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  r
## W = 0.99254, p-value = 4.682e-05
## 
## 
##  studentized Breusch-Pagan test
## 
## data:  m
## BP = 86.64, df = 9, p-value = 7.671e-15
## 
## 
##  Breusch-Pagan test
## 
## data:  m
## BP = 101.39, df = 9, p-value < 2.2e-16
## 
##  lag Autocorrelation D-W Statistic p-value
##    1       0.3540996      1.284039       0
##  Alternative hypothesis: rho != 0

Para la inferencia cuando algún supuesto no se cumple, se calculan los errores estándar robustos HC3, con los paquetes lmtest y sandwich (Zeileis, 2004; Zeileis y Hothorn, 2002), y se comparan los tres modelos con AIC y BIC.

# Errores estandar robustos (HC3)
coeftest(modelo3, vcov. = vcovHC(modelo3, type = "HC3"))
## 
## t test of coefficients:
## 
##                     Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)      -55.6890571  22.2674858 -2.5009 0.0125432 *  
## cement             0.1372653   0.0068356 20.0810 < 2.2e-16 ***
## slag               0.1165318   0.0086018 13.5474 < 2.2e-16 ***
## flyash             0.0938971   0.0094738  9.9112 < 2.2e-16 ***
## water             -0.1160249   0.0342693 -3.3857 0.0007371 ***
## superplasticizer   0.1067381   0.0752897  1.4177 0.1565841    
## coarseaggregate    0.0319016   0.0078138  4.0827 4.799e-05 ***
## fineaggregate      0.0381406   0.0087169  4.3755 1.336e-05 ***
## logage_c           8.6115779   0.1966768 43.7854 < 2.2e-16 ***
## logage_c2         -0.6524412   0.1420640 -4.5926 4.924e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Comparacion de modelos
AIC(modelo, modelo2, modelo3)
##         df      AIC
## modelo  10 7757.924
## modelo2 11 6960.565
## modelo3 11 6960.565
BIC(modelo, modelo2, modelo3)
##         df      BIC
## modelo  10 7807.297
## modelo2 11 7014.875
## modelo3 11 7014.875

Con el logaritmo de la edad centrado, el mayor VIF del modelo es 7.58, frente a 7.49 en el modelo original. Los p-valores de las pruebas de normalidad (Lilliefors, Anderson-Darling y Shapiro-Wilk) son 0.0008948, < 1e-04 y < 1e-04; los de Breusch-Pagan (studentizado y original) son < 1e-04 y < 1e-04, y el de Durbin-Watson es < 1e-04. El AIC pasa de 7757.9 en el modelo original a 6960.6 en modelo3, y el BIC pasa de 7807.3 a 7014.9; el modelo mejorado ajusta mejor según ambos criterios. Para los supuestos que sigan sin cumplirse, la tabla con errores estándar HC3 permite hacer la inferencia sobre los coeficientes sin depender de la varianza constante.

En conclusión, la transformación logarítmica de la edad con centrado es la recomendación que se aplica frente a los supuestos y a la multicolinealidad, y modelo3 es la base para la selección de variables del apartado 8.

8. Selección del modelo

8.1 Explicación del procedimiento

Seleccionar un modelo consiste en encontrar el subconjunto de variables que explique bien la resistencia sin incluir variables que aporten poco. Para ello se calculan cinco criterios sobre todos los subconjuntos posibles de las nueve variables del modelo mejorado (logage_c, logage_c2 y las siete variables de la mezcla), que son \(2^9 - 1 = 511\) modelos sin contar el que solo tiene el intercepto. Se mantiene la jerarquía: un modelo con logage_c2 solo se considera si también incluye logage_c. Después se aplican los algoritmos Forward Selection, Backward Elimination y Stepwise.

8.2 Fórmulas

Para un modelo con \(p\) parámetros (incluido el intercepto), donde \(\hat{\sigma}^2\) es el cuadrado medio residual del modelo completo y \(h_{ii}\) es el apalancamiento de la observación \(i\):

\[R^2_{aj} = 1 - \frac{SSE_p/(n - p)}{SST/(n - 1)}, \qquad C_p = \frac{SSE_p}{\hat{\sigma}^2} - (n - 2p)\]

\[AIC = -2\ln(L) + 2\,m, \qquad BIC = -2\ln(L) + m\,\ln(n), \qquad PRESS = \sum_{i=1}^{n}\left(\frac{e_i}{1 - h_{ii}}\right)^2\]

donde \(m\) es el número de parámetros que estima lm (los coeficientes y la varianza). Un mejor modelo tiene mayor \(R^2_{aj}\), menor AIC, BIC y PRESS, y un \(C_p\) pequeño y cercano a \(p\). El BIC penaliza la complejidad más que el AIC. El AIC corresponde a Akaike (1974) y el PRESS a Allen (1974).

8.3 Código y resultados: todas las regresiones posibles

datos2 <- data.frame(csMPa, logage_c, logage_c2, cement, slag, flyash, water,
                     superplasticizer, coarseaggregate, fineaggregate)
vars   <- names(datos2)[-1]
full   <- lm(csMPa ~ ., data = datos2)
sigma2 <- summary(full)$sigma^2
press  <- function(m) sum((residuals(m) / (1 - hatvalues(m)))^2)

tabla <- NULL
for (r in 1:length(vars)) {
  for (cmb in combn(vars, r, simplify = FALSE)) {
    m <- lm(reformulate(cmb, "csMPa"), data = datos2)
    p <- length(coef(m))
    tabla <- rbind(tabla, data.frame(
      variables = paste(cmb, collapse = " + "), p = p,
      R2aj = summary(m)$adj.r.squared,
      Cp   = sum(residuals(m)^2) / sigma2 - (n - 2 * p),
      AIC = AIC(m), BIC = BIC(m), PRESS = press(m)))
  }
}
nrow(tabla)                                   # 511 modelos
## [1] 511
# Jerarquia: si esta logage_c2 debe estar logage_c
tabla <- tabla[vapply(strsplit(tabla$variables, " + ", fixed = TRUE),
                      function(x) !("logage_c2" %in% x) || ("logage_c" %in% x),
                      logical(1)), ]
tabla <- tabla[order(tabla$BIC), ]

Los diez mejores modelos según el BIC y el mejor modelo según cada criterio son los siguientes.

knitr::kable(head(tabla, 10), digits = 3, row.names = FALSE,
             caption = "Tabla 4. Los diez mejores modelos según el BIC")
Tabla 4. Los diez mejores modelos según el BIC
variables p R2aj Cp AIC BIC PRESS
logage_c + logage_c2 + cement + slag + flyash + water + coarseaggregate + fineaggregate 9 0.821 10.814 6961.402 7010.776 51999.81
logage_c + logage_c2 + cement + slag + flyash + water + superplasticizer + coarseaggregate + fineaggregate 10 0.821 10.000 6960.565 7014.875 51995.36
logage_c + logage_c2 + cement + slag + flyash + water 7 0.817 32.935 6983.377 7022.876 53069.18
logage_c + logage_c2 + cement + slag + flyash + water + fineaggregate 8 0.818 30.860 6981.368 7025.804 52967.39
logage_c + logage_c2 + cement + slag + flyash + superplasticizer + coarseaggregate + fineaggregate 9 0.819 26.004 6976.587 7025.960 52760.88
logage_c + logage_c2 + cement + slag + flyash + water + coarseaggregate 8 0.817 34.198 6984.653 7029.088 53138.02
logage_c + logage_c2 + cement + slag + flyash + water + superplasticizer 8 0.817 34.699 6985.145 7029.581 53213.83
logage_c + cement + slag + flyash + water + coarseaggregate + fineaggregate 8 0.817 35.725 6986.152 7030.588 53255.94
logage_c + logage_c2 + cement + slag + flyash + water + superplasticizer + fineaggregate 9 0.817 32.859 6983.366 7032.739 53118.67
logage_c + cement + slag + flyash + water + superplasticizer + coarseaggregate + fineaggregate 9 0.817 34.233 6984.720 7034.093 53224.58 
tabla[which.max(tabla$R2aj), ]
##                                                                                                      variables
## 511 logage_c + logage_c2 + cement + slag + flyash + water + superplasticizer + coarseaggregate + fineaggregate
##      p     R2aj Cp      AIC      BIC    PRESS
## 511 10 0.821498 10 6960.565 7014.875 51995.36
tabla[which.min(tabla$Cp), ]
##                                                                                                      variables
## 511 logage_c + logage_c2 + cement + slag + flyash + water + superplasticizer + coarseaggregate + fineaggregate
##      p     R2aj Cp      AIC      BIC    PRESS
## 511 10 0.821498 10 6960.565 7014.875 51995.36
tabla[which.min(tabla$AIC), ]
##                                                                                                      variables
## 511 logage_c + logage_c2 + cement + slag + flyash + water + superplasticizer + coarseaggregate + fineaggregate
##      p     R2aj Cp      AIC      BIC    PRESS
## 511 10 0.821498 10 6960.565 7014.875 51995.36
tabla[which.min(tabla$BIC), ]
##                                                                                   variables
## 504 logage_c + logage_c2 + cement + slag + flyash + water + coarseaggregate + fineaggregate
##     p      R2aj       Cp      AIC      BIC    PRESS
## 504 9 0.8211808 10.81403 6961.402 7010.776 51999.81
tabla[which.min(tabla$PRESS), ]
##                                                                                                      variables
## 511 logage_c + logage_c2 + cement + slag + flyash + water + superplasticizer + coarseaggregate + fineaggregate
##      p     R2aj Cp      AIC      BIC    PRESS
## 511 10 0.821498 10 6960.565 7014.875 51995.36

8.4 Código y resultados: Forward Selection, Backward Elimination y Stepwise

En Forward Selection se parte del modelo sin variables y en cada paso entra la que más mejora el criterio. En Backward Elimination se parte del modelo completo y en cada paso sale la que menos aporta. Stepwise combina ambos movimientos. Forward y Backward usan el AIC; Stepwise usa el BIC.

nulo <- lm(csMPa ~ 1, data = datos2)
modelo_forward  <- step(nulo, scope = formula(full), direction = "forward", k = 2, trace = 0)
modelo_backward <- step(full, direction = "backward", k = 2, trace = 0)
modelo_stepwise <- step(nulo, scope = formula(full), direction = "both", k = log(n), trace = 0)

modelos <- list(modelo_forward, modelo_backward, modelo_stepwise)
data.frame(Metodo = c("Forward AIC", "Backward AIC", "Stepwise BIC"),
           R2_ajustado = sapply(modelos, function(m) summary(m)$adj.r.squared),
           AIC = sapply(modelos, AIC),
           BIC = sapply(modelos, BIC))
##         Metodo R2_ajustado      AIC      BIC
## 1  Forward AIC   0.8214980 6960.565 7014.875
## 2 Backward AIC   0.8214980 6960.565 7014.875
## 3 Stepwise BIC   0.8169725 6983.377 7022.876

Las variables que conserva cada método son:

attr(terms(modelo_forward), "term.labels")
## [1] "logage_c"         "cement"           "water"            "slag"            
## [5] "flyash"           "logage_c2"        "fineaggregate"    "coarseaggregate" 
## [9] "superplasticizer"
attr(terms(modelo_backward), "term.labels")
## [1] "logage_c"         "logage_c2"        "cement"           "slag"            
## [5] "flyash"           "water"            "superplasticizer" "coarseaggregate" 
## [9] "fineaggregate"
attr(terms(modelo_stepwise), "term.labels")
## [1] "logage_c"  "cement"    "water"     "slag"      "flyash"    "logage_c2"

8.5 Comparación con el criterio de decisión

Tras aplicar la jerarquía quedan 383 de los 511 modelos. Cada criterio elige el siguiente modelo:

  • \(R^2_{aj}\) (0.8215): logage_c + logage_c2 + cement + slag + flyash + water + superplasticizer + coarseaggregate + fineaggregate.
  • \(C_p\) (10): logage_c + logage_c2 + cement + slag + flyash + water + superplasticizer + coarseaggregate + fineaggregate.
  • AIC (6960.6): logage_c + logage_c2 + cement + slag + flyash + water + superplasticizer + coarseaggregate + fineaggregate.
  • BIC (7010.8): logage_c + logage_c2 + cement + slag + flyash + water + coarseaggregate + fineaggregate.
  • PRESS (5.19954^{4}): logage_c + logage_c2 + cement + slag + flyash + water + superplasticizer + coarseaggregate + fineaggregate.

Forward y Backward (con criterio AIC) conservan 9 y 9 variables, respectivamente, y Stepwise con criterio BIC conserva 6.

8.6 Modelo final

Como modelo final se toma el que resulta del método Stepwise con criterio BIC. Se revisan su resumen, el VIF, los supuestos y los errores estándar robustos HC3, y se evalúa su capacidad de predicción con validación cruzada de 10 particiones, comparándola con la del modelo original.

modelo_final <- modelo_stepwise
summary(modelo_final)
## 
## Call:
## lm(formula = csMPa ~ logage_c + cement + water + slag + flyash + 
##     logage_c2, data = datos2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.1594  -4.6004  -0.1999   4.6005  28.8100 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 40.384164   2.529418  15.966  < 2e-16 ***
## logage_c     8.531173   0.191980  44.438  < 2e-16 ***
## cement       0.110951   0.002692  41.210  < 2e-16 ***
## water       -0.244556   0.011340 -21.565  < 2e-16 ***
## slag         0.085111   0.003114  27.328  < 2e-16 ***
## flyash       0.059442   0.004640  12.811  < 2e-16 ***
## logage_c2   -0.612878   0.128412  -4.773 2.08e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.147 on 1023 degrees of freedom
## Multiple R-squared:  0.818,  Adjusted R-squared:  0.817 
## F-statistic: 766.5 on 6 and 1023 DF,  p-value: < 2.2e-16
vif(modelo_final)
##  logage_c    cement     water      slag    flyash logage_c2 
##  1.053974  1.594790  1.181378  1.454517  1.776174  1.098162
supuestos(modelo_final, "Modelo final")

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  r
## D = 0.03625, p-value = 0.002832
## 
## 
##  Anderson-Darling normality test
## 
## data:  r
## A = 1.0567, p-value = 0.008892
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  r
## W = 0.99617, p-value = 0.01217
## 
## 
##  studentized Breusch-Pagan test
## 
## data:  m
## BP = 71.904, df = 6, p-value = 1.662e-13
## 
## 
##  Breusch-Pagan test
## 
## data:  m
## BP = 81.439, df = 6, p-value = 1.801e-15
## 
##  lag Autocorrelation D-W Statistic p-value
##    1       0.3555746      1.282598       0
##  Alternative hypothesis: rho != 0
coeftest(modelo_final, vcov. = vcovHC(modelo_final, type = "HC3"))
## 
## t test of coefficients:
## 
##               Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 40.3841638  2.9595418  13.6454 < 2.2e-16 ***
## logage_c     8.5311728  0.1998716  42.6833 < 2.2e-16 ***
## cement       0.1109513  0.0025716  43.1449 < 2.2e-16 ***
## water       -0.2445560  0.0140400 -17.4185 < 2.2e-16 ***
## slag         0.0851107  0.0033594  25.3347 < 2.2e-16 ***
## flyash       0.0594421  0.0043769  13.5808 < 2.2e-16 ***
## logage_c2   -0.6128781  0.1406411  -4.3577 1.447e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
set.seed(123)
folds <- sample(rep(1:10, length.out = n))

cv_rmse <- function(fm, d) {
  mean(sapply(1:10, function(i) {
    m <- lm(fm, data = d[folds != i, ])
    sqrt(mean((d$csMPa[folds == i] - predict(m, d[folds == i, ]))^2))
  }))
}

RMSE_final    <- cv_rmse(formula(modelo_final), datos2)
RMSE_original <- cv_rmse(csMPa ~ cement + slag + flyash + water + superplasticizer +
                           coarseaggregate + fineaggregate + age, Concreto)

data.frame(Modelo = c("Modelo original", "Modelo final"),
           RMSE = c(RMSE_original, RMSE_final))
##            Modelo      RMSE
## 1 Modelo original 10.496535
## 2    Modelo final  7.171034

8.7 Interpretación y conclusión

El BIC penaliza más las variables adicionales que el AIC, por lo que tiende a elegir modelos más pequeños. El modelo que elige el BIC tiene 8 variables regresoras y el que elige el \(R^2_{aj}\) tiene 9. Las variables adicionales del modelo del R² ajustado aportan poco a la explicación de la resistencia, de modo que el modelo más pequeño conserva casi la misma capacidad explicativa con menos componentes.

El modelo final contiene las variables logage_c, cement, water, slag, flyash, logage_c2, con un \(R^2\) ajustado de 0.817 y un VIF máximo de 1.78. En la validación cruzada, el error de predicción (RMSE) baja de 10.5 MPa en el modelo original a 7.17 MPa en el modelo final. Por lo tanto, el modelo final es mejor para predecir la resistencia del concreto que el modelo con las ocho variables originales.

Referencias

Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716-723.

Allen, D. M. (1974). The relationship between variable selection and data augmentation and a method for prediction. Technometrics, 16(1), 125-127.

Anderson, T. W., & Darling, D. A. (1952). Asymptotic theory of certain “goodness of fit” criteria based on stochastic processes. The Annals of Mathematical Statistics, 23(2), 193-212.

Breusch, T. S., & Pagan, A. R. (1979). A simple test for heteroscedasticity and random coefficient variation. Econometrica, 47(5), 1287-1294.

Cook, R. D. (1977). Detection of influential observation in linear regression. Technometrics, 19(1), 15-18.

Durbin, J., & Watson, G. S. (1950). Testing for serial correlation in least squares regression I. Biometrika, 37(3-4), 409-428.

Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3.ª ed.). Sage.

Gross, J., & Ligges, U. (2015). nortest: Tests for normality. R package. https://CRAN.R-project.org/package=nortest

Kaggle. (s. f.). Concrete compressive strength [Conjunto de datos]. https://www.kaggle.com/

Koenker, R. (1981). A note on studentizing a test for heteroscedasticity. Journal of Econometrics, 17(1), 107-112.

Lilliefors, H. W. (1967). On the Kolmogorov-Smirnov test for normality with mean and variance unknown. Journal of the American Statistical Association, 62(318), 399-402.

Shapiro, S. S., & Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3-4), 591-611.

Yeh, I.-C. (1998). Modeling of strength of high-performance concrete using artificial neural networks. Cement and Concrete Research, 28(12), 1797-1808.

Zeileis, A. (2004). Econometric computing with HC and HAC covariance matrix estimators. Journal of Statistical Software, 11(10), 1-17.

Zeileis, A., & Hothorn, T. (2002). Diagnostic checking in regression relationships. R News, 2(3), 7-10.