Comment les erreurs-types des coefficients sont-elles calculées dans une régression?

115

Pour ma propre compréhension, je suis intéressé par la réplication manuelle du calcul des erreurs types des coefficients estimés, car, par exemple, le résultat de la lm()fonction est fourni R, mais je n’ai pas pu le localiser. Quelle est la formule / mise en œuvre utilisée?

ako
la source
8
bonne question, beaucoup de gens connaissent la régression du point d'algèbre linéaire de vue, où vous résoudre l'équation linéaire XXβ=Xy et obtenir la réponse pour la bêta. Pas clair pourquoi nous avons l'erreur type et l'hypothèse derrière elle.
Haitao Du

Réponses:

123

Le modèle linéaire est écrit comme Y désigne le vecteur des réponses, est le vecteur des paramètres d'effets fixes, est la matrice de plan correspondant dont les colonnes sont les valeurs des variables explicatives et est le vecteur des erreurs aléatoires.

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

Il est bien connu qu'une estimation de est donnée par (voir, par exemple, l'article de Wikipédia ) D'où [rappel: , pour un vecteur aléatoire et pour une matrice non aléatoire ]ß = ( X ' X ) - 1 X ' y . Var ( β ) = ( X ' X ) - 1 X 'β

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

afin que où peut être obtenu par l'erreur quadratique moyenne (MSE) dans la table ANOVA. σ 2

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

Exemple avec une régression linéaire simple 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 
#----------------------

Lorsqu'il n'y a qu'une seule variable explicative, le modèle se réduit à et pour que et les formules deviennent plus transparentes. Par exemple, l’erreur type de la pente estimée est X = ( 1 x 1 1 x 21 x n ) ,

yi=a+bxi+ϵi,i=1,,n
( X ' X ) - 1 = 1
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
la source
Merci pour la réponse complète. Donc, je suppose que la dernière formule ne tient pas dans le cas multivarié?
ako
1
Non, la toute dernière formule ne fonctionne que pour la matrice X spécifique du modèle linéaire simple. Dans le cas multivarié, vous devez utiliser la formule générale donnée ci-dessus.
ocram
4
+1, une petite question, comment vient ? Var(β^)
avocat le
6
@loganecolss: Cela vient du fait que , pour un vecteur aléatoire et une matrice non aléatoire . X AVar(AX)=AVar(X)AXA
ocram
4
Notez que ce sont les bonnes réponses pour le calcul manuel, mais l'implémentation réelle utilisée dans lm.fit/ summary.lmest un peu différente, pour la stabilité et l'efficacité ...
Ben Bolker
26

Les formules pour celles-ci peuvent être trouvées dans n'importe quel texte intermédiaire sur les statistiques, en particulier dans Sheather (2009, chapitre 5) , à partir duquel l'exercice suivant est également réalisé (page 138).

Le code R suivant calcule les estimations de coefficients et leurs erreurs types manuellement

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

qui produit la sortie

                         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

Comparez à la sortie de lm():

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

qui produit la sortie:

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
la source
Belle astuce avec la solve()fonction. Ce serait un peu plus long sans l'algèbre matricielle. Existe-t-il un moyen succinct de réaliser cette ligne spécifique avec uniquement des opérateurs de base?
ako
1
@AkselO Il existe l'expression de formulaire fermé bien connue de l'estimateur OLS, , que vous pouvez calculer en calculant explicitement l'inverse de la matrice ( (comme l'a fait @ ocram), mais cela devient délicat avec des matrices mal conditionnées. β^=(XX)1XY(XX)
tchakravarty
0

Une partie de la réponse d'Ocram est fausse. Réellement:

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

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

Et le commentaire de la première réponse montre qu'il est nécessaire d'expliquer davantage la variance du coefficient:

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


Modifier

Merci, ignoré le chapeau de cette version bêta. La déduction ci-dessus est . Le résultat correct est:wronglywrong

1.(Pour obtenir cette équation, définissez la dérivée de premier ordre de sur sur zéro, pour maximiser )β^=(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

J'espère que ça aide.

Linzhe Nie
la source
1
La dérivation de l'estimateur OLS pour le vecteur bêta, , se trouve dans n'importe quel manuel de régression décent. À la lumière de cela, pouvez-vous fournir la preuve que ce doit être place? β =(X'X)-1X'Y-(X'X)-1X'εβ^=(XX)1XYβ^=(XX)1Xy(XX)1Xϵ
gung
4
Votre n'est même pas un estimateur, car n'est pas observable! eβ^ϵ
whuber
Ceci peut également être visionné dans cette vidéo: youtube.com/watch?v=jyBtfhQsf44
StatsStudent