1 Datos

Supongamos que conocemos:

# y = a + bx + u, u ~ N(0,1)
a = 2
b = 0.5
sigma = 1
x = 1:20

Podemos simular la variable y:

set.seed(2) # para obtener siempre los mismos resultados
u = rnorm(20, 0, sigma)
# datos
y = a + b*x + u
# dibujamos los datos obtenidos (x,y)
plot(x,y, xlim = c(0,20), ylim = c(0,15), pch = 19)
# añadimos la recta "teórica" entre x e y
abline(a, b, lty = 2, lwd = 1, col = "red")

2 Precisión de los estimadores mediante probabilidad

Ya hemos visto que

\[ \hat \beta \sim N(\beta, (X^TX)^{-1} \sigma^2) \]

Este resultado se deriva matemáticamente a partir de la distribución

\[ \epsilon \sim N(0, \sigma^2 I) \]

Estimamos a y b utilizando x e y:

m = lm(y ~ x)
(beta_e = coef(m))

## (Intercept)           x 
##   2.0610497   0.5128011

Podemos dibujar la recta estimada y la teórica:

plot(x,y, xlim = c(0,20), ylim = c(0,15), pch = 19)
abline(a, b, lty = 2, lwd = 1, col = "red")
abline(beta_e[1], beta_e[2], lwd = 1, col = "blue")
legend(0, 15, legend=c("teorica", "estimada"), col=c("red", "blue"), lty=c(2,1))

La varianza de los estimadores (o el standard error) es una manera de calcular su precisión:

(beta_var = diag(vcov(m)))

## (Intercept)           x 
## 0.253927377 0.001769529

Otra manera es mediante los intervalos de confianza:

confint(m)

##                 2.5 %   97.5 %
## (Intercept) 1.0023697 3.119730
## x           0.4244242 0.601178

3 Precisión de los estimadores mediante simulaciones

Otra manera de obtener los resultados del apartado anterior es utilizando simulaciones en lugar de matemáticas. En este caso lo podemos hacer porque conocemos a, b y \(\sigma\):

beta_e_sim = matrix(0, nrow = 500, ncol = 2)
for (i in 1:500){
  ui = rnorm(20, 0, sigma)
  yi = a + b*x + ui
  mi = lm(yi ~ x)
  beta_e_sim[i,] = coef(mi)
}

Al final hemos obtenido 500 valores de \(\hat a\) y $ b$. Podemos dibujar el histograma:

par(mfrow = c(1,2))
hist(beta_e_sim[,1], freq = F, xlab = "a estimada", main = "")
hist(beta_e_sim[,1], freq = F, xlab = "b estimada", main = "")

EL valor final que elegimos para \(\hat a\) y \(\hat b\) es la media de los datos:

apply(beta_e_sim, 2, mean)

## [1] 2.0211790 0.4990019

La varianza de los estimadores es la varianza de los datos correspondientes:

apply(beta_e_sim, 2, var)

## [1] 0.21604665 0.00151284

El intervalo de confianza se puede calcular a partir de los cuantiles de la distribución:

apply(beta_e_sim, 2, quantile, probs = c(0.025,0.975))

##           [,1]      [,2]
## 2.5%  1.053055 0.4207785
## 97.5% 2.859275 0.5744381

4 Precisión de los estimadores mediante bootstrap

El problema con el método de simulaciones es que se tiene que conocer {a, b, \(\sigma\)}, y eso nunca ocurre. En la práctica solo se conoce {\(x_i, y_i\)}. Una alternativa es utilizar el método bootstrap:

# residuos
e = resid(m)
# muestreamos los residuos CON REPOSICION
beta_e_b = matrix(0, nrow = 500, ncol = 2)
for (i in 1:500){
  eb = sample(e, rep = T)
  yb = beta_e[1] + beta_e[2]*x + eb
  mb = lm(yb ~ x)
  beta_e_b[i,] = coef(mb)
}

Podemos calcular la varianza de los estimadores como:

apply(beta_e_b, 2, var)

## [1] 0.23774930 0.00159955

Y los intervalos de confianza:

apply(beta_e_b, 2, quantile, probs = c(0.025,0.975))

##           [,1]      [,2]
## 2.5%  1.048589 0.4388716
## 97.5% 2.929364 0.5975489

El método bootstrap cobra importancia cuando la distribución en el muestreo teórica no se puede calcular matemáticamente.

5 Ejemplo

Vamos a calcular la precisión de los parámetros del modelo siguiente utilizando probabilidad y bootstrap:

d = faraway::gala
m1 = lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, data = d)
beta_e = coef(m1)

# BOOTSTRAP
# residuos
e = resid(m1)
# matrix de las X
X = model.matrix(m1)
# muestreamos los residuos CON REPOSICION
beta_e_b = matrix(0, nrow = 500, ncol = 6)
for (i in 1:500){
  eb = sample(e, rep = T)
  yb = X %*% beta_e + eb
  mb = lm(yb ~ Area + Elevation + Nearest + Scruz + Adjacent, data = d)
  beta_e_b[i,] = coef(mb)
}

• Varianza de los parémetros estimados

# con probabilidad
diag(vcov(m1))

##  (Intercept)         Area    Elevation      Nearest        Scruz 
## 3.668833e+02 5.027618e-04 2.879697e-03 1.111203e+00 4.639813e-02 
##     Adjacent 
## 3.132967e-04

# con bootstrap
apply(beta_e_b, 2, var)

## [1] 2.857035e+02 4.113618e-04 2.472941e-03 9.982770e-01 4.040625e-02
## [6] 2.799286e-04

• intervalos de confianza

# con probabilidad
confint(m1)

##                   2.5 %      97.5 %
## (Intercept) -32.4641006 46.60054205
## Area         -0.0702158  0.02233912
## Elevation     0.2087102  0.43021935
## Nearest      -2.1664857  2.18477363
## Scruz        -0.6850926  0.20404416
## Adjacent     -0.1113362 -0.03827344

# con bootstrap
t(apply(beta_e_b, 2, quantile, probs = c(0.025,0.975)))

##              2.5%       97.5%
## [1,] -24.42202493 38.76288849
## [2,]  -0.06102367  0.02540897
## [3,]   0.22562407  0.42286039
## [4,]  -1.70054626  2.31979901
## [5,]  -0.61829317  0.15621835
## [6,]  -0.10266255 -0.03516742