Как стандартные ошибки коэффициентов рассчитываются в регрессии?

115

Для моего собственного понимания я заинтересован в том, чтобы вручную повторить вычисление стандартных ошибок оценочных коэффициентов, поскольку, например, они поставляются с выходными данными lm()функции R, но не смогли ее определить. Какая формула / реализация используется?

АКО
источник
8
Хороший вопрос, многие знают регрессию с точки зрения линейной алгебры, где вы решаете линейное уравнение и получаете ответ для бета-версии. Непонятно, почему у нас есть стандартная ошибка и предположение за этим. XXβ=Xy
Haitao Du

Ответы:

123

Линейная модель записывается как где обозначает вектор ответов, - вектор параметров фиксированных эффектов, - соответствующая матрица дизайна, столбцы которой являются значениями объясняющих переменных, и - вектор случайных ошибок. у β Х & epsi ;

|y=Xβ+ϵϵN(0,σ2I),
yβXϵ

Хорошо известно, что оценка дается (см., Например, статью в Википедии ) Следовательно, [напоминание: для некоторого случайного вектора и некоторой неслучайной матрицы ]& beta ; = ( Х ' х ) - 1 х ' у . Вар ( β ) = ( X ' X ) - 1 Х 'β

β^=(XX)1Xy.
вар ( Х ) = A × Var ( X ) × '
Var(β^)=(XX)1Xσ2IX(XX)1=σ2(XX)1,
Var(AX)=A×Var(X)×AXA

так что где можно получить с помощью ошибки среднего квадрата (MSE) в таблице ANOVA.

Var^(β^)=σ^2(XX)1,
σ^2

Пример с простой линейной регрессией в 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 
#----------------------

Когда есть одна объясняющая переменная, модель сводится к и чтобы и формулы становятся более прозрачными. Например, стандартная ошибка предполагаемого уклона:

yi=a+bxi+ϵi,i=1,,n
X=(1x11x21xn),β=(ab)
(XX)1=1nxi2(xi)2(xi2xixin)
Var^(b^)=[σ^2(XX)1]22=nσ^2nxi2(xi)2.
> num <- n * anova(mod)[[3]][2]
> denom <- n * sum(x^2) - sum(x)^2
> sqrt(num / denom)
[1] 0.09725302
ocram
источник
Спасибо за подробный ответ. Итак, я так понимаю, последняя формула не выполняется в многомерном случае?
АКО
1
Нет, самая последняя формула работает только для конкретной X-матрицы простой линейной модели. В многомерном случае вы должны использовать общую формулу, приведенную выше.
Октябрь
4
+1, быстрый вопрос, как получается ? Var(β^)
Авокадо
6
@loganecolss: Он исходит из того , что для некоторого случайного вектора и некоторой неслучайной матрицы . X AVar(AX)=AVar(X)AXA
октября
4
обратите внимание, что это правильные ответы для ручного расчета, но фактическая реализация, используемая в lm.fit/ summary.lmнемного отличается, для стабильности и эффективности ...
Бен Болкер,
26

Формулы для них можно найти в любом промежуточном тексте по статистике, в частности, вы можете найти их в Sheather (2009, глава 5) , откуда также взято следующее упражнение (стр. 138).

Следующий код R вычисляет оценки коэффициентов и их стандартные ошибки вручную

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

который производит вывод

                         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

Сравните с выводом из lm():

# using lm()
names(dfData)
summary(lm(Price ~ InMichelin + Food + Decor + Service, data = dfData))

который производит вывод:

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 
tchakravarty
источник
Хороший трюк с solve()функцией. Это было бы немного дольше без матричной алгебры. Есть ли краткий способ выполнения этой конкретной строки только с базовыми операторами?
АКО
1
@AkselO Для оценки OLS существует хорошо известное выражение в закрытой форме, , который можно вычислить, явно вычислив обратную матрицу ( (как это сделал @ ocram), но это сложно сделать с плохо обусловленными матрицами. β^=(XX)1XY(XX)
Чакраварти
0

Часть ответа Окрама неверна. Фактически:

β^=(XX)1Xy(XX)1Xϵ.

E(β^)=(XX)1Xy.

А комментарий к первому ответу показывает, что требуется более подробное объяснение дисперсии коэффициента:

Var(β^)=E(β^E(β^))2=Var((XX)1Xϵ)=(XX)1Xσ2IX(XX)1=σ2(XX)1


редактировать

Спасибо, я проигнорировал шляпу в этой бета-версии. Вышеприведенный вывод . Правильный результат:wronglywrong

1.(Чтобы получить это уравнение, установите производную первого порядка от в равную нулю, для максимизации )β^=(XX)1Xy.SSRβSSR

2.E(β^|X)=E((XX)1X(Xβ+ϵ)|X)=β+((XX)1X)E(ϵ|X)=β.

3.Var(β^)=E(β^E(β^|X))2=Var((XX)1Xϵ)=(XX)1Xσ2IX(XX)1=σ2(XX)1

Надеюсь, это поможет.

Линже Не
источник
1
Вывод оценки OLS для бета-вектора, , можно найти в любом учебнике по достойной регрессии. В свете этого, можете ли вы предоставить доказательство того, что это должно быть вместо? β =(Х'х)-1х'у-(Х'х)-1Х'εβ^=(XX)1XYβ^=(XX)1Xy(XX)1Xϵ
Кун
4
Ваша даже не оценщик, потому что не наблюдается! ; & epsiβ^ϵ
whuber
Это также можно посмотреть в этом видео: youtube.com/watch?v=jyBtfhQsf44
StatsStudent