Pourquoi lm () de R renvoie-t-il des estimations de coefficient différentes de celles de mon manuel?

13

Contexte

J'essaie de comprendre le premier exemple d'un cours sur l'ajustement de modèles (cela peut donc sembler ridiculement simple). J'ai fait les calculs à la main et ils correspondent à l'exemple, mais quand je les répète en R, les coefficients du modèle sont désactivés. Je pensais que la différence peut être due au manuel utilisant la variance de la population ( σ2 ) tandis que R peut utiliser la variance de l'échantillon ( S2 ), mais je ne vois pas où ceux-ci sont utilisés dans les calculs. Par exemple, si vous l' lm()utilisez var()quelque part, la section d'aide sur les var()notes:

Le dénominateur n - 1 est utilisé, ce qui donne un estimateur non biaisé de la (co) variance pour les observations iid.

J'ai regardé le code pour les deux lm()et je ne m'en lm.fit()sers pas var(), mais lm.fit()passe ces données au code C compilé ( z <- .Call(C_Cdqrls, x, y, tol, FALSE)) auquel je n'ai pas accès.

Question

Quelqu'un peut-il expliquer pourquoi R donne des résultats différents? Même s'il existe une différence dans l'utilisation de la variance échantillon vs population, pourquoi les estimations des coefficients diffèrent-elles?

Les données

Ajustez une ligne pour prédire la taille des chaussures à partir du niveau scolaire.

# model data
mod.dat <- read.table(
    text = 'grade shoe
                1    1
                2    5
                4    9'
    , header = T);

# mean
mod.mu  <- mean(mod.dat$shoe);
# variability 
mod.var <- sum((mod.dat$shoe - mod.mu)^2)

# model coefficients from textbook
mod.m  <- 8/3;
mod.b  <- -1;

# predicted values  ( 1.666667 4.333333 9.666667 )
mod.man.pred       <- mod.dat$grade * mod.m + mod.b;
# residuals         ( -0.6666667  0.6666667 -0.6666667 )
mod.man.resid      <- (mod.dat$shoe - mod.man.pred)
# residual variance ( 1.333333 )
mod.man.unexpl.var <- sum(mod.man.resid^2);
# r^2               ( 0.9583333 )
mod.man.expl.var   <- 1 - mod.man.unexpl.var / mod.var;

# but lm() gives different results:
summary(lm(shoe ~ grade, data = mod.dat))
Call:
lm(formula = shoe ~ grade, data = mod.dat)

Residuals:
      1       2       3 
-0.5714  0.8571 -0.2857 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -1.0000     1.3093  -0.764    0.585
grade         2.5714     0.4949   5.196    0.121

Residual standard error: 1.069 on 1 degrees of freedom
Multiple R-squared:  0.9643,    Adjusted R-squared:  0.9286 
F-statistic:    27 on 1 and 1 DF,  p-value: 0.121

Éditer

Comme l' a montré Ben Bolker , il semble que les enseignants font parfois des erreurs. Il semble que les calculs de R soient corrects. Morale de l'histoire: ne croyez pas quelque chose simplement parce qu'un enseignant dit que c'est vrai. Vérifiez par vous-même!

post-hoc
la source
2
Revérifiez mod.m=8/3. Parce que si vous définissez mod.m=2.5714, ils semblent être identiques.
Stat
2
Les coefficients mod.m = 8/3 et mod.b = -1 ne sont calculés nulle part dans les commentaires pour autant que je sache, donc ce n'est pas évident. Comme le commente @Stat ci-dessus, l'erreur semble être dans le calcul du mod.m.
Juho Kokkala
2
Il est important de garder à l'esprit que n'importe qui peut faire des erreurs - votre professeur, vous, les répondeurs ici, les programmeurs R - n'importe qui. Donc, lorsque vous essayez de comprendre où les erreurs peuvent se trouver lorsque les choses ne sont pas d'accord, considérez combien d'autres personnes vérifient chaque chose. Dans le cas de la lmfonction dans R, des dizaines de milliers de personnes ont littéralement vérifié les résultats en les comparant à d'autres choses, et la sortie de lmest vérifiée par rapport à des exemples connus chaque fois que quelque chose change dans le code. Avec les réponses ici, au moins quelques personnes sont susceptibles de vérifier (votre question a été consultée 29 fois).
Glen_b -Reinstate Monica
1
@Glen_b Votre point est en fait la raison pour laquelle je suis venu ici pour demander. Je ne pouvais pas comprendre comment R pouvait se tromper sur un calcul aussi basique, mais je ne pouvais pas comprendre pourquoi ils étaient différents. J'ai fouillé le code source. Mais au final, l'erreur a été au dernier endroit que je pensais chercher, principalement parce que la partie calcul est à la limite de mes connaissances. J'ai cependant beaucoup appris de la réponse!
post-hoc
2
Oui, il est important d'essayer de comprendre pourquoi ils diffèrent; il est logique de demander ici si vous ne pouvez pas le résoudre. J'essayais de suggérer pourquoi le dernier endroit que vous avez envisagé aurait plutôt pu être l'un des premiers endroits à chercher. J'ai moi-même été surpris par des modifications de dernière minute «simplifiantes» aux exemples à une ou deux reprises.
Glen_b -Reinstate Monica

Réponses:

25

Il semble que l'auteur ait fait une erreur mathématique quelque part.

Si vous augmentez l'écart de somme des carrés

S=((b+m)-1)2+((b+2m)-5)2+((b+4m)-9)2
S=b2+2bm+m2+12b2m+b2+4bm+4m2+2510b20m+b2+8bm+16m2+8118b72m

3b2+14bm+21m2+10730b94m

Sbm

dS/db=6b+14m303b+7m15=0
dS/dm=14b+42m947b+21m47=0

Résoudre

b=(157m)/30=7(157m)/3+21m474735=(49/3+21)mm=(4735)/(2149/3)=18/7

R dit que c'est bien 2,571429 ...

Sur la base de ce lien, cela semble provenir d'un cours Coursera ...? Peut-être qu'il y a eu une mauvaise transcription des données quelque part?

(yy¯)(xx¯)(xx¯)2

g <- c(1,2,4)
g0 <- g - mean(g)
s <- c(1,5,9)
s0 <- s- mean(s)
sum(g0*s0)/(sum(g0^2))
## [1] 2.571429

{1,11/3,9}{1,5,9}

Ben Bolker
la source
2
Sensationnel. Oui, tu as raison. Cela vient d'un cours de Coursera et c'est de la vidéo, pas de transcription. Je suppose donc qu'il l'a simplifié pour rendre les calculs plus simples pour la vidéo et ne s'attendait pas à ce que quelqu'un essaie de le répéter. Il se trouve que c'est la première vidéo que j'ai vue, j'ai donc essayé de suivre. Il est clair que je dois me perfectionner en matière de mathématiques. Je pense que j'ai trouvé l'erreur. Le terme constant, dont vous dites qu'il n'a pas d'importance, est probablement la valeur correcte qui, à travers ses calculs. Je vais parcourir votre réponse encore quelques fois pour m'enseigner. J'apprécie vraiment cela!
post-hoc le
Je ne pense pas que le terme constant détruira les calculs. Cela n'affectera pas les estimations de la pente et de l'ordonnée à l'origine (il disparaît lorsque nous prenons la dérivée), seulement les estimations de la SSQ résiduelle / écart-type.
Ben Bolker