Para mi propio entendimiento, estoy interesado en el manual de replicar el cálculo de los errores estándar de los coeficientes estimados como, por ejemplo, vienen con la salida de la lm() de la función en R, pero no he sido capaz de precisar. ¿Cuál es la fórmula / aplicación utiliza?
Respuestas
¿Demasiados anuncios?
El modelo lineal se escribe como $$ \left| \begin{array}{l} \mathbf{y} = \mathbf{X} \mathbf{\beta} + \mathbf{\epsilon} \\ \mathbf{\epsilon} \sim N(0, \sigma^2 \mathbf{I}), \end{array} \right.$$ donde $\mathbf{y}$ denota el vector de respuestas, $\mathbf{\beta}$ es el vector de efectos fijos de los parámetros, $\mathbf{X}$ es el correspondiente diseño de la matriz cuyas columnas son los valores de las variables explicativas, y $\mathbf{\epsilon}$ es el vector de errores aleatorios.
Es bien sabido que un estimado de $\mathbf{\beta}$ está dada por (consulte, por ejemplo, para el artículo de la wikipedia) $$\hat{\mathbf{\beta}} = (\mathbf{X}^{\prime} \mathbf{X})^{-1} \mathbf{X}^{\prime} \mathbf{y}.$$ Por lo tanto $$ \textrm{Var}(\hat{\mathbf{\beta}}) = (\mathbf{X}^{\prime} \mathbf{X})^{-1} \mathbf{X}^{\prime} \;\sigma^2 \mathbf{I} \; \mathbf{X} (\mathbf{X}^{\prime} \mathbf{X})^{-1} = \sigma^2 (\mathbf{X}^{\prime} \mathbf{X})^{-1}, $$ así que $$ \widehat{\textrm{Var}}(\hat{\mathbf{\beta}}) = \hat{\sigma}^2 (\mathbf{X}^{\prime} \mathbf{X})^{-1}, $$ donde $\hat{\sigma}^2$ se puede obtener mediante el Error cuadrático medio (MSE) en la tabla ANOVA.
Ejemplo con una simple regresión lineal en R
#------generate one data set with epsilon ~ N(0, 0.25)------
seed <- 1152 #seed
n <- 100 #nb of observations
a <- 5 #intercept
b <- 2.7 #slope
set.seed(seed)
epsilon <- rnorm(n, mean=0, sd=sqrt(0.25))
x <- sample(x=c(0, 1), size=n, replace=TRUE)
y <- a + b * x + epsilon
#-----------------------------------------------------------
#------using lm------
mod <- lm(y ~ x)
#--------------------
#------using the explicit formulas------
X <- cbind(1, x)
betaHat <- solve(t(X) %*% X) %*% t(X) %*% y
var_betaHat <- anova(mod)[[3]][2] * solve(t(X) %*% X)
#---------------------------------------
#------comparison------
#estimate
> mod$coef
(Intercept) x
5.020261 2.755577
> c(betaHat[1], betaHat[2])
[1] 5.020261 2.755577
#standard error
> summary(mod)$coefficients[, 2]
(Intercept) x
0.06596021 0.09725302
> sqrt(diag(var_betaHat))
x
0.06596021 0.09725302
#----------------------
Cuando hay una sola variable explicativa, el modelo se reduce a $$y_i = a + bx_i + \epsilon_i, \qquad i = 1, \dotsc, n$$ y $$\mathbf{X} = \left( \begin{array}{cc} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_n \end{array} \right) \qquad \mathbf{\beta} = \left( \begin{array}{c} a\\b \end{array} \right)$$ así que $$(\mathbf{X}^{\prime} \mathbf{X})^{-1} = \frac{1}{n\sum x_i^2 - (\sum x_i)^2} \left( \begin{array}{cc} \sum x_i^2 & -\sum x_i \\ -\sum x_i & n \end{array} \right)$$ y las fórmulas se vuelven más transparentes. Por ejemplo, el error estándar de la estimación de la pendiente es $$\sqrt{\widehat{\textrm{Var}}(\hat{b})} = \sqrt{[\hat{\sigma}^2 (\mathbf{X}^{\prime} \mathbf{X})^{-1}]_{22}} = \sqrt{\frac{n \hat{\sigma}^2}{n\sum x_i^2 - (\sum x_i)^2}}.$$
> num <- n * anova(mod)[[3]][2]
> denom <- n * sum(x^2) - sum(x)^2
> sqrt(num / denom)
[1] 0.09725302
Jamie Brennan
Puntos
86
Las fórmulas para estos se pueden encontrar en cualquier intermedio de texto en las estadísticas, en particular, usted puede encontrarlos en Sheather (2009, Capítulo 5), desde donde el siguiente ejercicio es también de la toma (página 138).
El siguiente código R calcula el coeficiente de estimaciones y sus errores estándar manualmente
dfData <- as.data.frame(
read.csv("http://www.stat.tamu.edu/~sheather/book/docs/datasets/MichelinNY.csv",
header=T))
# using direct calculations
vY <- as.matrix(dfData[, -2])[, 5] # dependent variable
mX <- cbind(constant = 1, as.matrix(dfData[, -2])[, -5]) # design matrix
vBeta <- solve(t(mX)%*%mX, t(mX)%*%vY) # coefficient estimates
dSigmaSq <- sum((vY - mX%*%vBeta)^2)/(nrow(mX)-ncol(mX)) # estimate of sigma-squared
mVarCovar <- dSigmaSq*chol2inv(chol(t(mX)%*%mX)) # variance covariance matrix
vStdErr <- sqrt(diag(mVarCovar)) # coeff. est. standard errors
print(cbind(vBeta, vStdErr)) # output
el cual se produce la salida
vStdErr
constant -57.6003854 9.2336793
InMichelin 1.9931416 2.6357441
Food 0.2006282 0.6682711
Decor 2.2048571 0.3929987
Service 3.0597698 0.5705031
Comparar a la salida de lm():
# using lm()
names(dfData)
summary(lm(Price ~ InMichelin + Food + Decor + Service, data = dfData))
el cual se produce la salida:
Call:
lm(formula = Price ~ InMichelin + Food + Decor + Service, data = dfData)
Residuals:
Min 1Q Median 3Q Max
-20.898 -5.835 -0.755 3.457 105.785
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -57.6004 9.2337 -6.238 3.84e-09 ***
InMichelin 1.9931 2.6357 0.756 0.451
Food 0.2006 0.6683 0.300 0.764
Decor 2.2049 0.3930 5.610 8.76e-08 ***
Service 3.0598 0.5705 5.363 2.84e-07 ***
---
Signif. codes: 0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1
Residual standard error: 13.55 on 159 degrees of freedom
Multiple R-squared: 0.6344, Adjusted R-squared: 0.6252
F-statistic: 68.98 on 4 and 159 DF, p-value: < 2.2e-16
