Cómo usar este archivo. Compile con Knit en RStudio o con rmarkdown::render("taller_rlm.Rmd"). Para la versión del estudiante (sin soluciones): rmarkdown::render("taller_rlm.Rmd", params = list(solucionario = FALSE), output_file = "taller_rlm_estudiante.html"). Paquetes requeridos: car y lmtest (install.packages(c("car", "lmtest"))). El resto usa R base (≥ 4.1).

library(car)      # vif(), qqPlot(), influencePlot()
library(lmtest)   # bptest(), dwtest()

1 Presentación

Este taller desarrolla, en una sesión de 4 horas, el ajuste, la inferencia y el diagnóstico de modelos de regresión lineal múltiple (RLM) aplicados a datos químicos, usando exclusivamente R.

Objetivos de aprendizaje. Al finalizar, el estudiante será capaz de:

Formular un modelo de RLM en notación escalar y matricial, e identificar sus supuestos.
Estimar los coeficientes por mínimos cuadrados ordinarios (MCO), con álgebra matricial en R y con lm().
Construir e interpretar la tabla ANOVA, la prueba F global, las pruebas t individuales, \(R^2\) y \(R^2_{adj}\).
Calcular intervalos de confianza para los coeficientes y para la respuesta media, e intervalos de predicción.
Diagnosticar el modelo: normalidad y homocedasticidad de residuos, multicolinealidad (VIF) y observaciones influyentes.
Linealizar modelos químicos no lineales y reconocer los límites de la extrapolación.

Organización sugerida.

Bloque	Contenido	Tiempo
1	Fundamento teórico (secciones 2 a 4)	60 min
2	Ejemplos resueltos 1 y 2	60 min
3	Problemas 1 a 3 (en parejas)	75 min
4	Problemas 4 y 5, discusión y cierre	45 min

Convenciones. \(\alpha = 0.05\); \(n\) = observaciones; \(k\) = predictores; \(p = k + 1\) = parámetros.

2 Fundamento teórico

2.1 El modelo

\[ y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_k x_{ik} + \varepsilon_i, \qquad i = 1,\dots,n \]

“Lineal” se refiere a los parámetros, no a las variables: un modelo con \(x^2\), \(1/T\) o \(\ln x\) sigue siendo lineal (Montgomery et al., 2012, cap. 3).

Parámetros y estimadores. \(\beta_j\) es el parámetro poblacional (desconocido); \(b_j\) (o \(\hat\beta_j\)) es su estimación a partir de la muestra. \(b_0\) es el intercepto y cada \(b_j\) (\(j \ge 1\)) es un coeficiente parcial: el cambio estimado en \(y\) por unidad de \(x_j\) manteniendo constantes los demás predictores.

2.2 Notación matricial y mínimos cuadrados

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \mathbf{X}=\begin{pmatrix}1 & x_{11} & \cdots & x_{1k}\\ \vdots & \vdots & & \vdots\\ 1 & x_{n1} & \cdots & x_{nk}\end{pmatrix} \]

MCO minimiza \(S(\boldsymbol\beta) = (\mathbf y-\mathbf X\boldsymbol\beta)^{\mathsf T}(\mathbf y-\mathbf X\boldsymbol\beta)\), lo que conduce a las ecuaciones normales:

\[ \mathbf{X}^{\mathsf T}\mathbf{X}\,\mathbf{b} = \mathbf{X}^{\mathsf T}\mathbf{y} \;\Longrightarrow\; \mathbf{b} = (\mathbf{X}^{\mathsf T}\mathbf{X})^{-1}\mathbf{X}^{\mathsf T}\mathbf{y}, \qquad \hat{\mathbf y} = \mathbf H\mathbf y,\quad \mathbf H = \mathbf X(\mathbf X^{\mathsf T}\mathbf X)^{-1}\mathbf X^{\mathsf T} \]

\[ s^2 = MS_E = \frac{SS_E}{n-k-1}, \qquad \widehat{\mathrm{Var}}(\mathbf b) = s^2(\mathbf X^{\mathsf T}\mathbf X)^{-1}, \qquad se(b_j) = s\sqrt{C_{jj}} \]

2.3 Supuestos

Linealidad: \(E(\varepsilon)=0\); el modelo está bien especificado.
Homocedasticidad: \(\mathrm{Var}(\varepsilon_i)=\sigma^2\) constante.
Independencia de los errores.
Normalidad: \(\boldsymbol\varepsilon \sim N(\mathbf 0, \sigma^2\mathbf I)\) (necesaria para pruebas e intervalos exactos).

Con 1 a 3, \(\mathbf b\) es el mejor estimador lineal insesgado (teorema de Gauss-Markov; Kutner et al., 2005, cap. 6).

Demostración en R: MCO “a mano” frente a `lm()`

set.seed(1)
x1 <- runif(10, 0, 10); x2 <- runif(10, 0, 5)
y  <- 2 + 1.5 * x1 - 0.8 * x2 + rnorm(10, sd = 0.5)
X  <- cbind(1, x1, x2)                       # matriz de diseño
b  <- solve(t(X) %*% X, t(X) %*% y)          # ecuaciones normales
cbind(matricial = drop(b), lm = coef(lm(y ~ x1 + x2)))

#>    matricial       lm
#>      2.58112  2.58112
#> x1   1.39705  1.39705
#> x2  -0.76023 -0.76023

3 Inferencia

3.1 Tabla ANOVA y prueba F global

\[ SS_T = SS_R + SS_E,\qquad F_0 = \frac{SS_R/k}{SS_E/(n-k-1)} \sim F_{k,\,n-k-1}\ \text{bajo } H_0:\beta_1=\cdots=\beta_k=0 \]

Fuente	gl	SC	CM	\(F_0\)
Regresión	\(k\)	\(SS_R\)	\(SS_R/k\)	\(MS_R/MS_E\)
Residual	\(n-k-1\)	\(SS_E\)	\(SS_E/(n-k-1)\)
Total	\(n-1\)	\(SS_T\)

En R, anova(modelo) da la descomposición secuencial (tipo I) por predictor. La prueba F global aparece al final de summary(modelo) o con anova(lm(y ~ 1), modelo).

3.2 Pruebas t, R² y R² ajustado

\[ t_0 = \frac{b_j-\beta_{j,0}}{se(b_j)} \sim t_{n-k-1}, \qquad R^2 = 1-\frac{SS_E}{SS_T}, \qquad R^2_{adj} = 1-\frac{SS_E/(n-k-1)}{SS_T/(n-1)} \]

La prueba t es parcial: evalúa \(x_j\) dado que los demás ya están en el modelo. \(R^2\) nunca baja al añadir predictores; para comparar modelos use \(R^2_{adj}\), \(s\) y AIC.

3.3 Intervalos

\[ b_j \pm t_{\alpha/2,\,n-p}\,se(b_j);\qquad \hat y_0 \pm t_{\alpha/2,\,n-p}\sqrt{s^2\,\mathbf x_0^{\mathsf T}(\mathbf X^{\mathsf T}\mathbf X)^{-1}\mathbf x_0}\ \text{(media)};\qquad \hat y_0 \pm t_{\alpha/2,\,n-p}\sqrt{s^2\left(1+\mathbf x_0^{\mathsf T}(\mathbf X^{\mathsf T}\mathbf X)^{-1}\mathbf x_0\right)}\ \text{(predicción)} \]

En R: confint(m), predict(m, nuevo, interval = "confidence") y predict(m, nuevo, interval = "prediction").

4 Diagnóstico

Verificación	Gráfico en R	Prueba en R	Señal de problema
Linealidad	`plot(m, which = 1)`; `car::residualPlots(m)`	—	Curvatura
Homocedasticidad	`plot(m, which = 3)`	`lmtest::bptest(m)` (Breusch y Pagan, 1979)	Embudo; p < 0.05
Normalidad	`car::qqPlot(m)`	`shapiro.test(rstudent(m))` (Shapiro y Wilk, 1965)	p < 0.05
Independencia	`plot(residuals(m), type = "b")`	`lmtest::dwtest(m)` (Durbin y Watson, 1950)	d lejos de 2
Multicolinealidad	—	`car::vif(m)`	VIF > 5–10
Influencia	`plot(m, which = 4:5)`; `car::influencePlot(m)`	`cooks.distance(m)`, `hatvalues(m)`, `rstudent(m)`	\(D_i>1\); \(h_{ii}>2p/n\); \(\|t_i\|>3\)

\[ VIF_j = \frac{1}{1-R_j^2}, \qquad D_i = \frac{r_i^2}{p}\cdot\frac{h_{ii}}{1-h_{ii}} \quad(\text{Cook, 1977}) \]

Una observación influyente no se elimina automáticamente: se revisa la bitácora, se repite la medición y, si se excluye, se reportan ambos ajustes.

Función auxiliar que se usará en todo el taller:

diagnostico <- function(m) {
  n <- nobs(m); p <- length(coef(m))
  tab <- round(data.frame(h = hatvalues(m), t_est = rstudent(m), D_cook = cooks.distance(m)), 3)
  tab$alerta <- with(tab, ifelse(D_cook > 1 | abs(t_est) > 3 | h > 2 * p / n, "revisar", ""))
  list(shapiro = shapiro.test(rstudent(m)),
       breusch_pagan = bptest(m),
       influencia = tab)
}

5 Ejemplo resuelto 1: rendimiento de una esterificación de Fischer

En la síntesis de acetato de isoamilo (ácido acético + alcohol isoamílico, catalizada por H₂SO₄) se estudió el rendimiento \(Y\) (%) en función de la temperatura \(T\) (°C), la carga de catalizador \(C\) (% m/m) y la razón molar alcohol:ácido \(R\). Datos didácticos simulados a partir de un modelo conocido con error aleatorio.

ester <- data.frame(
  corrida = 1:15,
  Temp = c(60, 60, 60, 60, 70, 70, 70, 70, 80, 80, 80, 80, 70, 70, 70),
  Cat  = c(1, 3, 1, 3, 1, 3, 2, 2, 1, 3, 1, 3, 2, 2, 2),
  Raz  = c(2, 2, 6, 6, 4, 4, 2, 6, 2, 2, 6, 6, 4, 4, 4),
  Y    = c(50.1, 60.2, 56.8, 70.4, 62.1, 69.0, 60.6, 70.0,
           62.0, 70.4, 72.0, 80.0, 65.2, 65.2, 65.6))
ester

corrida	Temp	Cat	Raz	Y
1	60	1	2	50.1
2	60	3	2	60.2
3	60	1	6	56.8
4	60	3	6	70.4
5	70	1	4	62.1
6	70	3	4	69.0
7	70	2	2	60.6
8	70	2	6	70.0
9	80	1	2	62.0
10	80	3	2	70.4
11	80	1	6	72.0
12	80	3	6	80.0
13	70	2	4	65.2
14	70	2	4	65.2
15	70	2	4	65.6

Se usa Temp en lugar de T porque en R T es un alias de TRUE.

5.1 Exploración

pairs(ester[, c("Temp", "Cat", "Raz", "Y")], pch = 19, col = "grey30",
      main = "Matriz de dispersión")

round(cor(ester[, c("Temp", "Cat", "Raz", "Y")]), 3)

#>       Temp   Cat  Raz     Y
#> Temp 1.000 0.000 0.00 0.617
#> Cat  0.000 1.000 0.00 0.553
#> Raz  0.000 0.000 1.00 0.540
#> Y    0.617 0.553 0.54 1.000

Las correlaciones entre predictores son cero: el diseño es ortogonal.

5.2 Ajuste e inferencia

m1 <- lm(Y ~ Temp + Cat + Raz, data = ester)
summary(m1)

#> 
#> Call:
#> lm(formula = Y ~ Temp + Cat + Raz, data = ester)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2.5342 -0.2879 -0.0542  0.4696  1.6658 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   5.6892     3.0751    1.85    0.091 .  
#> Temp          0.5863     0.0411   14.26  1.9e-08 ***
#> Cat           4.7000     0.3678   12.78  6.1e-08 ***
#> Raz           2.2950     0.1839   12.48  7.8e-08 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.16 on 11 degrees of freedom
#> Multiple R-squared:  0.979,  Adjusted R-squared:  0.974 
#> F-statistic:  174 on 3 and 11 DF,  p-value: 1.5e-09

anova(lm(Y ~ 1, data = ester), m1)   # prueba F global

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
14	721.409	NA	NA	NA	NA
11	14.877	3	706.53	174.13	0

confint(m1)

#>                2.5 %   97.5 %
#> (Intercept) -1.07898 12.45732
#> Temp         0.49575  0.67675
#> Cat          3.89057  5.50943
#> Raz          1.89028  2.69972

qf(0.95, 3, 11); qt(0.975, 11)        # valores críticos

#> [1] 3.5874

#> [1] 2.201

Interpretación. \(F_0 = 174.1\) (\(p = 1.5\times10^{-9}\)), \(R^2 = 0.979\), \(R^2_{adj} = 0.974\), \(s = 1.16\) %. Con \(C\) y \(R\) constantes, cada °C aumenta el rendimiento 0.59 puntos (IC 95 %: 0.50 a 0.68); 1 % m/m más de catalizador lo eleva 4.7 puntos y una unidad de razón molar 2.3 puntos, coherente con el desplazamiento del equilibrio por exceso de alcohol (Le Châtelier). El intercepto no tiene sentido físico (0 °C está fuera del dominio).

5.3 Verificación matricial

X <- model.matrix(m1); y <- ester$Y
XtX_inv <- solve(t(X) %*% X)
b  <- XtX_inv %*% t(X) %*% y
e  <- y - X %*% b
s2 <- sum(e^2) / (nrow(X) - ncol(X))
data.frame(b = drop(b), se = sqrt(diag(s2 * XtX_inv)))

	b	se
(Intercept)	5.68917	3.07506
Temp	0.58625	0.04112
Cat	4.70000	0.36776
Raz	2.29500	0.18388

5.4 Diagnóstico

par(mfrow = c(2, 2)); plot(m1); par(mfrow = c(1, 1))

vif(m1)

#> Temp  Cat  Raz 
#>    1    1    1

diagnostico(m1)

#> $shapiro
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  rstudent(m)
#> W = 0.791, p-value = 0.0029
#> 
#> 
#> $breusch_pagan
#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  m
#> BP = 6.33, df = 3, p-value = 0.097
#> 
#> 
#> $influencia
#>        h  t_est D_cook  alerta
#> 1  0.392 -0.057  0.001        
#> 2  0.392  0.695  0.082        
#> 3  0.392 -4.943  1.256 revisar
#> 4  0.392  2.103  0.543        
#> 5  0.167  1.481  0.099        
#> 6  0.167 -0.943  0.045        
#> 7  0.167 -0.105  0.001        
#> 8  0.167  0.093  0.000        
#> 9  0.392  0.127  0.003        
#> 10 0.392 -0.966  0.151        
#> 11 0.392  1.041  0.173        
#> 12 0.392 -0.488  0.041        
#> 13 0.067 -0.091  0.000        
#> 14 0.067 -0.091  0.000        
#> 15 0.067  0.250  0.001

VIF = 1 para los tres factores (diseño ortogonal).
Shapiro-Wilk y Breusch-Pagan no rechazan normalidad ni homocedasticidad al 5 %.
La corrida 3 tiene \(t_3 = -4.94\) y \(D_3 = 1.26 > 1\); su apalancamiento (\(h_{33}=0.392\)) es menor que \(2p/n = 0.533\): influye por su residuo, no por su posición.

Análisis de sensibilidad (no sustituye la revisión de la bitácora):

m1b <- update(m1, subset = corrida != 3)
cbind(con_3 = coef(m1), sin_3 = coef(m1b))

#>               con_3   sin_3
#> (Intercept) 5.68917 9.61192
#> Temp        0.58625 0.53418
#> Cat         4.70000 4.28342
#> Raz         2.29500 2.50329

c(s_con = sigma(m1), s_sin = sigma(m1b))

#>  s_con  s_sin 
#> 1.1630 0.6573

5.5 Estimación y predicción

nuevo <- data.frame(Temp = 75, Cat = 2.5, Raz = 5)
predict(m1, nuevo, interval = "confidence")

#>      fit    lwr    upr
#> 1 72.883 71.898 73.867

predict(m1, nuevo, interval = "prediction")

#>      fit   lwr    upr
#> 1 72.883 70.14 75.625

6 Ejemplo resuelto 2: relación estructura-propiedad en alcanos (modelo de Wiener)

Wiener (1947) relacionó el punto de ebullición normal de los alcanos con dos descriptores del grafo molecular sin hidrógenos:

W (índice de Wiener): suma de las distancias topológicas entre todos los pares de carbonos.
p (número de polaridad): número de pares de carbonos separados por exactamente tres enlaces.

6.1 Cálculo de W y p en R base

Cada molécula se describe por su lista de enlaces C–C; las distancias se obtienen con el algoritmo de Floyd-Warshall.

descriptores <- function(n, enlaces) {
  D <- matrix(Inf, n, n); diag(D) <- 0
  for (e in enlaces) { D[e[1], e[2]] <- 1; D[e[2], e[1]] <- 1 }
  for (k in 1:n) for (i in 1:n) for (j in 1:n)
    D[i, j] <- min(D[i, j], D[i, k] + D[k, j])
  d <- D[upper.tri(D)]
  c(W = sum(d), p = sum(d == 3))
}
cadena <- function(n) lapply(1:(n - 1), function(i) c(i, i + 1))

mols <- list(
  "n-pentano"            = list(5, cadena(5)),
  "2-metilbutano"        = list(5, c(cadena(4), list(c(2, 5)))),
  "2,2-dimetilpropano"   = list(5, list(c(1, 2), c(2, 3), c(2, 4), c(2, 5))),
  "n-hexano"             = list(6, cadena(6)),
  "2-metilpentano"       = list(6, c(cadena(5), list(c(2, 6)))),
  "3-metilpentano"       = list(6, c(cadena(5), list(c(3, 6)))),
  "2,2-dimetilbutano"    = list(6, c(cadena(4), list(c(2, 5), c(2, 6)))),
  "2,3-dimetilbutano"    = list(6, c(cadena(4), list(c(2, 5), c(3, 6)))),
  "n-heptano"            = list(7, cadena(7)),
  "2-metilhexano"        = list(7, c(cadena(6), list(c(2, 7)))),
  "3-metilhexano"        = list(7, c(cadena(6), list(c(3, 7)))),
  "2,2-dimetilpentano"   = list(7, c(cadena(5), list(c(2, 6), c(2, 7)))),
  "2,3-dimetilpentano"   = list(7, c(cadena(5), list(c(2, 6), c(3, 7)))),
  "2,4-dimetilpentano"   = list(7, c(cadena(5), list(c(2, 6), c(4, 7)))),
  "3,3-dimetilpentano"   = list(7, c(cadena(5), list(c(3, 6), c(3, 7)))),
  "3-etilpentano"        = list(7, c(cadena(5), list(c(3, 6), c(6, 7)))),
  "2,2,3-trimetilbutano" = list(7, c(cadena(4), list(c(2, 5), c(2, 6), c(3, 7)))))

desc <- t(sapply(mols, function(m) descriptores(m[[1]], m[[2]])))
# Puntos de ebullición normales (°C): CRC Handbook of Chemistry and Physics (Haynes, 2011)
alcanos <- data.frame(
  compuesto = names(mols),
  nC  = sapply(mols, `[[`, 1),
  W   = desc[, "W"], p = desc[, "p"],
  Teb = c(36.06, 27.85, 9.50, 68.73, 60.26, 63.27, 49.73, 57.98,
          98.38, 90.05, 91.95, 79.20, 89.78, 80.49, 86.06, 93.50, 80.86),
  row.names = NULL)
alcanos

compuesto	nC	W	p	Teb
n-pentano	5	20	2	36.06
2-metilbutano	5	18	2	27.85
2,2-dimetilpropano	5	16	0	9.50
n-hexano	6	35	3	68.73
2-metilpentano	6	32	3	60.26
3-metilpentano	6	31	4	63.27
2,2-dimetilbutano	6	28	3	49.73
2,3-dimetilbutano	6	29	4	57.98
n-heptano	7	56	4	98.38
2-metilhexano	7	52	4	90.05
3-metilhexano	7	50	5	91.95
2,2-dimetilpentano	7	46	4	79.20
2,3-dimetilpentano	7	46	6	89.78
2,4-dimetilpentano	7	48	4	80.49
3,3-dimetilpentano	7	44	6	86.06
3-etilpentano	7	48	6	93.50
2,2,3-trimetilbutano	7	42	6	80.86

Verificación: n-pentano, \(W = 4(1)+3(2)+2(3)+1(4) = 20\) y \(p = 2\).

6.2 Correlaciones y comparación de modelos

round(cor(alcanos[, c("nC", "W", "p", "Teb")]), 3)

#>        nC     W     p   Teb
#> nC  1.000 0.961 0.851 0.953
#> W   0.961 1.000 0.753 0.966
#> p   0.851 0.753 1.000 0.865
#> Teb 0.953 0.966 0.865 1.000

formulas <- c("Teb ~ nC", "Teb ~ nC + p", "Teb ~ nC + W", "Teb ~ W + p", "Teb ~ nC + W + p")
comparacion <- do.call(rbind, lapply(formulas, function(f) {
  m <- lm(as.formula(f), data = alcanos)
  data.frame(modelo = f, R2 = summary(m)$r.squared, R2adj = summary(m)$adj.r.squared,
             s = sigma(m), AIC = AIC(m),
             VIF_max = if (length(coef(m)) > 2) max(vif(m)) else NA)
}))
comparacion$R2 <- round(comparacion$R2, 3); comparacion$R2adj <- round(comparacion$R2adj, 3)
comparacion$s <- round(comparacion$s, 2); comparacion$AIC <- round(comparacion$AIC, 1)
comparacion$VIF_max <- round(comparacion$VIF_max, 1)
comparacion

modelo	R2	R2adj	s	AIC	VIF_max
Teb ~ nC	0.908	0.902	8.01	122.9	NA
Teb ~ nC + p	0.918	0.907	7.81	122.8	3.6
Teb ~ nC + W	0.942	0.933	6.60	117.1	13.0
Teb ~ W + p	0.977	0.974	4.14	101.2	2.3
Teb ~ nC + W + p	0.983	0.980	3.64	97.6	25.4

coef(lm(Teb ~ nC + W + p, data = alcanos))

#> (Intercept)          nC           W           p 
#>     47.1168    -13.1498      2.0882      6.7324

El modelo completo tiene el mejor \(R^2_{adj}\), pero \(b(nC) = -13.1\) °C: a \(W\) y \(p\) constantes, agregar un carbono “bajaría” el punto de ebullición, lo que carece de sentido físico. Es multicolinealidad (\(r(nC,W)=0.961\), VIF = 25.4). Se elige Teb ~ W + p por parsimonia, estabilidad e interpretabilidad.

6.3 Modelo seleccionado y diagnóstico

m2 <- lm(Teb ~ W + p, data = alcanos)
summary(m2); confint(m2); vif(m2)

#> 
#> Call:
#> lm(formula = Teb ~ W + p, data = alcanos)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -7.833 -2.108 -0.169  2.262  8.599 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   -6.363      3.266   -1.95  0.07173 .  
#> W              1.481      0.125   11.83  1.1e-08 ***
#> p              4.886      0.951    5.14  0.00015 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 4.14 on 14 degrees of freedom
#> Multiple R-squared:  0.977,  Adjusted R-squared:  0.974 
#> F-statistic:  298 on 2 and 14 DF,  p-value: 3.35e-12

#>                2.5 % 97.5 %
#> (Intercept) -13.3676 0.6422
#> W             1.2125 1.7495
#> p             2.8468 6.9262

#>      W      p 
#> 2.3112 2.3112

d2 <- diagnostico(m2)
d2$shapiro

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  rstudent(m)
#> W = 0.945, p-value = 0.38

cbind(compuesto = alcanos$compuesto, d2$influencia)[d2$influencia$alerta != "", ]

	compuesto	h	t_est	D_cook	alerta
3	2,2-dimetilpropano	0.403	-3.125	1.354	revisar

El 2,2-dimetilpropano (neopentano) es influyente: \(t=-3.13\), \(h=0.403 > 2p/n = 0.353\), \(D = 1.35\). El dato es correcto; es el isómero más esférico y el extremo del dominio de \(W\). Se conserva y se informa la limitación.

6.4 Extrapolación a C₈

c8 <- list("n-octano"              = list(8, cadena(8)),
           "2,2,4-trimetilpentano" = list(8, c(cadena(5), list(c(2, 6), c(2, 7), c(4, 8)))))
nuevos <- data.frame(t(sapply(c8, function(m) descriptores(m[[1]], m[[2]]))))
pred <- predict(m2, nuevos, interval = "prediction")
data.frame(compuesto = names(c8), nuevos, Teb_real = c(125.6, 99.2), round(pred, 1))

	compuesto	W	p	Teb_real	fit	lwr	upr
n-octano	n-octano	84	5	125.6	142.5	128.3	156.6
2,2,4-trimetilpentano	2,2,4-trimetilpentano	66	5	99.2	115.8	104.8	126.8

plot(alcanos$Teb, fitted(m2), pch = 19, col = "grey55",
     xlim = c(0, 160), ylim = c(0, 160),
     xlab = "Teb observada (°C)", ylab = "Teb predicha (°C)",
     main = "Modelo W + p: ajuste C5–C7 y extrapolación a C8")
abline(0, 1, lty = 2, col = "grey40")
arrows(c(125.6, 99.2), pred[, "lwr"], c(125.6, 99.2), pred[, "upr"],
       angle = 90, code = 3, length = 0.05, col = "steelblue", lwd = 2)
points(c(125.6, 99.2), pred[, "fit"], pch = 19, col = "steelblue", cex = 1.4)
text(c(125.6, 99.2), pred[, "fit"], names(c8), pos = 2, cex = 0.8)
legend("topleft", c("C5–C7 (ajuste)", "C8 extrapolado ± IP 95 %"),
       pch = 19, col = c("grey55", "steelblue"), bty = "n")

Ambos valores reales quedan fuera del intervalo de predicción: la relación Teb–W no es lineal en un intervalo amplio de tamaños y el modelo solo vale dentro de su dominio.

7 Problemas de práctica

Datos de los problemas 1, 3 y 4: didácticos simulados. Use \(\alpha = 0.05\).

7.1 Problema 1. Cálculo matricial

La señal \(Y\) (mV) de un sensor se registró para seis disoluciones con analito \(X_1\) (mmol/L) e interferente \(X_2\) (mmol/L).

sensor <- data.frame(X1 = 1:6, X2 = c(2, 1, 4, 3, 6, 5), Y = c(5, 6, 10, 10, 15, 15))
sensor

X1	X2	Y
1	2	5
2	1	6
3	4	10
4	3	10
5	6	15
6	5	15

Construya \(\mathbf X\) con cbind() y calcule \(\mathbf X^{\mathsf T}\mathbf X\), \(\mathbf X^{\mathsf T}\mathbf y\) y \(\det(\mathbf X^{\mathsf T}\mathbf X)\).
Obtenga \(\mathbf b\) con solve() y exprese los coeficientes como fracciones (MASS::fractions()).
Calcule \(\hat{\mathbf y}\), residuos, \(SS_E\), \(s^2\) y \(R^2\) sin usar lm().
Calcule \(se(b_1)\) y pruebe \(H_0:\beta_1=0\). Verifique todo con lm().
Construya la tabla ANOVA y realice la prueba F global.

7.2 Problema 2. Completar una tabla ANOVA

Se modeló la DQO (mg/L) de 25 muestras de efluente de un beneficio de café con cuatro predictores. El software reportó solo \(SS_T = 480.0\) y \(SS_E = 60.0\).

Complete la tabla ANOVA y realice la prueba F (obtenga el valor crítico con qf()).
Calcule \(R^2\), \(R^2_{adj}\) y \(s\).
Al agregar un quinto predictor, \(R^2\) sube a 0.880. Calcule el nuevo \(R^2_{adj}\) y la prueba F parcial. ¿Conviene incluirlo?

7.3 Problema 3. Calibración multicomponente: dicromato y permanganato

Absorbancias de mezclas de K₂Cr₂O₇ (\(c_D\)) y KMnO₄ (\(c_M\)), en mmol/L, a 440 y 545 nm (celda de 1.00 cm). Por aditividad de la ley de Beer: \(A_\lambda = \varepsilon_{D,\lambda} b\, c_D + \varepsilon_{M,\lambda} b\, c_M\).

mezclas <- data.frame(
  cD   = c(0, 0.5, 1, 1.5, 2, 0, 0.5, 1, 1.5, 2, 0.5, 1.5),
  cM   = c(0.1, 0.1, 0.1, 0.1, 0.1, 0.3, 0.3, 0.3, 0.3, 0.3, 0.2, 0.2),
  A440 = c(0.014, 0.199, 0.381, 0.563, 0.750, 0.029, 0.217, 0.407, 0.584, 0.768, 0.209, 0.578),
  A545 = c(0.238, 0.240, 0.249, 0.257, 0.255, 0.706, 0.706, 0.714, 0.717, 0.729, 0.473, 0.491))

Ajuste A440 ~ cD + cM y A545 ~ cD + cM. Reporte coeficientes, errores estándar, \(R^2\) y \(s\).
¿Es el intercepto distinto de cero? ¿Qué significa químicamente?
Exprese las pendientes como absortividades molares \(\varepsilon\) (L mol⁻¹ cm⁻¹).
Una muestra da \(A_{440}=0.512\) y \(A_{545}=0.583\). Determine \(c_D\) y \(c_M\) resolviendo el sistema con solve().

7.4 Problema 4. Ley de velocidad y energía de activación

Velocidades iniciales \(r\) (mol L⁻¹ s⁻¹) de A + B → productos. Modelo: \(r = k_0\, e^{-E_a/RT}[A]^a[B]^b\) con error multiplicativo.

cinetica <- data.frame(
  TK = rep(c(298.15, 308.15, 318.15, 328.15), each = 4),
  A  = rep(c(0.05, 0.10, 0.05, 0.10), 4),
  B  = rep(c(0.05, 0.05, 0.20, 0.20), 4),
  r  = c(0.000926, 0.001950, 0.001940, 0.003620, 0.001950, 0.003870, 0.004040, 0.007880,
         0.004150, 0.007480, 0.008570, 0.016000, 0.008070, 0.015600, 0.015500, 0.032000))

Linealice el modelo y ajústelo con lm() usando I(1/TK), log(A) y log(B).
Estime \(E_a\) (kJ/mol) con su IC 95 % (\(R = 8.314\) J mol⁻¹ K⁻¹).
Pruebe \(H_0: a = 1\) y \(H_0: b = 0.5\) (sugerencia: car::linearHypothesis() o el estadístico t manual).
Estime \(k_0\) con sus unidades. ¿Por qué su incertidumbre es tan grande?

7.5 Problema 5. Multicolinealidad y extrapolación oculta

Si \(r(x_1,x_2) = 0.95\), calcule el VIF de cada predictor.
Las regresiones auxiliares de un modelo dan \(R^2_1=0.20\), \(R^2_2=0.75\) y \(R^2_3=0.92\). Calcule los VIF.
Explique por qué \(b(nC)\) resultó negativo en el modelo completo del Ejemplo 2.
Con m2, ¿es confiable predecir Teb para \(W=56\) y \(p=0\)? Calcule \(h_{00}=\mathbf x_0^{\mathsf T}(\mathbf X^{\mathsf T}\mathbf X)^{-1}\mathbf x_0\) y compárelo con max(hatvalues(m2)).

8 Soluciones

Valores calculados con R. Diferencias en el último decimal por redondeo manual son aceptables.

8.1 Solución 1

X <- cbind(1, sensor$X1, sensor$X2); y <- sensor$Y
XtX <- t(X) %*% X; Xty <- t(X) %*% y
XtX; Xty; det(XtX)

#>      [,1] [,2] [,3]
#> [1,]    6   21   21
#> [2,]   21   91   88
#> [3,]   21   88   91

#>      [,1]
#> [1,]   61
#> [2,]  252
#> [3,]  251

#> [1] 576

b <- solve(XtX, Xty); MASS::fractions(drop(b))      # 89/48, 65/48, 49/48

#> [1] 89/48 65/48 49/48

yhat <- drop(X %*% b); e <- y - yhat
SSE <- sum(e^2); SST <- sum((y - mean(y))^2); s2 <- SSE / (6 - 3)
c(SSE = SSE, s2 = s2, R2 = 1 - SSE / SST)

#>     SSE      s2      R2 
#> 0.41667 0.13889 0.99541

se_b1 <- sqrt(s2 * solve(XtX)[2, 2]); t0 <- b[2] / se_b1
c(se_b1 = se_b1, t0 = t0, t_crit = qt(0.975, 3), p = 2 * pt(-abs(t0), 3))

#>     se_b1        t0    t_crit         p 
#> 0.1591172 8.5104977 3.1824463 0.0034075

mP1 <- lm(Y ~ X1 + X2, data = sensor); summary(mP1)

#> 
#> Call:
#> lm(formula = Y ~ X1 + X2, data = sensor)
#> 
#> Residuals:
#>         1         2         3         4         5         6 
#> -2.50e-01  4.17e-01 -1.30e-15 -3.33e-01  2.50e-01 -8.33e-02 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)   
#> (Intercept)    1.854      0.360    5.15   0.0142 * 
#> X1             1.354      0.159    8.51   0.0034 **
#> X2             1.021      0.159    6.42   0.0077 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.373 on 3 degrees of freedom
#> Multiple R-squared:  0.995,  Adjusted R-squared:  0.992 
#> F-statistic:  326 on 2 and 3 DF,  p-value: 0.000311

anova(lm(Y ~ 1, data = sensor), mP1); qf(0.95, 2, 3)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
5	90.83333	NA	NA	NA	NA
3	0.41667	2	90.417	325.5	0.00031

#> [1] 9.5521

\(b_0 = 89/48 = 1.854\), \(b_1 = 65/48 = 1.354\), \(b_2 = 49/48 = 1.021\); \(SS_E = 0.4167\), \(s^2 = 0.1389\), \(R^2 = 0.9954\). \(t_0(b_1) = 8.51 > 3.182\) (\(p = 0.0034\)). \(F_0 = 325.5 > F_{0.05;2,3} = 9.552\). El interferente también es significativo (\(p = 0.0077\)): el sensor no es selectivo.

8.2 Solución 2

n <- 25; k <- 4; SST <- 480; SSE <- 60; SSR <- SST - SSE
tabla <- data.frame(gl = c(k, n - k - 1, n - 1), SC = c(SSR, SSE, SST),
                    CM = c(SSR / k, SSE / (n - k - 1), NA),
                    row.names = c("Regresión", "Residual", "Total"))
tabla$F0 <- c(tabla$CM[1] / tabla$CM[2], NA, NA); tabla

	gl	SC	CM	F0
Regresión	4	420	105	35
Residual	20	60	3	NA
Total	24	480	NA	NA

c(F_crit = qf(0.95, 4, 20), p = pf(35, 4, 20, lower.tail = FALSE))

#>     F_crit          p 
#> 2.8661e+00 9.0804e-09

c(R2 = SSR / SST, R2adj = 1 - (SSE / 20) / (SST / 24), s = sqrt(SSE / 20))

#>     R2  R2adj      s 
#> 0.8750 0.8500 1.7321

SSE5 <- SST * (1 - 0.880)
c(R2adj_5 = 1 - (SSE5 / 19) / (SST / 24),
  F_parcial = (SSE - SSE5) / (SSE5 / 19), F_crit = qf(0.95, 1, 19))

#>   R2adj_5 F_parcial    F_crit 
#>   0.84842   0.79167   4.38075

\(F_0 = 35.0\) (\(p = 9.1\times10^{-9}\)); \(R^2 = 0.875\), \(R^2_{adj} = 0.850\), \(s = 1.73\) mg/L. Con el quinto predictor, \(R^2_{adj} = 0.848 < 0.850\) y \(F_{parcial} = 0.79 < 4.381\): no conviene incluirlo.

8.3 Solución 3

m440 <- lm(A440 ~ cD + cM, data = mezclas); m545 <- lm(A545 ~ cD + cM, data = mezclas)
summary(m440)$coefficients; summary(m545)$coefficients

#>              Estimate Std. Error  t value   Pr(>|t|)
#> (Intercept) 0.0038015  0.0023121   1.6442 1.3456e-01
#> cD          0.3681818  0.0012088 304.5728 2.2567e-19
#> cM          0.0980000  0.0089650  10.9313 1.6972e-06

#>              Estimate Std. Error  t value   Pr(>|t|)
#> (Intercept) 0.0031955  0.0026873   1.1891 2.6483e-01
#> cD          0.0114545  0.0014050   8.1527 1.9025e-05
#> cM          2.3330000  0.0104197 223.9022 3.5977e-18

c(R2_440 = summary(m440)$r.squared, s_440 = sigma(m440),
  R2_545 = summary(m545)$r.squared, s_545 = sigma(m545))

#>   R2_440    s_440   R2_545    s_545 
#> 0.999903 0.002835 0.999821 0.003295

confint(m440)["(Intercept)", ]; confint(m545)["(Intercept)", ]

#>      2.5 %     97.5 % 
#> -0.0014289  0.0090319

#>      2.5 %     97.5 % 
#> -0.0028836  0.0092746

epsilon <- 1000 * rbind(A440 = coef(m440)[-1], A545 = coef(m545)[-1]); round(epsilon, 1)

#>         cD   cM
#> A440 368.2   98
#> A545  11.5 2333

K <- rbind(coef(m440)[-1], coef(m545)[-1])
solve(K, c(0.512, 0.583) - c(coef(m440)[1], coef(m545)[1]))

#>      cD      cM 
#> 1.31586 0.24206

Los interceptos no difieren de cero (\(p = 0.135\) y \(0.265\)): sin sesgo del blanco ni desviación de la ley de Beer. \(\varepsilon_{D,440} = 368\), \(\varepsilon_{M,440} = 98\), \(\varepsilon_{D,545} = 11.5\), \(\varepsilon_{M,545} = 2333\) L mol⁻¹ cm⁻¹. Muestra desconocida: \(c_D = 1.32\) mmol/L y \(c_M = 0.242\) mmol/L.

8.4 Solución 4

m4 <- lm(log(r) ~ I(1/TK) + log(A) + log(B), data = cinetica)
summary(m4)

#> 
#> Call:
#> lm(formula = log(r) ~ I(1/TK) + log(A) + log(B), data = cinetica)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.05726 -0.01998 -0.00382  0.01859  0.06189 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  2.07e+01   2.44e-01    84.9  < 2e-16 ***
#> I(1/TK)     -6.94e+03   7.31e+01   -95.0  < 2e-16 ***
#> log(A)       9.59e-01   2.41e-02    39.8  4.1e-14 ***
#> log(B)       5.10e-01   1.20e-02    42.3  2.0e-14 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.0334 on 12 degrees of freedom
#> Multiple R-squared:  0.999,  Adjusted R-squared:  0.999 
#> F-statistic: 4.13e+03 on 3 and 12 DF,  p-value: <2e-16

R <- 8.314
Ea <- -coef(m4)[2] * R / 1000; Ea_IC <- -rev(confint(m4)[2, ]) * R / 1000
c(Ea_kJ = unname(Ea), IC_inf = unname(Ea_IC[1]), IC_sup = unname(Ea_IC[2]))

#>  Ea_kJ IC_inf IC_sup 
#> 57.711 56.387 59.034

linearHypothesis(m4, "log(A) = 1"); linearHypothesis(m4, "log(B) = 0.5")

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
13	0.01657	NA	NA	NA	NA
12	0.01339	1	0.00318	2.8535	0.11696

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
13	0.01414	NA	NA	NA	NA
12	0.01339	1	0.00075	0.6678	0.42974

cf <- summary(m4)$coefficients
c(t_a = (cf[3, 1] - 1) / cf[3, 2], t_b = (cf[4, 1] - 0.5) / cf[4, 2], t_crit = qt(0.975, 12))

#>      t_a      t_b   t_crit 
#> -1.68924  0.81719  2.17881

c(k0 = exp(coef(m4)[1]), exp(confint(m4)[1, ]))

#> k0.(Intercept)          2.5 %         97.5 % 
#>      992659407      583444866     1688887426

par(mfrow = c(2, 2)); plot(m4); par(mfrow = c(1, 1))

\(E_a = 57.7\) kJ/mol (IC 95 %: 56.4 a 59.0). \(a = 0.959\) (\(t = -1.69\)) y \(b = 0.510\) (\(t = 0.82\)): compatibles con \(a = 1\) y \(b = 1/2\) (orden global 3/2). \(k_0 \approx 9.9\times10^{8}\) L\(^{0.5}\) mol\(^{-0.5}\) s\(^{-1}\); su incertidumbre es grande porque \(1/T = 0\) está muy lejos de los datos. Para el docente: los datos se generaron con \(E_a = 58.0\) kJ/mol, \(a = 1\) y \(b = 0.5\).

8.5 Solución 5

1 / (1 - 0.95^2)

#> [1] 10.256

1 / (1 - c(0.20, 0.75, 0.92))

#> [1]  1.25  4.00 12.50

X2 <- model.matrix(m2); x0 <- c(1, 56, 0)
c(h00 = drop(t(x0) %*% solve(t(X2) %*% X2) %*% x0), h_max = max(hatvalues(m2)))

#>     h00   h_max 
#> 1.90488 0.40346

predict(m2, data.frame(W = 56, p = 0), interval = "prediction")

#>      fit    lwr  upr
#> 1 76.572 61.445 91.7

VIF = 10.3 para ambos. b) VIF = 1.25, 4.0 y 12.5: \(x_3\) es problemático. c) \(r(nC,W) = 0.961\): el coeficiente de \(nC\) compensa a \(W\) y se vuelve inestable (VIF = 25.4); “más carbonos con \(W\) y \(p\) fijos” es una situación que casi no existe en los datos. d) No: cada valor está en su rango, pero la combinación no existe (W = 56 solo para n-heptano, con p = 4; p = 0 solo para neopentano, con W = 16). \(h_{00} = 1.90\), casi cinco veces el máximo de los datos (0.40): extrapolación oculta.

9 Referencias

Breusch, T. S., & Pagan, A. R. (1979). A simple test for heteroscedasticity and random coefficient variation. Econometrica, 47(5), 1287–1294. https://doi.org/10.2307/1911963
Brereton, R. G. (2003). Chemometrics: Data analysis for the laboratory and chemical plant. Wiley.
Cook, R. D. (1977). Detection of influential observation in linear regression. Technometrics, 19(1), 15–18. https://doi.org/10.1080/00401706.1977.10489493
Draper, N. R., & Smith, H. (1998). Applied regression analysis (3.ª ed.). Wiley. https://doi.org/10.1002/9781118625590
Durbin, J., & Watson, G. S. (1950). Testing for serial correlation in least squares regression I. Biometrika, 37(3/4), 409–428. https://doi.org/10.2307/2332391
Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3.ª ed.). SAGE. (Paquete car.)
Haynes, W. M. (Ed.). (2011). CRC handbook of chemistry and physics (92.ª ed.). CRC Press.
Kutner, M. H., Nachtsheim, C. J., Neter, J., & Li, W. (2005). Applied linear statistical models (5.ª ed.). McGraw-Hill/Irwin.
Miller, J. N., & Miller, J. C. (2018). Statistics and chemometrics for analytical chemistry (7.ª ed.). Pearson.
Montgomery, D. C., Peck, E. A., & Vining, G. G. (2012). Introduction to linear regression analysis (5.ª ed.). Wiley.
R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/
Shapiro, S. S., & Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3–4), 591–611. https://doi.org/10.1093/biomet/52.3-4.591
Wiener, H. (1947). Structural determination of paraffin boiling points. Journal of the American Chemical Society, 69(1), 17–20. https://doi.org/10.1021/ja01193a005
Xie, Y., Allaire, J. J., & Grolemund, G. (2018). R Markdown: The definitive guide. Chapman & Hall/CRC. https://bookdown.org/yihui/rmarkdown/
Zeileis, A., & Hothorn, T. (2002). Diagnostic checking in regression relationships. R News, 2(3), 7–10. (Paquete lmtest.)

#> Compilado con R version 4.6.1 (2026-06-24 ucrt) | car 3.1.5 | lmtest 0.9.40

Taller de Regresión Lineal Múltiple con R

Teoría, ejemplos resueltos y problemas de práctica con datos químicos

Dr. Pedro González Beermann — Escuela de Química, UNACHI

27/10/2026

1 Presentación

2 Fundamento teórico

2.1 El modelo

2.2 Notación matricial y mínimos cuadrados

2.3 Supuestos

Demostración en R: MCO “a mano” frente a `lm()`

3 Inferencia

3.1 Tabla ANOVA y prueba F global

3.2 Pruebas t, R² y R² ajustado

3.3 Intervalos

4 Diagnóstico

5 Ejemplo resuelto 1: rendimiento de una esterificación de Fischer

5.1 Exploración

5.2 Ajuste e inferencia

5.3 Verificación matricial

5.4 Diagnóstico

5.5 Estimación y predicción

6 Ejemplo resuelto 2: relación estructura-propiedad en alcanos (modelo de Wiener)

6.1 Cálculo de W y p en R base

6.2 Correlaciones y comparación de modelos

6.3 Modelo seleccionado y diagnóstico

6.4 Extrapolación a C₈

7 Problemas de práctica

7.1 Problema 1. Cálculo matricial

7.2 Problema 2. Completar una tabla ANOVA

7.3 Problema 3. Calibración multicomponente: dicromato y permanganato

7.4 Problema 4. Ley de velocidad y energía de activación

7.5 Problema 5. Multicolinealidad y extrapolación oculta

8 Soluciones

8.1 Solución 1

8.2 Solución 2

8.3 Solución 3

8.4 Solución 4

8.5 Solución 5

9 Referencias

Taller de Regresión Lineal Múltiple con R

Teoría, ejemplos resueltos y problemas de práctica con datos químicos

Dr. Pedro González Beermann — Escuela de Química, UNACHI

27/10/2026

1 Presentación

2 Fundamento teórico

2.1 El modelo

2.2 Notación matricial y mínimos cuadrados

2.3 Supuestos

Demostración en R: MCO “a mano” frente a lm()

3 Inferencia

3.1 Tabla ANOVA y prueba F global

3.2 Pruebas t, R² y R² ajustado

3.3 Intervalos

4 Diagnóstico

5 Ejemplo resuelto 1: rendimiento de una esterificación de Fischer

5.1 Exploración

5.2 Ajuste e inferencia

5.3 Verificación matricial

5.4 Diagnóstico

5.5 Estimación y predicción

6 Ejemplo resuelto 2: relación estructura-propiedad en alcanos (modelo de Wiener)

6.1 Cálculo de W y p en R base

6.2 Correlaciones y comparación de modelos

6.3 Modelo seleccionado y diagnóstico

6.4 Extrapolación a C₈

7 Problemas de práctica

7.1 Problema 1. Cálculo matricial

7.2 Problema 2. Completar una tabla ANOVA

7.3 Problema 3. Calibración multicomponente: dicromato y permanganato

7.4 Problema 4. Ley de velocidad y energía de activación

7.5 Problema 5. Multicolinealidad y extrapolación oculta

8 Soluciones

8.1 Solución 1

8.2 Solución 2

8.3 Solución 3

8.4 Solución 4

8.5 Solución 5

9 Referencias

Demostración en R: MCO “a mano” frente a `lm()`