¿Por qué R's lm () devuelve estimaciones de coeficientes diferentes a las de mi libro de texto?

13

Antecedentes

Estoy tratando de entender el primer ejemplo en un curso sobre modelos de adaptación (por lo que esto puede parecer ridículamente simple). He hecho los cálculos a mano y coinciden con el ejemplo, pero cuando los repito en R, los coeficientes del modelo están desactivados. Pensé que la diferencia puede deberse a que el libro de texto usa la varianza de la población ( ), mientras que R puede estar usando la varianza de la muestra ( S 2 ), pero no puedo ver dónde se usan en los cálculos. Por ejemplo, si se usa en alguna parte, la sección de ayuda sobre notas:σ2S2lm()var()var()

Se utiliza el denominador n - 1 que proporciona un estimador imparcial de la (co) varianza para las observaciones de iid.

He mirado el código para ambos lm()y lm.fit()como tampoco hacer uso devar() , pero lm.fit()pasa esos datos al código C compilado ( z <- .Call(C_Cdqrls, x, y, tol, FALSE)) al que no tengo acceso.

Pregunta

¿Alguien puede explicar por qué R está dando resultados diferentes? Incluso si hay una diferencia en el uso de la varianza muestra versus población, ¿por qué difieren las estimaciones de coeficientes?

Datos

Ajuste una línea para predecir el tamaño del zapato de grado en la escuela.

# 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

Editar

Como ha demostrado Ben Bolker , a veces parece que los maestros cometen errores. Parece que los cálculos de R son correctos. Moraleja de la historia: no creas algo solo porque un maestro dice que es verdad. ¡Compruébalo por ti mismo!

post-hoc
fuente
2
Retención doble mod.m=8/3. Porque si configuras mod.m=2.5714, entonces parecen ser idénticos.
Estadísticas
2
Los coeficientes mod.m = 8/3 y mod.b = -1 no se calculan en ninguna parte de los comentarios, por lo que entiendo, por lo que no es obvio. Como comenta @Stat arriba, el error parece estar en el mod de computación.
Juho Kokkala
2
Es importante tener en cuenta que cualquiera puede cometer errores: su maestro, usted, los que responden aquí, los programadores de R, cualquiera. Entonces, cuando trate de averiguar dónde pueden estar los errores cuando las cosas no están de acuerdo, considere cuántas otras personas están revisando cada cosa. En el caso de la lmfunción en R, literalmente decenas de miles de personas han verificado los resultados comparándolos con otras cosas, y la salida de lmse compara con ejemplos conocidos cada vez que algo cambia en el código. Con las respuestas aquí, es probable que al menos algunas personas verifiquen (su pregunta se ha examinado 29 veces).
Glen_b -Reinstalar Monica
1
@Glen_b Tu punto es en realidad la razón por la que vine a preguntar. No podía entender cómo R podía estar equivocado en un cálculo tan básico, pero no podía entender por qué eran diferentes. Evento husmeó alrededor del código fuente. Pero al final, el error fue en el último lugar en el que pensé mirar, principalmente porque la parte de cálculo está en los límites de mi conocimiento. Sin embargo, ¡aprendí mucho de la respuesta!
post-hoc
2
Sí, es importante tratar de descubrir por qué difieren; tiene sentido preguntar aquí si no puedes resolverlo. Estaba tratando de sugerir por qué el último lugar que consideró podría haber sido uno de los primeros lugares para buscar. He sido atrapado haciendo cambios 'simplificadores' de último minuto en ejemplos en una o dos ocasiones.
Glen_b: reinstala a Monica el

Respuestas:

25

Parece que el autor cometió un error matemático en alguna parte.

Si expande la desviación de la suma de cuadrados

S=((si+metro)-1)2+((si+2metro)-5 5)2+((si+4 4metro)-9 9)2
S=si2+2simetro+metro2+1-2si-2metro+si2+4 4simetro+4 4metro2+25-10si-20metro+si2+8simetro+dieciséismetro2+81-18 añossi-72metro

3si2+14simetro+21metro2+107-30si-94metro

Ssimetro

reS/ /resi=6 6si+14metro-303si+7 7metro-15=0 0
reS/ /remetro=14si+42metro-947 7si+21metro-47=0 0

Resolver

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

R says this is indeed 2.571429 ...

Based on this link this seems to be from a Coursera course ... ? Maybe there was a mis-transcription of the data somewhere?

The other, independent way to do this calculation is to know that the estimated regression slope is equal to the sum of cross products ((yy¯)(xx¯)) divided by the sum of squares ((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

If think if the shoe sizes were {1,11/3,9} instead of {1,5,9} then the slope would come out to 8/3 ...

Ben Bolker
fuente
2
Wow. Yes, you are right. It's from a Coursera course and it's from the video, not transcription. So I'm guessing he simplified it to make the calculations simpler for the video and didn't expect anyone to try and repeat it. It just happened to be the first video that I saw so I tried to follow along. It's clear that I need to upskill when it comes to maths. I think found the error though. The constant term, which you say doesn't matter, is probably the correct value which through off his calculations. I'll look through your answer a few more times to teach myself. I really appreciate it!
post-hoc
I don't think the constant term will throw off the calculations. It won't affect the estimates of the slope and intercept (it disappears when we take the derivative), only the estimates of the residual SSQ/standard deviation.
Ben Bolker