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

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: $\sigma^2$ $S^2$ lm()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

Podwójna kontrola mod.m=8/3. Ponieważ jeśli ustawisz mod.m=2.5714, będą one wyglądać identycznie.

— Stat

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

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

@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

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

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 + 2 m) - 5)^{2} + ((b + 4 m) - 9)^{2}

$S = ((b+m)-1)^2+ ((b+2m)-5)^2 + ((b+4m)-9)^2$ otrzymasz

\begin{aligned} S = & b^{2} + 2 b m + m^{2} + 1 - 2 b - 2 m \\ + & b^{2} + 4 b m + 4 m^{2} + 25 - 10 b - 20 m \\ + & b^{2} + 8 b m + 16 m^{2} + 81 - 18 b - 72 m \end{aligned}

$\begin{split} S = & b^2+2 b m+ m^2 + 1 - 2 b - 2 m \\ + & b^2+4 b m+ 4 m^2 + 25 - 10 b -20 m \\ + & b^2+8 b m+16 m^2 + 81 - 18 b -72 m \end{split}$

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).

3 b^{2} + 14 b m + 21 m^{2} + 107 - 30 b - 94 m

$3 b^2 + 14 b m + 21 m^2 + 107 - 30 b - 94 m$

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

d S / d b = 6 b + 14 m - 30 \to 3 b + 7 m - 15 = 0

$dS/db = 6 b + 14 m -30 \to 3 b +7 m-15 = 0$

d S / d m = 14 b + 42 m - 94 \to 7 b + 21 m - 47 = 0

$dS/dm = 14 b +42 m -94 \to 7 b + 21 m -47 = 0$

Rozwiązać

\begin{aligned} b & = (15 - 7 m) / 3 \\ 0 & = 7 (15 - 7 m) / 3 + 21 m - 47 \\ 47 - 35 & = (- 49 / 3 + 21) m \\ m & = (47 - 35) / (21 - 49 / 3) = 18 / 7 \end{aligned}

$\begin{split} b & = (15-7m)/3 \\ 0 & = 7 (15-7m)/3 + 21 m-47 \\ 47 - 35 & = (-49/3 + 21) m \\ m & = (47-35)/(21-49/3) = 18/7 \end{split}$

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 ( ). $\sum (y-\bar y) (x-\bar x)$ $\sum (x-\bar 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

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

— Ben Bolker
źródło

Ł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