El peso de una persona no solo depende de su altura, sino también de su edad. Suponemos que ambas relaciones son lineales y aditivas, es decir:

En la realidad, personas con la misma altura y edad no tienen exactamente el mismo peso, por lo que añadimos un término de error. El modelo se escribe como: \[\begin{equation} \text{Peso}_i = \beta_0 + \beta_1 \cdot \text{Altura}_i + \beta_2 \cdot \text{Edad}_i + \varepsilon_i \label{eq:model} \end{equation}\]

El objetivo es encontrar los coeficientes \(\beta_0\), \(\beta_1\) y \(\beta_2\) que se ajustan a los datos. Usamos el criterio de : minimizar la suma de los errores al cuadrado. \[\begin{equation} E(\beta_0,\beta_1,\beta_2) = \sum_{i=1}^n (\text{peso}_i - \beta_0 - \beta_1 \text{altura}_i - \beta_2 \text{edad}_i)^2 \label{eq:error} \end{equation}\]

Para minimizar \(E\) calculamos las derivadas parciales respecto a cada parámetro y las igualamos a cero (el gradiente se anula en el mínimo): \[\begin{align} \frac{\partial E}{\partial \beta_0} &= -2 \sum_{i=1}^n (\text{peso}_i - \beta_0 - \beta_1 \text{altura}_i - \beta_2 \text{edad}_i) = 0 \label{eq:d0}\\ \frac{\partial E}{\partial \beta_1} &= -2 \sum_{i=1}^n (\text{peso}_i - \beta_0 - \beta_1 \text{altura}_i - \beta_2 \text{edad}_i) \cdot \text{altura}_i = 0 \label{eq:d1}\\ \frac{\partial E}{\partial \beta_2} &= -2 \sum_{i=1}^n (\text{peso}_i - \beta_0 - \beta_1 \text{altura}_i - \beta_2 \text{edad}_i) \cdot \text{edad}_i = 0 \label{eq:d2} \end{align}\]

Simplificando se obtienen las para el caso de dos predictores: \[\begin{align} \sum \text{peso}_i &= n\beta_0 + \beta_1\sum \text{altura}_i + \beta_2\sum \text{edad}_i \label{eq:n1}\\ \sum \text{peso}_i\text{altura}_i &= \beta_0\sum \text{altura}_i + \beta_1\sum \text{altura}_i^2 + \beta_2\sum \text{altura}_i\text{edad}_i \label{eq:n2}\\ \sum \text{peso}_i\text{edad}_i &= \beta_0\sum \text{edad}_i + \beta_1\sum \text{altura}_i\text{edad}_i + \beta_2\sum \text{edad}_i^2 \label{eq:n3} \end{align}\]

Este sistema lineal de tres ecuaciones con tres incógnitas se puede escribir elegantemente en notación matricial. Definiendo: \[ \mathbf{X} = \begin{bmatrix} 1 & \text{altura}_1 & \text{edad}_1 \\ 1 & \text{altura}_2 & \text{edad}_2 \\ \vdots & \vdots & \vdots \\ 1 & \text{altura}_n & \text{edad}_n \end{bmatrix}, \quad \boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{bmatrix}, \quad \mathbf{y} = \begin{bmatrix} \text{peso}_1 \\ \text{peso}_2 \\ \vdots \\ \text{peso}_n \end{bmatrix} \] el modelo es \(\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\) y el error cuadrático \(E = (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\). Derivando e igualando a cero se obtienen las ecuaciones normales: \[ \mathbf{X}^\top \mathbf{X}\,\boldsymbol{\beta} = \mathbf{X}^\top \mathbf{y} \quad\Rightarrow\quad \hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y} \] que es la generalización directa de la solución univariada.

Supongamos que la relación teórica (sin error) es: \[ \text{Peso}_{\text{teórico}} = -110 + 1.0 \cdot \text{Altura} + 0.2 \cdot \text{Edad} \]

Generamos una población de \(n = 10\,000\) individuos.

set.seed(123)
n <- 10000

# Predictores
altura <- rnorm(n, mean = 170, sd = 15)
edad   <- runif(n, min = 18, max = 80)

# Error
error <- rnorm(n, mean = 0, sd = 5)

# Coeficientes verdaderos
beta0_true <- -110
beta1_true <- 1.0
beta2_true <- 0.2

# Peso observado
peso <- beta0_true + beta1_true * altura + beta2_true * edad + error

# Resumen de las variables
summary(altura)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   112.3   160.0   169.8   170.0   180.1   227.7
summary(edad)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.02   33.61   48.86   49.04   64.85   80.00
summary(peso)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.25   58.97   69.62   69.68   80.62  131.87
# Correlaciones
cor(cbind(altura, edad, peso))
##             altura        edad      peso
## altura  1.00000000 -0.01012044 0.9254573
## edad   -0.01012044  1.00000000 0.2054538
## peso    0.92545729  0.20545384 1.0000000

Al tener dos predictores, la nube de puntos es tridimensional. Para obtener una idea gráfica podemos colorear los puntos según la edad en un diagrama de dispersión de altura vs peso, o bien examinar cada predictor por separado.

paleta <- colorRampPalette(c("blue", "red"))(100)
color_edad <- paleta[cut(edad, breaks = 100, labels = FALSE)]

plot(altura, peso,
     col = color_edad, pch = 19, cex = 0.5,
     xlab = "Altura (cm)", ylab = "Peso (kg)",
     main = "Altura vs Peso (color = Edad)")
legend("topleft", legend = c("Joven (18)", "Mayor (80)"),
       fill = c("blue", "red"), cex = 0.8)

par(mfrow = c(1,2))
plot(altura, peso, pch = 19, col = rgb(0,0,1,0.3),
     xlab = "Altura (cm)", ylab = "Peso (kg)", main = "Peso vs Altura")
abline(lm(peso ~ altura), col = "red", lwd = 2)
plot(edad, peso, pch = 19, col = rgb(0,0,1,0.3),
     xlab = "Edad (años)", ylab = "Peso (kg)", main = "Peso vs Edad")
abline(lm(peso ~ edad), col = "red", lwd = 2)

par(mfrow = c(1,1))

(Ejecute el código anterior en R para generar las figuras; las rectas de regresión simple no capturan el efecto conjunto.)

Estimamos los coeficientes \(\beta_0, \beta_1, \beta_2\) por mínimos cuadrados usando la función de R y también de forma matricial para ilustrar la fórmula \(\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}\).

# Ajuste con lm()
modelo <- lm(peso ~ altura + edad)
summary(modelo)
## 
## Call:
## lm(formula = peso ~ altura + edad)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.3883  -3.4273  -0.0237   3.3815  21.4896 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.093e+02  5.902e-01 -185.19   <2e-16 ***
## altura       9.975e-01  3.356e-03  297.21   <2e-16 ***
## edad         1.926e-01  2.798e-03   68.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.027 on 9997 degrees of freedom
## Multiple R-squared:  0.9026, Adjusted R-squared:  0.9026 
## F-statistic: 4.633e+04 on 2 and 9997 DF,  p-value: < 2.2e-16
# Cálculo matricial directo
X <- cbind(1, altura, edad)
beta_hat_mat <- solve(t(X) %*% X, t(X) %*% peso)
rownames(beta_hat_mat) <- c("beta0", "beta1", "beta2")

cat("Coeficientes estimados (matricial):\n")
## Coeficientes estimados (matricial):
print(beta_hat_mat)
##               [,1]
## beta0 -109.2983637
## beta1    0.9974829
## beta2    0.1925806
cat("\nValores teóricos: beta0 = -110, beta1 = 1, beta2 = 0.2\n")
## 
## Valores teóricos: beta0 = -110, beta1 = 1, beta2 = 0.2

Los tres coeficientes estimados se acercan mucho a los teóricos. El pequeño error es puramente muestral y disminuye al aumentar \(n\).

Podemos comprobar la calidad del ajuste observando los valores ajustados frente a los reales y los residuos.

peso_ajustado <- fitted(modelo)
plot(peso_ajustado, peso,
     pch = 19, col = rgb(0,0.5,0,0.3),
     xlab = "Peso predicho (kg)", ylab = "Peso observado (kg)",
     main = "Ajuste del modelo múltiple")
abline(a = 0, b = 1, col = "red", lwd = 2)

plot(peso_ajustado, resid(modelo),
     pch = 19, col = rgb(0.8,0,0,0.3),
     xlab = "Peso predicho", ylab = "Residuo",
     main = "Gráfico de residuos")
abline(h = 0, col = "blue", lwd = 2)

Los residuos se distribuyen alrededor de cero sin un patrón claro, lo que indica que el modelo lineal con dos predictores captura razonablemente la relación.

Hemos extendido el modelo de regresión lineal simple (un predictor) al caso de dos predictores (altura y edad) usando el mismo principio de mínimos cuadrados. La derivación mediante el gradiente produce un sistema de ecuaciones normales cuya solución matricial es \(\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}\). Con una simulación en R recuperamos los parámetros teóricos con gran precisión, y los diagnósticos visuales confirman el buen ajuste. Este esquema se generaliza directamente a cualquier número de predictores, dando lugar a los modelos de regresión lineal múltiple.