1 Introducción

1.1 ¿Qué es la Econometría?

La Econometría es la rama de la economía que combina teoría económica, matemáticas (particularmente álgebra lineal y cálculo) y estadística para medir relaciones económicas, contrastar teorías y pronosticar variables de interés a partir de datos observados. Mientras la teoría económica postula relaciones cualitativas (por ejemplo, “a mayor educación, mayor salario”), la econometría busca cuantificar esa relación: ¿en cuánto aumenta el salario esperado por cada año adicional de educación, manteniendo lo demás constante?

Todo modelo econométrico parte de una especificación teórica, se estima con datos y se somete a diagnóstico y validación. Este documento recorre ese camino completo: comenzamos por las herramientas matemáticas que sostienen la estimación (álgebra lineal), pasamos por los fundamentos probabilísticos, desarrollamos el modelo de regresión lineal simple y múltiple, estudiamos las propiedades de sus estimadores, la inferencia estadística asociada y el diagnóstico de los supuestos clásicos, cerrando con variables cualitativas y ejercicios integradores.

1.2 ¿Por qué R?

R es el lenguaje estándar en la enseñanza y práctica de la econometría académica por varias razones: es de código abierto, tiene una sintaxis matricial nativa (ideal para álgebra lineal), cuenta con un ecosistema estadístico maduro (lm(), car, lmtest, sandwich, entre otros) y permite reproducibilidad completa: cualquier lector puede ejecutar el mismo código y obtener exactamente los mismos resultados.

1.3 Datos reales utilizados en este documento

Para que los ejemplos sean lo más realistas posible, este documento usa exclusivamente conjuntos de datos reales, ya incluidos en R o en paquetes estándar de uso académico (ninguno es simulado):

Dataset Contenido Fuente / época
cars Velocidad y distancia de frenado de automóviles Ezekiel (1930), EE.UU.
swiss Fecundidad y determinantes socioeconómicos de 47 provincias Suiza, 1888 (Oficina de Estadística suiza)
longley PIB, empleo, fuerza armada, población EE.UU., 1947-1962 (Bureau of Labor Statistics)
state.x77 Ingreso, criminalidad, esperanza de vida, educación 50 estados de EE.UU., ~1970 (U.S. Census)
Boston (paquete MASS) Precios de vivienda y determinantes urbanos Boston, EE.UU., 1970 (Harrison & Rubinfeld, 1978)
mtcars Consumo de combustible y características de automóviles Motor Trend US, 1974

1.4 Paquetes utilizados

install.packages(c("MASS", "car", "lmtest", "sandwich"))
library(MASS)      # dataset Boston
library(car)       # vif() - factor de inflacion de varianza
library(lmtest)    # bptest(), dwtest() - tests de diagnostico
library(sandwich)  # errores estandar robustos a heterocedasticidad

2 Álgebra lineal aplicada a la econometría

El modelo de regresión lineal múltiple se escribe y se resuelve enteramente en notación matricial. Antes de llegar ahí, repasamos las operaciones de álgebra lineal que R maneja de forma nativa y que sostienen toda la estimación econométrica posterior.

2.1 Escalares, vectores y matrices en R

Un escalar es un número; un vector es un arreglo unidimensional de números; una matriz es un arreglo bidimensional (filas \(\times\) columnas). En notación económica, si observamos \(n\) individuos y \(k\) variables explicativas, organizamos los datos en una matriz \(X\) de dimensión \(n\times k\).

# Vector: ingresos mensuales (S/.) de 5 personas (dato ilustrativo)
ingreso <- c(1200, 1800, 950, 2500, 1600)
educacion <- c(11, 16, 8, 18, 14)   # anios de educacion

ingreso
## [1] 1200 1800  950 2500 1600
class(ingreso)
## [1] "numeric"
length(ingreso)
## [1] 5
# Matriz: combinamos ambos vectores como columnas
X <- cbind(intercepto = 1, educacion = educacion)
X
##      intercepto educacion
## [1,]          1        11
## [2,]          1        16
## [3,]          1         8
## [4,]          1        18
## [5,]          1        14
dim(X)      # 5 filas (observaciones), 2 columnas (variables)
## [1] 5 2

2.2 Operaciones matriciales fundamentales

Las operaciones clave que usaremos repetidamente son la transposición (\(X'\) o \(X^\top\)), el producto matricial (\(AB\), distinto de la multiplicación elemento a elemento) y la suma.

A <- matrix(c(2, 1, 0, 3), nrow = 2, byrow = TRUE)
B <- matrix(c(1, 4, 2, 1), nrow = 2, byrow = TRUE)
A; B
##      [,1] [,2]
## [1,]    2    1
## [2,]    0    3
##      [,1] [,2]
## [1,]    1    4
## [2,]    2    1
A + B                 # suma elemento a elemento
##      [,1] [,2]
## [1,]    3    5
## [2,]    2    4
A %*% B                # producto MATRICIAL (ojo: no es A*B)
##      [,1] [,2]
## [1,]    4    9
## [2,]    6    3
A * B                  # producto elemento a elemento (Hadamard) - distinto de %*%
##      [,1] [,2]
## [1,]    2    4
## [2,]    0    3
t(A)                   # transpuesta de A
##      [,1] [,2]
## [1,]    2    0
## [2,]    1    3

Nota importante: en R, A %*% B es el producto matricial (regla fila-por-columna), mientras que A * B multiplica elemento a elemento. Confundir ambos operadores es el error más común al traducir álgebra lineal a código, y en econometría lleva a resultados completamente incorrectos en la estimación de \(\hat\beta\).

2.3 Determinante, rango y matriz inversa

El determinante de una matriz cuadrada indica si es invertible (\(\det(A)\neq 0\)) y está ligado geométricamente al “volumen” que genera la transformación lineal. El rango es el número de filas (o columnas) linealmente independientes. La matriz inversa \(A^{-1}\) satisface \(AA^{-1}=I\), y es la pieza central del estimador MCO: \(\hat\beta = (X'X)^{-1}X'y\).

M <- matrix(c(4, 2, 7, 6), nrow = 2, byrow = TRUE)
M
##      [,1] [,2]
## [1,]    4    2
## [2,]    7    6
det(M)                  # determinante
## [1] 10
solve(M)                # matriz inversa
##      [,1] [,2]
## [1,]  0.6 -0.2
## [2,] -0.7  0.4
M %*% solve(M)           # verificacion: M %*% M^-1 = I (identidad)
##      [,1]          [,2]
## [1,]    1 -1.110223e-16
## [2,]    0  1.000000e+00
qr(M)$rank               # rango de la matriz
## [1] 2

Cuando el determinante es exactamente cero (o muy cercano a cero), la matriz es singular (no invertible) o casi singular. En econometría esto ocurre cuando existe multicolinealidad perfecta entre las variables explicativas -por ejemplo, si una columna de \(X\) es combinación lineal exacta de otra-, y hace que \((X'X)^{-1}\) no exista, impidiendo estimar \(\hat\beta\) por MCO. Lo verificamos en la práctica más adelante (Sección 8.1) con el dataset longley.

# Ejemplo de matriz singular: la segunda columna es 2 veces la primera
S <- matrix(c(1, 2, 3, 6), nrow = 2, byrow = TRUE)
S
##      [,1] [,2]
## [1,]    1    2
## [2,]    3    6
det(S)                     # exactamente 0 -> matriz singular
## [1] 0
tryCatch(solve(S), error = function(e) cat("Error al invertir:", conditionMessage(e), "\n"))
## Error al invertir: Lapack routine dgesv: system is exactly singular: U[2,2] = 0

2.4 Sistemas de ecuaciones lineales (\(Ax = b\))

Muchos problemas económicos se reducen a resolver un sistema lineal \(Ax=b\). Un ejemplo clásico es encontrar el equilibrio de mercado: dadas una función de demanda \(Q^d = a - bP\) y una función de oferta \(Q^s = c + dP\), el equilibrio (\(Q^d=Q^s\)) es la solución de un sistema lineal en \((P^*, Q^*)\).

# Demanda: Q = 100 - 2P  ->  Q + 2P = 100
# Oferta:  Q = -20 + 3P  ->  Q - 3P = -20
A_mercado <- matrix(c(1, 2,
                       1, -3), nrow = 2, byrow = TRUE)
b_mercado <- c(100, -20)

equilibrio <- solve(A_mercado, b_mercado)
names(equilibrio) <- c("Q*", "P*")
equilibrio
## Q* P* 
## 52 24

El equilibrio de mercado es \(P^*=\) 24 y \(Q^*=\) 52. Esta misma lógica de resolver \(Ax=b\) es exactamente la que R utiliza internamente al estimar \(\hat\beta=(X'X)^{-1}X'y\) en la regresión lineal.

2.5 Autovalores y autovectores

Un autovalor \(\lambda\) y su autovector asociado \(v\) de una matriz \(A\) satisfacen \(Av = \lambda v\): la transformación lineal solo “escala” al vector \(v\), sin rotarlo. En econometría, los autovalores de \(X'X\) son la base del Análisis de Componentes Principales (PCA) y del diagnóstico de multicolinealidad (número de condición).

Sigma <- matrix(c(4, 2,
                   2, 3), nrow = 2, byrow = TRUE)   # matriz simetrica (p.ej., de covarianzas)
ev <- eigen(Sigma)
ev$values      # autovalores
## [1] 5.561553 1.438447
ev$vectors     # autovectores (columnas)
##            [,1]       [,2]
## [1,] -0.7882054  0.6154122
## [2,] -0.6154122 -0.7882054
# Verificacion: Sigma %*% v = lambda * v para el primer autovector
Sigma %*% ev$vectors[,1]
##           [,1]
## [1,] -4.383646
## [2,] -3.422648
ev$values[1] * ev$vectors[,1]
## [1] -4.383646 -3.422648

Ambos resultados coinciden, confirmando la definición \(Av=\lambda v\). El número de condición de una matriz, \(\kappa(A) = \lambda_{max}/\lambda_{min}\), es un diagnóstico directo de multicolinealidad: valores muy grandes (\(\kappa>30\) es una regla empírica común) señalan que \(X'X\) está cerca de ser singular.

kappa_Sigma <- max(ev$values) / min(ev$values)
kappa_Sigma
## [1] 3.866359

2.6 Formas cuadráticas y matrices definidas positivas

Una forma cuadrática es una expresión \(x'Ax\). Una matriz simétrica \(A\) es definida positiva si \(x'Ax>0\) para todo \(x\neq 0\) -equivalentemente, si todos sus autovalores son positivos. Esta propiedad es la que garantiza que la matriz de varianzas-covarianzas de un vector aleatorio sea válida, y que \((X'X)\) sea invertible cuando las columnas de \(X\) son linealmente independientes.

# Todos los autovalores de Sigma son positivos -> definida positiva
all(ev$values > 0)
## [1] TRUE
# La matriz X'X (con X del ejemplo inicial) tambien debe ser definida positiva
XtX <- t(X) %*% X
XtX
##            intercepto educacion
## intercepto          5        67
## educacion          67       961
eigen(XtX)$values
## [1] 965.672767   0.327233

2.7 Ejercicios de álgebra lineal

  1. (Operaciones básicas) Defina en R las matrices \(A=\begin{pmatrix}3&1\\2&4\end{pmatrix}\) y \(B=\begin{pmatrix}0&2\\1&1\end{pmatrix}\). Calcule \(A+B\), \(AB\) (producto matricial), \(A'\) y \(\det(A)\).

  2. (Sistema de ecuaciones) Una economía cerrada tiene dos sectores. El sector 1 necesita 0.2 unidades de su propio producto y 0.3 unidades del sector 2 por cada unidad producida; el sector 2 necesita 0.4 del sector 1 y 0.1 de sí mismo. Si la demanda final es \(d=(100, 50)'\), resuelva el modelo insumo-producto de Leontief \(x = (I-A)^{-1}d\) para hallar la producción total \(x\) de cada sector.

  3. (Autovalores) Calcule los autovalores de \(\begin{pmatrix}5&0\\0&2\end{pmatrix}\) sin usar eigen() (a mano, usando \(\det(A-\lambda I)=0\)) y luego verifique su resultado con eigen() en R.

  4. (Rango y singularidad) Construya una matriz \(3\times 3\) donde la tercera columna sea la suma de las dos primeras. Verifique con qr()$rank que su rango es 2, no 3, y explique la implicancia económica si esas columnas fueran variables explicativas de una regresión.

Ver soluciones 
# 1.
A <- matrix(c(3,1,2,4), nrow=2, byrow=TRUE)
B <- matrix(c(0,2,1,1), nrow=2, byrow=TRUE)
A + B
##      [,1] [,2]
## [1,]    3    3
## [2,]    3    5
A %*% B
##      [,1] [,2]
## [1,]    1    7
## [2,]    4    8
t(A)
##      [,1] [,2]
## [1,]    3    2
## [2,]    1    4
det(A)
## [1] 10
# 2. Modelo de Leontief
A_leontief <- matrix(c(0.2, 0.4,
                        0.3, 0.1), nrow = 2, byrow = TRUE)
d <- c(100, 50)
I2 <- diag(2)
x_total <- solve(I2 - A_leontief) %*% d
x_total
##          [,1]
## [1,] 183.3333
## [2,] 116.6667
# 3. Autovalores de una matriz diagonal = sus elementos diagonales: 5 y 2
eigen(matrix(c(5,0,0,2), nrow=2))$values
## [1] 5 2
# 4.
C <- matrix(c(1,2,3, 4,5,9, 7,8,15), nrow=3, byrow=TRUE)  # col3 = col1+col2
qr(C)$rank   # rango 2: hay multicolinealidad perfecta entre las 3 columnas
## [1] 2
Si las tres columnas fueran variables explicativas de una regresión, \(X'X\) sería singular y lm() eliminaría automáticamente una variable (o R devolvería coeficientes NA), porque no hay suficiente variación independiente para separar el efecto de cada una.

3 Fundamentos de probabilidad y estadística en R

3.1 Variables aleatorias y distribuciones

Una variable aleatoria asigna un número a cada resultado de un experimento. En econometría trabajamos constantemente con la distribución Normal (para errores y estimadores en muestras grandes), la t de Student (inferencia con varianza desconocida), la Chi-cuadrado y la F (contrastes de hipótesis conjuntas).

estilo_base(mfrow = c(1,2), mar = c(4.2, 4.2, 3, 1))
x <- seq(-4, 4, 0.01)
plot(x, dnorm(x), type = "l", col = col_azul, lwd = 2,
     main = "Densidad Normal estandar", xlab = "x", ylab = "f(x)")
plot(x, dt(x, df = 5), type = "l", col = col_rojo, lwd = 2,
     main = "Densidad t-Student (5 g.l.)", xlab = "x", ylab = "f(x)")

par(mfrow = c(1,1))

3.2 Momentos: media, varianza, covarianza, correlación

# Usamos datos reales: swiss (Suiza, 1888)
data(swiss)
head(swiss, 3)
##              Fertility Agriculture Examination Education Catholic
## Courtelary        80.2        17.0          15        12     9.96
## Delemont          83.1        45.1           6         9    84.84
## Franches-Mnt      92.5        39.7           5         5    93.40
##              Infant.Mortality
## Courtelary               22.2
## Delemont                 22.2
## Franches-Mnt             20.2
mean(swiss$Fertility)          # media
## [1] 70.14255
var(swiss$Fertility)           # varianza muestral
## [1] 156.0425
sd(swiss$Fertility)            # desviacion estandar
## [1] 12.4917
cov(swiss$Fertility, swiss$Education)     # covarianza
## [1] -79.72951
cor(swiss$Fertility, swiss$Education)     # correlacion de Pearson
## [1] -0.6637889

La correlación entre fecundidad y educación es -0.664: negativa y moderadamente fuerte, consistente con la teoría demográfica clásica (a mayor nivel educativo, menor tasa de fecundidad).

M_cor <- cor(swiss)
round(M_cor, 2)
##                  Fertility Agriculture Examination Education Catholic
## Fertility             1.00        0.35       -0.65     -0.66     0.46
## Agriculture           0.35        1.00       -0.69     -0.64     0.40
## Examination          -0.65       -0.69        1.00      0.70    -0.57
## Education            -0.66       -0.64        0.70      1.00    -0.15
## Catholic              0.46        0.40       -0.57     -0.15     1.00
## Infant.Mortality      0.42       -0.06       -0.11     -0.10     0.18
##                  Infant.Mortality
## Fertility                    0.42
## Agriculture                 -0.06
## Examination                 -0.11
## Education                   -0.10
## Catholic                     0.18
## Infant.Mortality             1.00

3.3 Distribuciones muestrales: t, F y Chi-cuadrado en la práctica

qt(0.975, df = 44)      # cuantil t al 97.5% con 44 g.l. (usado en IC al 95%)
## [1] 2.015368
qf(0.95, df1 = 2, df2 = 44)   # cuantil F al 95%
## [1] 3.209278
qchisq(0.95, df = 5)          # cuantil chi-cuadrado al 95%
## [1] 11.0705

3.4 Ley de los Grandes Números y Teorema Central del Límite

El Teorema Central del Límite (TCL) es la razón por la cual, incluso si los errores individuales no son normales, la distribución del estimador MCO en muestras grandes se aproxima a una Normal. Lo verificamos con una simulación: promediamos muestras de una distribución claramente no-normal (uniforme) y observamos que la distribución de esos promedios sí se vuelve Normal.

set.seed(2025)
medias_muestrales <- replicate(5000, mean(runif(30, min = 0, max = 10)))

estilo_base()
hist(medias_muestrales, breaks = 40, col = col_fill, border = "white",
     main = "TCL: distribucion de medias muestrales (n=30, dist. Uniforme)",
     xlab = "Media muestral", ylab = "Frecuencia")
curve(dnorm(x, mean = mean(medias_muestrales), sd = sd(medias_muestrales)) *
        5000 * diff(hist(medias_muestrales, plot = FALSE, breaks = 40)$breaks)[1],
      add = TRUE, col = col_rojo, lwd = 2)

Aunque la variable original es Uniforme (plana, nada parecida a una campana), la distribución de las 5000 medias muestrales sí sigue de cerca la curva Normal (línea roja) superpuesta. Este es el TCL en acción, y es el fundamento que permite construir intervalos de confianza y pruebas t sobre \(\hat\beta\) incluso sin normalidad exacta de los errores.

3.5 Ejercicios

  1. Usando los datos swiss, calcule la correlación entre Agriculture y Fertility. Interprete el signo económico.
  2. Simule 10 000 lanzamientos de un dado (sample(1:6, 10000, replace=TRUE)) y calcule la media y varianza muestral. Compárelas con los valores teóricos (\(E[X]=3.5\), \(Var(X)=35/12\)).
  3. Repita la simulación del TCL de esta sección pero partiendo de una distribución Exponencial (rexp()) en vez de Uniforme. ¿Sigue funcionando el TCL?
Ver soluciones 
# 1.
cor(swiss$Agriculture, swiss$Fertility)
## [1] 0.3530792
# Positiva: provincias con mayor proporcion de hombres en agricultura
# tienden a tener mayor fecundidad (asociado a menor industrializacion/urbanizacion).

# 2.
set.seed(1)
dado <- sample(1:6, 10000, replace = TRUE)
mean(dado); 3.5
## [1] 3.5138
## [1] 3.5
var(dado); 35/12
## [1] 2.926902
## [1] 2.916667
# 3.
medias_exp <- replicate(5000, mean(rexp(30, rate = 1)))
estilo_base()
hist(medias_exp, breaks = 40, col = col_fill_rojo, border = "white",
     main = "TCL con distribucion Exponencial subyacente", xlab = "Media muestral",
     ylab = "Frecuencia")

# Si, el TCL sigue funcionando: la distribucion de las medias se aproxima a
# una Normal incluso partiendo de una Exponencial (fuertemente asimetrica).

4 El modelo de regresión lineal simple

4.1 Especificación del modelo

El modelo de regresión lineal simple relaciona una variable dependiente \(Y\) con una única variable explicativa \(X\): \[Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i, \qquad i=1,\dots,n,\] donde \(\beta_0\) es el intercepto, \(\beta_1\) es la pendiente (el efecto marginal de \(X\) sobre \(Y\)) y \(\varepsilon_i\) es un término de error no observado que captura todo lo que influye en \(Y_i\) además de \(X_i\).

4.2 Estimación por Mínimos Cuadrados Ordinarios (MCO)

El método MCO elige \(\hat\beta_0,\hat\beta_1\) que minimizan la suma de residuos al cuadrado, \(\sum_i \hat\varepsilon_i^2 = \sum_i (Y_i-\hat\beta_0-\hat\beta_1 X_i)^2\). Derivando e igualando a cero se obtienen las fórmulas cerradas: \[\hat\beta_1 = \frac{\sum_i (X_i-\bar X)(Y_i-\bar Y)}{\sum_i (X_i-\bar X)^2} = \frac{\operatorname{Cov}(X,Y)}{\operatorname{Var}(X)}, \qquad \hat\beta_0 = \bar Y - \hat\beta_1 \bar X.\]

4.3 Aplicación con datos reales: cars

El dataset cars (Ezekiel, 1930) contiene la velocidad de 50 automóviles (mph) y su distancia de frenado (pies) -un ejemplo real y clásico de relación física con componente aleatorio.

data(cars)
head(cars)
##   speed dist
## 1     4    2
## 2     4   10
## 3     7    4
## 4     7   22
## 5     8   16
## 6     9   10
# Estimacion manual con las formulas cerradas
x <- cars$speed; y <- cars$dist
beta1_hat <- cov(x, y) / var(x)
beta0_hat <- mean(y) - beta1_hat * mean(x)
c(beta0 = beta0_hat, beta1 = beta1_hat)
##      beta0      beta1 
## -17.579095   3.932409
# Verificacion con lm()
modelo_simple <- lm(dist ~ speed, data = cars)
coef(modelo_simple)
## (Intercept)       speed 
##  -17.579095    3.932409

Las dos estimaciones coinciden exactamente (3.932), confirmando que lm() internamente resuelve el mismo problema de minimización.

estilo_base()
plot(cars$speed, cars$dist, pch = 19, col = col_gris,
     xlab = "Velocidad (mph)", ylab = "Distancia de frenado (pies)",
     main = "Regresion simple: velocidad vs. distancia de frenado")
abline(modelo_simple, col = col_rojo, lwd = 2.4)
legend("topleft", legend = c("Datos observados","Recta MCO ajustada"),
       col = c(col_gris, col_rojo), pch = c(19, NA), lty = c(NA,1), lwd = c(NA,2.4),
       bty = "n", bg = "white", inset = 0.02, cex = 0.9)

Interpretación: por cada milla por hora adicional de velocidad, la distancia de frenado esperada aumenta en 3.93 pies, en promedio.

4.4 Bondad de ajuste: \(R^2\)

\[R^2 = 1 - \frac{\sum_i \hat\varepsilon_i^2}{\sum_i (Y_i-\bar Y)^2} = 1 - \frac{SCR}{SCT}\] mide la proporción de la variabilidad de \(Y\) explicada por el modelo.

resumen_simple <- summary(modelo_simple)
resumen_simple$r.squared
## [1] 0.6510794
# Verificacion manual
y_hat <- fitted(modelo_simple)
SCR <- sum((y - y_hat)^2)
SCT <- sum((y - mean(y))^2)
1 - SCR/SCT
## [1] 0.6510794

4.5 Ejercicios

  1. Usando cars, estime a mano (con las fórmulas cerradas) un modelo que explique dist en función de speed pero SIN intercepto (\(Y_i=\beta_1 X_i+\varepsilon_i\)). Pista: \(\hat\beta_1=\sum X_iY_i/\sum X_i^2\).
  2. Grafique los residuos (\(\hat\varepsilon_i\)) contra speed. ¿Observa algún patrón que sugiera no linealidad?
  3. Calcule el \(R^2\) de un modelo dist ~ speed usando solo la correlación: \(R^2 = \text{cor}(X,Y)^2\) en la regresión simple. Verifíquelo.
Ver soluciones 
# 1.
beta1_sin_intercepto <- sum(x*y) / sum(x^2)
beta1_sin_intercepto
## [1] 2.909132
lm(dist ~ speed - 1, data = cars)   # verificacion
## 
## Call:
## lm(formula = dist ~ speed - 1, data = cars)
## 
## Coefficients:
## speed  
## 2.909
# 2.
estilo_base()
plot(cars$speed, resid(modelo_simple), pch = 19, col = col_azul,
     xlab = "Velocidad (mph)", ylab = "Residuo",
     main = "Residuos vs. velocidad")
abline(h = 0, lty = 2, col = col_gris)

# Se aprecia un leve patron: los residuos tienden a ser mas dispersos
# (mayor varianza) para velocidades altas -> posible heterocedasticidad.

# 3.
cor(x,y)^2
## [1] 0.6510794
resumen_simple$r.squared    # coinciden
## [1] 0.6510794

5 El modelo de regresión lineal múltiple (enfoque matricial)

5.1 Especificación matricial

Cuando hay \(k\) variables explicativas, el modelo se escribe de forma compacta como \[y = X\beta + \varepsilon,\] donde \(y\) es el vector \(n\times 1\) de la variable dependiente, \(X\) es la matriz \(n\times(k+1)\) de regresores (incluyendo una columna de unos para el intercepto), \(\beta\) es el vector \((k+1)\times 1\) de coeficientes, y \(\varepsilon\) es el vector \(n\times 1\) de errores.

5.2 El estimador MCO matricial

Minimizando la suma de cuadrados \(\varepsilon'\varepsilon = (y-X\beta)'(y-X\beta)\) respecto a \(\beta\) se obtiene la condición de primer orden \(X'(y-X\beta)=0\), cuya solución es: \[\boxed{\hat\beta = (X'X)^{-1}X'y}\] Esta es exactamente la fórmula que programamos con solve() en la Sección 2.4 para resolver sistemas lineales -la regresión múltiple es, en esencia, resolver el sistema \(X'X\hat\beta = X'y\).

5.3 Aplicación con datos reales: swiss (Suiza, 1888)

Modelamos la fecundidad (Fertility) en función de la proporción de hombres en agricultura, el nivel educativo y la proporción católica.

data(swiss)
str(swiss)
## 'data.frame':    47 obs. of  6 variables:
##  $ Fertility       : num  80.2 83.1 92.5 85.8 76.9 76.1 83.8 92.4 82.4 82.9 ...
##  $ Agriculture     : num  17 45.1 39.7 36.5 43.5 35.3 70.2 67.8 53.3 45.2 ...
##  $ Examination     : int  15 6 5 12 17 9 16 14 12 16 ...
##  $ Education       : int  12 9 5 7 15 7 7 8 7 13 ...
##  $ Catholic        : num  9.96 84.84 93.4 33.77 5.16 ...
##  $ Infant.Mortality: num  22.2 22.2 20.2 20.3 20.6 26.6 23.6 24.9 21 24.4 ...
y <- swiss$Fertility
X <- as.matrix(cbind(Intercepto = 1, swiss[, c("Agriculture","Education","Catholic")]))
head(X)
##              Intercepto Agriculture Education Catholic
## Courtelary            1        17.0        12     9.96
## Delemont              1        45.1         9    84.84
## Franches-Mnt          1        39.7         5    93.40
## Moutier               1        36.5         7    33.77
## Neuveville            1        43.5        15     5.16
## Porrentruy            1        35.3         7    90.57
# Estimacion manual: beta_hat = (X'X)^-1 X'y
XtX <- t(X) %*% X
Xty <- t(X) %*% y
beta_hat_manual <- solve(XtX) %*% Xty
beta_hat_manual
##                   [,1]
## Intercepto  86.2250198
## Agriculture -0.2030377
## Education   -1.0721468
## Catholic     0.1452013
# Verificacion con lm()
modelo_multiple <- lm(Fertility ~ Agriculture + Education + Catholic, data = swiss)
coef(modelo_multiple)
## (Intercept) Agriculture   Education    Catholic 
##  86.2250198  -0.2030377  -1.0721468   0.1452013

Ambas estimaciones -la manual vía álgebra matricial y la de lm()- coinciden exactamente, confirmando que lm() resuelve internamente \(\hat\beta=(X'X)^{-1}X'y\).

Interpretación económica: manteniendo todo lo demás constante, un aumento de 1 punto porcentual en la proporción de hombres en agricultura se asocia a un aumento de -0.203 puntos en la tasa de fecundidad; un año adicional de educación promedio se asocia a una caída de 1.072 puntos.

5.4 Ejercicios

  1. Agregue la variable Examination (porcentaje de reclutas con nota alta en el examen militar) al modelo swiss y reestime manualmente \(\hat\beta=(X'X)^{-1}X'y\). Verifique con lm().
  2. Calcule los valores ajustados \(\hat y = X\hat\beta\) y los residuos \(\hat\varepsilon = y-\hat y\) usando solo álgebra matricial (sin fitted() ni resid()), y confirme que coinciden con los de lm().
  3. Demuestre en R que \(X'\hat\varepsilon = 0\) (los residuos son ortogonales a cada regresor) -una propiedad algebraica directa de la condición de primer orden del MCO.
Ver soluciones 
# 1.
X2 <- as.matrix(cbind(Intercepto=1, swiss[,c("Agriculture","Education","Catholic","Examination")]))
beta_hat2 <- solve(t(X2)%*%X2) %*% (t(X2)%*%y)
beta_hat2
##                   [,1]
## Intercepto  91.0554239
## Agriculture -0.2206455
## Education   -0.9616124
## Catholic     0.1244184
## Examination -0.2605824
coef(lm(Fertility ~ Agriculture + Education + Catholic + Examination, data = swiss))
## (Intercept) Agriculture   Education    Catholic Examination 
##  91.0554239  -0.2206455  -0.9616124   0.1244184  -0.2605824
# 2.
y_hat_manual <- X %*% beta_hat_manual
resid_manual <- y - y_hat_manual
head(cbind(manual = y_hat_manual, lm = fitted(modelo_multiple)))
##                             lm
## Courtelary   71.35382 71.35382
## Delemont     79.73758 79.73758
## Franches-Mnt 86.36550 86.36550
## Moutier      76.21257 76.21257
## Neuveville   62.05992 62.05992
## Porrentruy   84.70365 84.70365
max(abs(resid_manual - resid(modelo_multiple)))   # diferencia practicamente nula
## [1] 6.010747e-13
# 3.
round(t(X) %*% resid_manual, 8)   # vector (casi) cero para cada columna de X
##             [,1]
## Intercepto     0
## Agriculture    0
## Education      0
## Catholic       0

6 Propiedades del estimador MCO

6.1 Los supuestos de Gauss-Markov

Bajo los siguientes supuestos clásicos, el estimador MCO es el mejor estimador lineal insesgado (BLUE, Best Linear Unbiased Estimator):

  1. Linealidad en los parámetros: \(y=X\beta+\varepsilon\).
  2. Muestreo aleatorio de la población.
  3. No colinealidad perfecta entre regresores (rango completo de \(X\)).
  4. Exogeneidad estricta: \(E[\varepsilon_i \mid X] = 0\).
  5. Homocedasticidad: \(\operatorname{Var}(\varepsilon_i\mid X)=\sigma^2\) (constante).
  6. No autocorrelación: \(\operatorname{Cov}(\varepsilon_i,\varepsilon_j\mid X)=0\) para \(i\neq j\).

6.2 Insesgadez: \(E[\hat\beta]=\beta\)

Bajo los supuestos 1-4, \(\hat\beta\) es insesgado: en promedio, sobre repetidas muestras, el estimador acierta el valor verdadero de \(\beta\). Lo verificamos con una simulación Monte Carlo: generamos 2000 muestras artificiales de un modelo con \(\beta_0=5\), \(\beta_1=2\) conocidos, estimamos \(\hat\beta\) en cada una y promediamos.

set.seed(2025)
beta0_verdadero <- 5; beta1_verdadero <- 2
n <- 100
x_mc <- runif(n, 0, 10)

simular_beta1 <- function(){
  e <- rnorm(n, sd = 3)
  y_mc <- beta0_verdadero + beta1_verdadero * x_mc + e
  coef(lm(y_mc ~ x_mc))[2]
}

beta1_simulados <- replicate(2000, simular_beta1())

mean(beta1_simulados)     # deberia acercarse a beta1_verdadero = 2
## [1] 1.997347
beta1_verdadero
## [1] 2
estilo_base()
hist(beta1_simulados, breaks = 40, col = col_fill, border = "white",
     main = "Distribucion muestral de beta1_hat (2000 simulaciones Monte Carlo)",
     xlab = expression(hat(beta)[1]), ylab = "Frecuencia")
abline(v = beta1_verdadero, col = col_rojo, lwd = 2.5)
abline(v = mean(beta1_simulados), col = col_azul, lwd = 2, lty = 2)
legend("topright", legend = c("Valor verdadero (beta1=2)","Promedio de las estimaciones"),
       col = c(col_rojo, col_azul), lwd = c(2.5,2), lty = c(1,2), bty = "n",
       bg = "white", inset = 0.02, cex = 0.85)

El promedio de las 2000 estimaciones (1.9973) está prácticamente encima del valor verdadero (\(\beta_1=2\)), ilustrando la insesgadez de forma empírica.

6.3 Eficiencia: menor varianza entre los estimadores lineales insesgados

El teorema de Gauss-Markov garantiza que, entre todos los estimadores lineales e insesgados de \(\beta\), el MCO tiene la menor varianza. Lo comparamos con un estimador alternativo simple (la pendiente entre los dos puntos extremos de \(X\)), también insesgado pero menos eficiente.

estimador_extremos <- function(){
  e <- rnorm(n, sd = 3)
  y_mc <- beta0_verdadero + beta1_verdadero * x_mc + e
  i_min <- which.min(x_mc); i_max <- which.max(x_mc)
  (y_mc[i_max]-y_mc[i_min]) / (x_mc[i_max]-x_mc[i_min])
}

beta1_extremos <- replicate(2000, estimador_extremos())

cat("Varianza del estimador MCO:      ", round(var(beta1_simulados),5), "\n")
## Varianza del estimador MCO:       0.01084
cat("Varianza del estimador extremos: ", round(var(beta1_extremos),5), "\n")
## Varianza del estimador extremos:  0.19412

Ambos estimadores son insesgados (sus promedios se acercan a 2), pero el estimador MCO tiene varianza notablemente menor -es decir, es más preciso- confirmando en la práctica el teorema de Gauss-Markov.

6.4 Ejercicios

  1. Repita la simulación de insesgadez pero con errores heterocedásticos (\(\text{sd}=0.5\cdot x_{mc}\) en vez de constante). ¿Sigue siendo insesgado \(\hat\beta_1\)? ¿Qué esperaría que pase con su varianza?
  2. Aumente el tamaño muestral de \(n=100\) a \(n=1000\) en la simulación Monte Carlo original y observe cómo se reduce la varianza de \(\hat\beta_1\) (consistencia).
Ver soluciones 
# 1.
simular_beta1_heter <- function(){
  e <- rnorm(n, sd = 0.5*x_mc)   # varianza del error depende de x -> heterocedasticidad
  y_mc <- beta0_verdadero + beta1_verdadero * x_mc + e
  coef(lm(y_mc ~ x_mc))[2]
}
beta1_heter <- replicate(2000, simular_beta1_heter())
mean(beta1_heter)   # sigue siendo insesgado (la heterocedasticidad no rompe insesgadez)
## [1] 1.999201
var(beta1_heter)    # la varianza cambia; MCO ya no es el mas eficiente (deja de ser BLUE)
## [1] 0.01169596
# 2.
n_grande <- 1000
x_grande <- runif(n_grande, 0, 10)
simular_grande <- function(){
  e <- rnorm(n_grande, sd = 3)
  y_g <- beta0_verdadero + beta1_verdadero * x_grande + e
  coef(lm(y_g ~ x_grande))[2]
}
beta1_grande <- replicate(1000, simular_grande())
var(beta1_grande)   # notablemente menor que con n=100: el estimador es consistente
## [1] 0.001103967

7 Inferencia estadística en el modelo de regresión

7.1 Errores estándar y la prueba t individual

Cada coeficiente estimado \(\hat\beta_j\) tiene asociado un error estándar \(\text{ee}(\hat\beta_j)\) que mide su precisión. La prueba t contrasta \(H_0:\beta_j=0\) (la variable \(X_j\) no tiene efecto) usando el estadístico \[t = \frac{\hat\beta_j - 0}{\text{ee}(\hat\beta_j)} \sim t_{n-k-1} \text{ bajo } H_0.\]

7.2 Aplicación con datos reales: Boston (precios de vivienda)

data(Boston)
str(Boston[,1:6])
## 'data.frame':    506 obs. of  6 variables:
##  $ crim : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn   : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus: num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox  : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm   : num  6.58 6.42 7.18 7 7.15 ...
modelo_boston <- lm(medv ~ crim + rm + age + dis + tax, data = Boston)
summary(modelo_boston)
## 
## Call:
## lm(formula = medv ~ crim + rm + age + dis + tax, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.614  -2.911  -0.833   1.987  40.959 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -11.796901   3.236248  -3.645 0.000295 ***
## crim         -0.140933   0.037531  -3.755 0.000194 ***
## rm            7.731058   0.390879  19.779  < 2e-16 ***
## age          -0.079425   0.014272  -5.565 4.28e-08 ***
## dis          -0.944654   0.194106  -4.867 1.52e-06 ***
## tax          -0.011553   0.002144  -5.389 1.09e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.856 on 500 degrees of freedom
## Multiple R-squared:  0.5986, Adjusted R-squared:  0.5946 
## F-statistic: 149.2 on 5 and 500 DF,  p-value: < 2.2e-16

La variable medv es el valor mediano de la vivienda (en miles de USD, 1970), crim la tasa de criminalidad, rm el número medio de habitaciones, age la proporción de viviendas antiguas, dis la distancia a centros de empleo y tax la tasa impositiva.

Lectura de la salida: cada fila del bloque Coefficients reporta \(\hat\beta_j\), su error estándar, el estadístico t y el p-valor. Un p-valor menor a 0.05 indica que rechazamos \(H_0:\beta_j=0\) al 5% de significancia.

7.3 La prueba F de significancia global

La prueba F contrasta la hipótesis conjunta \(H_0:\beta_1=\beta_2=\dots=\beta_k=0\) (ninguna variable explica a \(Y\)) frente a \(H_1\): al menos una sí. El estadístico es \[F = \frac{(SCT-SCR)/k}{SCR/(n-k-1)} \sim F_{k,\,n-k-1} \text{ bajo } H_0.\]

resumen_boston <- summary(modelo_boston)
resumen_boston$fstatistic     # valor F, g.l. numerador, g.l. denominador
##    value    numdf    dendf 
## 149.1559   5.0000 500.0000
# Verificacion manual
y_b <- Boston$medv
SCT <- sum((y_b - mean(y_b))^2)
SCR <- sum(resid(modelo_boston)^2)
k <- 5; n_b <- nrow(Boston)
F_manual <- ((SCT-SCR)/k) / (SCR/(n_b-k-1))
F_manual
## [1] 149.1559

El estadístico F es enorme (muy por encima del valor crítico), por lo que rechazamos rotundamente \(H_0\): el conjunto de variables sí explica significativamente el precio de la vivienda.

7.4 Intervalos de confianza

confint(modelo_boston, level = 0.95)
##                    2.5 %       97.5 %
## (Intercept) -18.15522239 -5.438580527
## crim         -0.21467159 -0.067195121
## rm            6.96309004  8.499025884
## age          -0.10746483 -0.051384372
## dis          -1.32601674 -0.563290569
## tax          -0.01576487 -0.007341134

Un intervalo de confianza al 95% que no contiene el cero indica que el coeficiente es significativo a ese nivel -consistente con la prueba t.

estilo_base()
coefs <- coef(modelo_boston)[-1]     # sin intercepto
ic <- confint(modelo_boston)[-1,]
ord <- order(coefs)

plot(coefs[ord], 1:length(coefs), pch = 19, col = col_azul, cex=1.3,
     xlim = range(ic), yaxt = "n", ylab = "",
     xlab = "Coeficiente estimado (IC 95%)",
     main = "Coeficientes e intervalos de confianza - modelo Boston")
axis(2, at = 1:length(coefs), labels = names(coefs)[ord], las = 1, cex.axis = 0.85)
segments(ic[ord,1], 1:length(coefs), ic[ord,2], 1:length(coefs), col = col_azul, lwd = 2)
abline(v = 0, col = col_rojo, lty = 2, lwd = 1.5)

7.5 Ejercicios

  1. Estime un modelo medv ~ rm + lstat con los datos Boston y realice la prueba t sobre cada coeficiente. ¿Ambas variables son significativas al 5%?
  2. Calcule el estadístico F manualmente para el modelo del ejercicio 1 y compárelo con el de summary().
  3. Construya el intervalo de confianza al 90% (no al 95%) para los coeficientes del modelo modelo_boston.
Ver soluciones 
# 1.
m1 <- lm(medv ~ rm + lstat, data = Boston)
summary(m1)$coefficients   # ambas variables tienen p-valor practicamente 0: significativas
##               Estimate Std. Error     t value     Pr(>|t|)
## (Intercept) -1.3582728 3.17282778  -0.4280953 6.687649e-01
## rm           5.0947880 0.44446550  11.4627299 3.472258e-27
## lstat       -0.6423583 0.04373146 -14.6886992 6.669365e-41
# 2.
SCT1 <- sum((Boston$medv - mean(Boston$medv))^2)
SCR1 <- sum(resid(m1)^2)
F1 <- ((SCT1-SCR1)/2) / (SCR1/(nrow(Boston)-2-1))
F1
## [1] 444.3309
summary(m1)$fstatistic
##    value    numdf    dendf 
## 444.3309   2.0000 503.0000
# 3.
confint(modelo_boston, level = 0.90)
##                     5 %         95 %
## (Intercept) -17.1299370 -6.463865958
## crim         -0.2027812 -0.079085483
## rm            7.0869256  8.375190310
## age          -0.1029433 -0.055905886
## dis          -1.2645216 -0.624785737
## tax          -0.0150857 -0.008020302

8 Diagnóstico de los supuestos clásicos

8.1 Multicolinealidad: el caso longley

El dataset longley (EE.UU., 1947-1962) es el ejemplo histórico por excelencia de multicolinealidad severa: varias variables macroeconómicas (PIB, deflactor, población) se mueven virtualmente juntas a lo largo del tiempo.

data(longley)
str(longley)
## 'data.frame':    16 obs. of  7 variables:
##  $ GNP.deflator: num  83 88.5 88.2 89.5 96.2 ...
##  $ GNP         : num  234 259 258 285 329 ...
##  $ Unemployed  : num  236 232 368 335 210 ...
##  $ Armed.Forces: num  159 146 162 165 310 ...
##  $ Population  : num  108 109 110 111 112 ...
##  $ Year        : int  1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 ...
##  $ Employed    : num  60.3 61.1 60.2 61.2 63.2 ...
modelo_longley <- lm(Employed ~ GNP + GNP.deflator + Population + Armed.Forces, data = longley)
summary(modelo_longley)$coefficients
##                   Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept)  120.323684241 17.808528631  6.756520 3.130014e-05
## GNP            0.096598412  0.017023424  5.674441 1.435173e-04
## GNP.deflator  -0.136325926  0.090948799 -1.498930 1.620313e-01
## Population    -0.658925145  0.178083033 -3.700101 3.500993e-03
## Armed.Forces  -0.004689174  0.002639816 -1.776326 1.033103e-01

El Factor de Inflación de Varianza (VIF) cuantifica cuánto se infla la varianza de \(\hat\beta_j\) por la correlación de \(X_j\) con las demás regresoras: \(VIF_j = 1/(1-R_j^2)\), donde \(R_j^2\) es el R² de regresar \(X_j\) sobre el resto de regresores.

car::vif(modelo_longley)
##          GNP GNP.deflator   Population Armed.Forces 
##   186.474439    62.742064    99.947942     2.198175

Valores de VIF muy por encima de 10 (como ocurre aquí con GNP y GNP.deflator) son evidencia contundente de multicolinealidad severa: estas variables comparten tanta información que el modelo no puede separar de forma confiable el efecto individual de cada una, aunque el \(R^2\) global del modelo siga siendo alto.

X_longley <- model.matrix(modelo_longley)[,-1]   # sin intercepto
ev_longley <- eigen(cor(X_longley))$values
kappa_longley <- max(ev_longley) / min(ev_longley)
kappa_longley    # numero de condicion: muy por encima de 30 -> multicolinealidad severa
## [1] 909.6962

8.2 Heterocedasticidad: la prueba de Breusch-Pagan

El supuesto de homocedasticidad (\(\operatorname{Var}(\varepsilon_i)=\sigma^2\) constante) puede fallar en datos de corte transversal. La prueba de Breusch-Pagan contrasta \(H_0:\) homocedasticidad.

lmtest::bptest(modelo_boston)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_boston
## BP = 33.014, df = 5, p-value = 3.739e-06
estilo_base()
plot(fitted(modelo_boston), resid(modelo_boston), pch = 19, col = col_azul, cex = 0.7,
     xlab = "Valores ajustados", ylab = "Residuos",
     main = "Residuos vs. ajustados - diagnostico de heterocedasticidad")
abline(h = 0, col = col_rojo, lwd = 2, lty = 2)

Un p-valor pequeño en el test de Breusch-Pagan (como el obtenido aquí) rechaza \(H_0\) de homocedasticidad: la varianza del error no es constante. En ese caso, es recomendable reportar errores estándar robustos (Sandwich/White) en vez de los estándar clásicos.

# Errores estandar clasicos vs. robustos a heterocedasticidad (HC1, tipo White)
ee_clasicos <- sqrt(diag(vcov(modelo_boston)))
ee_robustos <- sqrt(diag(sandwich::vcovHC(modelo_boston, type = "HC1")))
cbind(clasicos = round(ee_clasicos,3), robustos = round(ee_robustos,3))
##             clasicos robustos
## (Intercept)    3.236    4.496
## crim           0.038    0.032
## rm             0.391    0.664
## age            0.014    0.010
## dis            0.194    0.132
## tax            0.002    0.002

8.3 Autocorrelación: la prueba de Durbin-Watson

En datos de series de tiempo, es común que los errores estén correlacionados en el tiempo (\(\operatorname{Cov}(\varepsilon_t,\varepsilon_{t-1})\neq 0\)). La prueba de Durbin-Watson usa el estadístico \(DW\approx 2(1-\hat\rho)\), con valores cercanos a 2 indicando ausencia de autocorrelación de primer orden.

# longley es una serie temporal (1947-1962): pertinente probar autocorrelacion
lmtest::dwtest(modelo_longley)
## 
##  Durbin-Watson test
## 
## data:  modelo_longley
## DW = 1.4058, p-value = 0.02034
## alternative hypothesis: true autocorrelation is greater than 0

8.4 Normalidad de los residuos

Aunque MCO no requiere normalidad para ser insesgado, sí la requiere (en muestras pequeñas) para que las pruebas t y F sean exactas. Se evalúa con un gráfico Q-Q y la prueba de Shapiro-Wilk.

estilo_base()
qqnorm(resid(modelo_boston), pch = 19, col = col_azul, cex = 0.6,
       main = "Grafico Q-Q de los residuos - modelo Boston",
       xlab = "Cuantiles teoricos", ylab = "Cuantiles muestrales")
qqline(resid(modelo_boston), col = col_rojo, lwd = 2)

shapiro.test(resid(modelo_boston))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(modelo_boston)
## W = 0.83569, p-value < 2.2e-16

Un p-valor pequeño en Shapiro-Wilk sugiere que los residuos se apartan de la normalidad -algo frecuente en muestras grandes como esta (\(n=506\)), donde incluso desviaciones menores resultan estadísticamente detectables. En muestras grandes, el TCL (Sección 3.4) compensa parcialmente este problema para la inferencia sobre \(\hat\beta\).

8.5 Ejercicios

  1. Calcule el VIF de un modelo medv ~ rm + lstat + age + dis con los datos Boston. ¿Hay evidencia de multicolinealidad?
  2. Aplique el test de Breusch-Pagan al modelo modelo_longley. ¿Se rechaza la homocedasticidad?
  3. Construya el gráfico Q-Q de los residuos del modelo modelo_simple (Sección 4) con los datos cars.
Ver soluciones 
# 1.
m_vif <- lm(medv ~ rm + lstat + age + dis, data = Boston)
car::vif(m_vif)   # valores moderados (bien por debajo de 10): sin multicolinealidad grave
##       rm    lstat      age      dis 
## 1.675700 2.493339 2.762474 2.287533
# 2.
lmtest::bptest(modelo_longley)   # con pocas observaciones (n=16), el test tiene poca potencia
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_longley
## BP = 4.7746, df = 4, p-value = 0.3112
# 3.
estilo_base()
qqnorm(resid(modelo_simple), pch = 19, col = col_azul, main = "Q-Q residuos: modelo cars",
       xlab = "Cuantiles teoricos", ylab = "Cuantiles muestrales")
qqline(resid(modelo_simple), col = col_rojo, lwd = 2)

9 Variables cualitativas (dummies)

9.1 Codificación de variables categóricas

Muchas variables económicas son cualitativas (sexo, región, tipo de contrato). Se incorporan a la regresión mediante variables dummy, que toman valor 1 si se cumple una condición y 0 en caso contrario.

9.2 Aplicación con datos reales: mtcars

data(mtcars)
mtcars$transmision <- factor(mtcars$am, labels = c("Automatica","Manual"))
table(mtcars$transmision)
## 
## Automatica     Manual 
##         19         13
modelo_dummy <- lm(mpg ~ wt + hp + transmision, data = mtcars)
summary(modelo_dummy)$coefficients
##                      Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)       34.00287512 2.642659337 12.866916 2.824030e-13
## wt                -2.87857541 0.904970538 -3.180850 3.574031e-03
## hp                -0.03747873 0.009605422 -3.901830 5.464023e-04
## transmisionManual  2.08371013 1.376420152  1.513862 1.412682e-01

R codifica automáticamente transmisionManual como una dummy (1 si es manual, 0 si automática, con “Automática” como categoría base). El coeficiente indica el efecto diferencial promedio de tener transmisión manual sobre el consumo (mpg, millas por galón), controlando por peso (wt) y potencia (hp).

estilo_base()
plot(mtcars$wt, mtcars$mpg, pch = 19,
     col = ifelse(mtcars$transmision=="Manual", col_rojo, col_azul),
     xlab = "Peso (miles de libras)", ylab = "Consumo (millas/galon)",
     main = "mpg vs. peso, por tipo de transmision")
legend("topright", legend = c("Manual","Automatica"), col = c(col_rojo,col_azul),
       pch = 19, bty = "n", bg = "white", inset = 0.02, cex = 0.9)

Interpretación: manteniendo constante el peso y la potencia, los autos con transmisión manual rinden en promedio 2.08 millas por galón más que los automáticos.

9.3 Ejercicios

  1. Cree una dummy a partir de cyl en mtcars que valga 1 si el motor tiene 8 cilindros y 0 en caso contrario. Inclúyala en un modelo mpg ~ wt + dummy_8cil.
  2. Compare el \(R^2\) de ese modelo con y sin la dummy. ¿Mejora el ajuste?
Ver soluciones 
# 1.
mtcars$dummy_8cil <- ifelse(mtcars$cyl == 8, 1, 0)
m_dummy8 <- lm(mpg ~ wt + dummy_8cil, data = mtcars)
summary(m_dummy8)$coefficients
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 35.066777  2.1073123 16.640522 2.243338e-16
## wt          -4.252326  0.7638749 -5.566784 5.258466e-06
## dummy_8cil  -2.960814  1.4829268 -1.996602 5.533186e-02
# 2.
summary(lm(mpg ~ wt, data = mtcars))$r.squared          # sin dummy
## [1] 0.7528328
summary(m_dummy8)$r.squared                              # con dummy: mejora el ajuste
## [1] 0.782703

10 Ejercicios integradores finales

Los siguientes ejercicios integran álgebra lineal, estimación, inferencia y diagnóstico, usando el dataset real state.x77 (50 estados de EE.UU., datos del censo, circa 1970).

estados <- as.data.frame(state.x77)
colnames(estados) <- make.names(colnames(estados))
head(estados, 3)
##         Population Income Illiteracy Life.Exp Murder HS.Grad Frost   Area
## Alabama       3615   3624        2.1    69.05   15.1    41.3    20  50708
## Alaska         365   6315        1.5    69.31   11.3    66.7   152 566432
## Arizona       2212   4530        1.8    70.55    7.8    58.1    15 113417
  1. Estime, usando álgebra matricial pura (solve(), sin lm()), un modelo que explique Life.Exp (esperanza de vida) en función de Illiteracy (analfabetismo), Murder (tasa de homicidios) e Income (ingreso). Verifique con lm().
  2. Realice la prueba F de significancia conjunta del modelo anterior.
  3. Calcule el VIF de las tres variables explicativas. ¿Hay multicolinealidad preocupante?
  4. Aplique el test de Breusch-Pagan. ¿Se sostiene el supuesto de homocedasticidad?
  5. Construya el intervalo de confianza al 95% para el coeficiente de Murder e interprételo económicamente.
Ver soluciones 
y_e <- estados$Life.Exp
X_e <- as.matrix(cbind(Intercepto=1, estados[,c("Illiteracy","Murder","Income")]))

# 1.
beta_e <- solve(t(X_e)%*%X_e) %*% (t(X_e)%*%y_e)
beta_e
##                     [,1]
## Intercepto 71.1638497306
## Illiteracy  0.0368691885
## Murder     -0.2736317404
## Income      0.0003810966
modelo_estados <- lm(Life.Exp ~ Illiteracy + Murder + Income, data = estados)
coef(modelo_estados)
##   (Intercept)    Illiteracy        Murder        Income 
## 71.1638497306  0.0368691885 -0.2736317404  0.0003810966
# 2.
summary(modelo_estados)$fstatistic
##    value    numdf    dendf 
## 26.91581  3.00000 46.00000
# 3.
car::vif(modelo_estados)
## Illiteracy     Murder     Income 
##   2.348640   2.006166   1.254405
# 4.
lmtest::bptest(modelo_estados)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_estados
## BP = 3.3758, df = 3, p-value = 0.3372
# 5.
confint(modelo_estados)["Murder",]
##      2.5 %     97.5 % 
## -0.3657207 -0.1815428
# Interpretacion: por cada homicidio adicional (por cada 100,000 habitantes),
# la esperanza de vida cae, en promedio, dentro de ese rango, manteniendo
# constantes el analfabetismo y el ingreso.

11 Conclusiones

Este documento recorrió el camino completo de la econometría básica en R: partiendo del álgebra lineal que sostiene toda estimación (vectores, matrices, determinantes, sistemas de ecuaciones y autovalores), pasando por los fundamentos probabilísticos, hasta construir, estimar e interpretar modelos de regresión lineal simple y múltiple -verificando en cada paso que la fórmula matricial \(\hat\beta=(X'X)^{-1}X'y\) coincide exactamente con lo que entrega lm()-. Se estudiaron las propiedades del estimador MCO mediante simulación Monte Carlo, la inferencia estadística asociada (pruebas t, F e intervalos de confianza) y el diagnóstico completo de los supuestos clásicos (multicolinealidad, heterocedasticidad, autocorrelación y normalidad), cerrando con variables cualitativas y ejercicios integradores.

Todos los ejemplos utilizaron datos reales: cars (1930), swiss (1888), longley (1947-1962), Boston (1970), mtcars (1974) y state.x77 (1970) -eligiendo deliberadamente conjuntos de datos históricos y ampliamente documentados en la literatura estadística y econométrica, de modo que cada resultado numérico presentado en este documento sea reproducible y contrastable por cualquier lector.

12 Referencias

  • Wooldridge, J. M. (2020). Introductory Econometrics: A Modern Approach (7th ed.). Cengage Learning.
  • Greene, W. H. (2018). Econometric Analysis (8th ed.). Pearson.
  • Gujarati, D. N., & Porter, D. C. (2010). Econometría (5ª ed.). McGraw-Hill.
  • Harrison, D., & Rubinfeld, D. L. (1978). Hedonic housing prices and the demand for clean air. Journal of Environmental Economics and Management, 5(1), 81-102. (Dataset Boston)
  • R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
  • Fox, J., & Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). Sage. (Paquete car)

Documento elaborado por Jeel Cueva. Todo el código es reproducible: basta con abrir el archivo .Rmd en RStudio y ejecutar “Knit to PDF”.