Dlaczego lm () R zwraca inne współczynniki niż mój podręcznik?

13

tło

Próbuję zrozumieć pierwszy przykład w kursie na temat dopasowywania modeli (więc może się to wydawać absurdalnie proste). Obliczenia wykonałem ręcznie i pasują one do przykładu, ale kiedy powtórzę je w R, współczynniki modelu są wyłączone. Myślałem, że różnica może wynikać z tego, że podręcznik używa wariancji populacji ( ), podczas gdy R może używać wariancji próbki ( ), ale nie widzę, gdzie są one stosowane w obliczeniach. Na przykład, jeśli używa gdzieś, sekcja pomocy na notatki:S 2σ2S2lm()var()var()

Stosuje się mianownik n - 1, który daje obiektywny estymator wariancji (ko) dla obserwacji iid.

Mam spojrzał na kod dla obu lm()i lm.fit()i nie wnosić wykorzystania var(), ale lm.fit()przekazuje te dane do skompilowanego kodu C ( z <- .Call(C_Cdqrls, x, y, tol, FALSE)), które nie mają dostępu.

Pytanie

Czy ktoś może wyjaśnić, dlaczego R daje różne wyniki? Nawet jeśli istnieje różnica w stosowaniu wariancji między próbą a populacją, dlaczego szacunki współczynników różnią się?

Dane

Dopasuj linię, aby przewidzieć rozmiar buta na podstawie klasy w szkole.

# 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

Edytować

Jak pokazał Ben Bolker , wygląda na to, że nauczyciele czasami popełniają błędy. Wygląda na to, że obliczenia R są prawidłowe. Morał tej historii: nie wierz w coś tylko dlatego, że nauczyciel mówi, że to prawda. Sprawdź to sam!

post-hoc
źródło
2
Podwójna kontrola mod.m=8/3. Ponieważ jeśli ustawisz mod.m=2.5714, będą one wyglądać identycznie.
Stat
2
Współczynniki mod.m = 8/3 i mod.b = -1 nie są obliczane nigdzie w komentarzach, o ile rozumiem, więc nie jest to oczywiste. Jak komentuje @Stat powyżej, błąd wydaje się związany z obliczaniem mod.m.
Juho Kokkala
2
Ważne jest, aby pamiętać, że każdy może popełniać błędy - twój nauczyciel, ty, tu odpowiadający, programiści R - każdy. Więc próbując dowiedzieć się, gdzie mogą leżeć błędy, gdy coś się nie zgadza, zastanów się, ilu innych ludzi sprawdza każdą rzecz. W przypadku lmfunkcji w R dosłownie dziesiątki tysięcy osób sprawdzało wyniki, porównując je z innymi rzeczami, a wyniki lmsprawdzane są na znanych przykładach za każdym razem, gdy cokolwiek zmienia się w kodzie. Z odpowiedziami tutaj co najmniej kilka osób może sprawdzić (twoje pytanie zostało sprawdzone 29 razy).
Glen_b
1
@Glen_b Twój punkt jest właściwie powodem, dla którego przybyłem tutaj, aby zapytać. Nie mogłem zrozumieć, w jaki sposób R może się mylić przy tak podstawowych obliczeniach, ale nie mogłem zrozumieć, dlaczego były różne. I zdarzało się węszyć wokół kodu źródłowego. Ale w końcu błąd był na ostatnim miejscu, o którym pomyślałem, głównie dlatego, że część rachunku różniczkowego leży na granicy mojej wiedzy. Wiele się jednak nauczyłem z odpowiedzi!
post-hoc
2
Tak, ważne jest, aby spróbować dowiedzieć się, dlaczego się różnią; warto zapytać tutaj, czy nie możesz tego wypracować. Próbowałem zasugerować, dlaczego ostatnie miejsce, które rozważałeś, mogło być jednym z pierwszych miejsc do obejrzenia. Zostałem przyłapany na tym, że w ostatniej chwili sam „upraszczałem” zmiany przykładów.
Glen_b

Odpowiedzi:

25

Wygląda na to, że autor popełnił gdzieś błąd matematyczny.

Jeśli powiększysz odchylenie sumy kwadratów

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

co zmniejsza się do co jest takie samo jak wyrażenie autora, z wyjątkiem stałego wyrażenia, które i tak nie ma znaczenia).

3b2+14bm+21m2+10730b94m

Teraz musimy spróbować to zminimalizować, ustawiając pochodne względem i na zero i rozwiązując układ. Sbm

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

Rozwiązać

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

R mówi, że to rzeczywiście 2.571429 ...

Na podstawie tego linku wydaje się, że pochodzi on z kursu Coursera ...? Może gdzieś była błędna transkrypcja danych?

Innym niezależnym sposobem wykonania tego obliczenia jest wiedza, że ​​oszacowane nachylenie regresji jest równe sumie iloczynów krzyżowych ( ) podzielonej przez sumę kwadratów ( ).(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

Jeśli pomyślimy, że rozmiary butów to zamiast wówczas nachylenie wyniesie 8/3 ...{ 1 , 5 , 9 }{1,11/3,9}{1,5,9}

Ben Bolker
źródło
2
Łał. Tak masz rację. Pochodzi z kursu Coursera i pochodzi z filmu, a nie z transkrypcji. Zgaduję więc, że uprościł to, aby ułatwić obliczenia wideo i nie spodziewał się, że ktoś spróbuje go powtórzyć. To był po prostu pierwszy film, który widziałem, więc starałem się śledzić. Oczywiste jest, że muszę podnieść umiejętności, jeśli chodzi o matematykę. Myślę jednak, że znalazłem błąd. Stały termin, który, jak mówisz, nie ma znaczenia, jest prawdopodobnie poprawną wartością, która wynika z jego obliczeń. Jeszcze kilka razy przejrzę twoją odpowiedź, by się nauczyć. Bardzo to doceniam!
post-hoc
Nie sądzę, aby stały termin obalił obliczenia. Nie wpłynie to na oszacowanie nachylenia i przecięcia (znika, gdy weźmiemy pochodną), tylko szacunki resztkowego SSQ / odchylenia standardowego.
Ben Bolker