Jak obliczane są standardowe błędy współczynników w regresji?

114

Dla własnego zrozumienia jestem zainteresowany ręczną replikacją obliczenia standardowych błędów szacowanych współczynników, ponieważ, na przykład, przychodzi z wyjściem lm()funkcji w R, ale nie byłem w stanie tego dokładnie określić. Jaka jest używana formuła / implementacja?

ako
źródło
8
dobre pytanie, wiele osób zna regresję z punktu widzenia algebry liniowej, gdzie rozwiązujesz równanie liniowe XXβ=Xy i otrzymujesz odpowiedź na beta. Nie jest jasne, dlaczego mamy za sobą standardowy błąd i założenie.
Haitao Du

Odpowiedzi:

122

Model liniowy jest zapisany jako gdzie y oznacza wektor odpowiedzi, β jest wektorem parametrów stałych efektów, X jest odpowiednią macierzą projektową, której kolumny są wartościami zmiennych objaśniających, a ϵ jest wektorem błędów losowych.

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

Jest dobrze wiadomo, że oszacowanie podaje się (patrz, na przykład, aby źródło Wikipedia ) β = ( X ' X ) - 1 X. ' y . Stąd Var ( β ) = ( X ' X ) - 1 X 'β

β^=(XX)1Xy.
[przypomnienie: Var ( A X ) = A × Var ( X ) × A , dla niektórych losowych wektorów X i niektórych nielosowych macierzy A ]
Var(β^)=(XX)1Xσ2IX(XX)1=σ2(XX)1,
Var(AX)=A×Var(X)×AXA

tak, że gdzie σ 2 można otrzymać przez błąd średni kwadratowy (MSE) w tabeli ANOVA.

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

Przykład z prostą regresją liniową w 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 
#----------------------

Gdy istnieje jedna zmienna objaśniająca, model zmniejsza się do

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
źródło
Dzięki za dokładną odpowiedź. Rozumiem, że ostatnia formuła nie dotyczy przypadku wielowymiarowego?
ako
1
Nie, ostatnia formuła działa tylko dla określonej macierzy X prostego modelu liniowego. W przypadku wielowymiarowym należy użyć ogólnej formuły podanej powyżej.
ocram
4
Var(β^)
6
Var(AX)=AVar(X)AXA
4
zauważ, że są to właściwe odpowiedzi do obliczeń ręcznych, ale rzeczywista implementacja zastosowana w lm.fit/ summary.lmjest nieco inna, dla stabilności i wydajności ...
Ben Bolker,
26

Formuły na nie można znaleźć w dowolnym tekście pośrednim dotyczącym statystyki, w szczególności można je znaleźć w Sheather (2009, rozdział 5) , z którego również pobierane jest następujące ćwiczenie (strona 138).

Poniższy kod R oblicza ręcznie oszacowania współczynników i ich standardowe błędy

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

co daje wynik

                         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

Porównaj z danymi wyjściowymi z lm():

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

co daje wynik:

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
źródło
Niezła sztuczka z solve()funkcją. Byłoby to nieco dłużej bez algebry macierzy. Czy istnieje zwięzły sposób wykonania tej konkretnej linii za pomocą podstawowych operatorów?
ako
1
β^=(XX)1XY(XX)
0

Część odpowiedzi Ocram jest błędna. Tak właściwie:

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

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

A komentarz do pierwszej odpowiedzi pokazuje, że potrzebne jest dalsze wyjaśnienie wariancji współczynnika:

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


Edytować

wronglywrong

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

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

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

Mam nadzieję, że to pomaga.

Linzhe Nie
źródło
1
β^=(XX)1XYβ^=(XX)1Xy(XX)1Xϵ
4
β^ϵ
Można to również obejrzeć w tym filmie: youtube.com/watch?v=jyBtfhQsf44
StatsStudent,