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 'β
β^=(X′X)−1X′y.
[przypomnienie:
Var ( A X ) = A × Var ( X ) × A ′ , dla niektórych losowych wektorów
X i niektórych nielosowych macierzy
A ]
Var(β^)=(X′X)−1X′σ2IX(X′X)−1=σ2(X′X)−1,
Var(AX)=A×Var(X)×A'XA
tak, że
gdzie σ 2 można otrzymać przez błąd średni kwadratowy (MSE) w tabeli ANOVA.
Varˆ(β^)=σ^2(X′X)−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=⎛⎝⎜⎜⎜⎜11⋮1x1x2⋮xn⎞⎠⎟⎟⎟⎟,β=(ab)
(X′X)−1=1n∑x2i−(∑xi)2(∑x2i−∑xi−∑xin)
Varˆ(b^)−−−−−−√=[σ^2(X′X)−1]22−−−−−−−−−−−√=nσ^2n∑x2i−(∑xi)2−−−−−−−−−−−−−−−√.
> num <- n * anova(mod)[[3]][2]
> denom <- n * sum(x^2) - sum(x)^2
> sqrt(num / denom)
[1] 0.09725302