Por que lm () de R retorna estimativas de coeficientes diferentes das do meu livro?

13

fundo

Estou tentando entender o primeiro exemplo de um curso sobre montagem de modelos (então isso pode parecer ridiculamente simples). Fiz os cálculos manualmente e eles correspondem ao exemplo, mas quando os repito em R, os coeficientes do modelo estão desativados. Eu pensei que a diferença pode ser devido ao livro didático usar variação populacional ( ), enquanto R pode estar usando variação amostral ( ), mas não consigo ver onde elas são usadas nos cálculos. Por exemplo, se usar em algum lugar, a seção de ajuda nas notas:σ2S2lm()var()var()

Utiliza-se o denominador n - 1, que fornece um estimador imparcial da (co) variância para as observações iid.

Eu olhei para o código para ambos lm()e lm.fit()e nem utilizo var(), mas lm.fit()passa esses dados para o código C compilado ( z <- .Call(C_Cdqrls, x, y, tol, FALSE)) ao qual não tenho acesso.

Questão

Alguém pode explicar por que R está dando resultados diferentes? Mesmo se houver uma diferença no uso da amostra versus variação populacional, por que as estimativas do coeficiente diferem?

Dados

Ajuste uma linha para prever o tamanho do sapato da série na escola.

# 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 Ben Bolker mostrou, às vezes parece que os professores cometem erros. Parece que os cálculos de R estão corretos. Moral da história: não acredite em algo apenas porque um professor diz que é verdade. Verifique você mesmo!

post-hoc
fonte
2
Verifique novamente mod.m=8/3. Porque se você definir mod.m=2.5714, eles parecerão idênticos.
Stat
2
Os coeficientes mod.m = 8/3 e mod.b = -1 não são computados em nenhum lugar dos comentários, pelo que entendi, por isso não é óbvio. Como o @Stat comenta acima, o erro parece estar no modo de computação.
Juho Kokkala
2
É importante ter em mente que qualquer um pode cometer erros - seu professor, você, atendentes aqui, os programadores do R - qualquer um. Portanto, ao tentar descobrir onde podem estar os erros quando as coisas discordam, considere quantas outras pessoas estão verificando cada coisa. No caso da lmfunção em R, literalmente dezenas de milhares de pessoas verificaram os resultados comparando-os com outras coisas, e a saída de lmé comparada com exemplos conhecidos sempre que algo muda no código. Com respostas aqui, pelo menos algumas pessoas provavelmente verificarão (sua pergunta foi examinada 29 vezes).
Glen_b -Reinstate Monica
1
@Glen_b Na verdade, seu ponto é a razão pela qual eu vim aqui para perguntar. Eu não conseguia entender como R poderia estar errado em um cálculo tão básico, mas não conseguia descobrir por que eles eram diferentes. Eu bisbilhotei o código fonte. Mas, no final, o erro estava no último lugar que pensei em procurar, principalmente porque a parte do cálculo está nos limites do meu conhecimento. Eu aprendi muito com a resposta!
post-hoc
2
Sim, é importante tentar descobrir por que eles diferem; faz sentido perguntar aqui se você não consegue resolver isso. Eu estava tentando sugerir por que o último lugar que você considerou poderia ter sido um dos primeiros lugares a procurar. Fui pego fazendo alterações de última hora "simplificadoras" em exemplos em uma ou duas ocasiões.
Glen_b -Reinstar Monica

Respostas:

25

Parece que o autor cometeu um erro matemático em algum lugar.

Se você expandir o desvio da soma dos quadrados

S = b 2 + 2 b m + m 2 + 1 - 2 b - 2 m + b 2 + 4 b

S=((b+m)-1)2+((b+2m)-5)2+((b+4m)-9)2
você obtém
S=b2+2bm+m2+1-2b-2m+b2+4bm+4m2+25-10b-20m+b2+8bm+16m2+81-18b-72m

que reduz para que é igual à expressão do autor, exceto o termo constante, que não importa de qualquer maneira).

3b2+14bm+21m2+107-30b-94m

Agora precisamos tentar minimizar isso definindo os derivados de com respeito à e a zero e resolvendo o sistema. Sbm

dS/db=6b+14m-303b+7m-15=0 0
dS/dm=14b+42.m-947b+21m-47=0 0

Resolver

b=(15-7m)/30 0=7(15-7m)/3+21m-4747-35=(-49./3+21)mm=(47-35)/(21-49./3)=18/7

R diz que este é realmente 2,571429 ...

Baseado neste link, isso parece ser de um curso Coursera ...? Talvez tenha havido uma transcrição incorreta dos dados em algum lugar?

A outra maneira independente de fazer esse cálculo é saber que a inclinação de regressão estimada é igual à soma dos produtos cruzados ( ) dividido pela soma dos quadrados ( ).(y-y¯)(x-x¯)(x-x¯)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

Se você pensa que se o tamanho do sapato fosse vez de , a inclinação chegaria a 8/3 ...{1,11/3,9}{1,5,9}

Ben Bolker
fonte
2
Uau. Sim você está certo. É de um curso Coursera e é do vídeo, não da transcrição. Então, acho que ele o simplificou para tornar os cálculos mais simples para o vídeo e não esperava que alguém tentasse repeti-lo. Por acaso foi o primeiro vídeo que eu vi, então tentei acompanhar. É claro que eu preciso aprimorar meus conhecimentos em matemática. Acho que encontrei o erro. O termo constante, que você diz que não importa, é provavelmente o valor correto que, através de seus cálculos. Examinarei sua resposta mais algumas vezes para me ensinar. Eu realmente gostei disso!
post-hoc
Não acho que o termo constante jogue fora os cálculos. Não afetará as estimativas da inclinação e da interceptação (desaparece quando usamos a derivada), apenas as estimativas do SSQ residual / desvio padrão.
precisa