Mengapa R's lm () mengembalikan estimasi koefisien yang berbeda dari buku teks saya?

13

Latar Belakang

Saya mencoba untuk memahami contoh pertama dalam kursus model pas (jadi ini mungkin tampak sangat sederhana). Saya sudah melakukan perhitungan dengan tangan dan mereka cocok dengan contoh, tetapi ketika saya ulangi dalam R, koefisien model tidak aktif. Saya pikir perbedaannya mungkin karena buku teks menggunakan varians populasi ( σ2 ) sedangkan R mungkin menggunakan varians sampel ( S2 ), tapi saya tidak bisa melihat di mana ini digunakan dalam perhitungan. Misalnya, jika lm()menggunakan var()suatu tempat, bagian bantuan pada var()catatan:

Penyebut n - 1 digunakan yang memberikan penaksir yang tidak bias dari varians (co) untuk pengamatan iid.

Saya telah melihat kode untuk keduanya lm()dan lm.fit()dan tidak menggunakan var(), tetapi lm.fit()melewati data itu untuk mengkompilasi kode C ( z <- .Call(C_Cdqrls, x, y, tol, FALSE)) yang saya tidak punya akses.

Pertanyaan

Adakah yang bisa menjelaskan mengapa R memberikan hasil yang berbeda? Bahkan jika ada perbedaan dalam menggunakan varians sampel vs populasi, mengapa estimasi koefisien berbeda?

Data

Sesuaikan garis untuk memprediksi ukuran sepatu dari kelas di sekolah.

# 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

Edit

Seperti yang ditunjukkan Ben Bolker , guru terkadang membuat kesalahan. Tampaknya perhitungan R sudah benar. Moral cerita: jangan percaya sesuatu hanya karena seorang guru mengatakan itu benar. Verifikasi sendiri!

r regression self-study lm pasca-hoc
sumber
2
Periksa ulang mod.m=8/3. Karena jika Anda mengatur mod.m=2.5714, maka mereka tampaknya identik.
Stat
2
Koefisien mod.m = 8/3 dan mod.b = -1 tidak dihitung di mana pun di komentar sejauh yang saya mengerti, jadi tidak jelas. Seperti yang dikomentari @Stat di atas, kesalahannya tampaknya ada dalam komputasi mod.m.
Juho Kokkala
2
Penting untuk diingat bahwa siapa pun dapat membuat kesalahan - guru Anda, Anda, penjawab di sini, programmer R - siapa pun. Jadi ketika mencoba mencari tahu di mana kesalahan mungkin terletak ketika hal-hal tidak setuju, pertimbangkan berapa banyak orang yang memeriksa setiap hal. Dalam kasus lmfungsi dalam R, secara harfiah puluhan ribu orang telah memeriksa hasilnya dengan membandingkannya dengan hal-hal lain, dan outputnya lmdiperiksa terhadap contoh-contoh yang diketahui setiap kali ada perubahan kode. Dengan jawaban di sini, setidaknya beberapa orang cenderung memeriksa (pertanyaan Anda telah dilihat sebanyak 29 kali).
Glen_b -Reinstate Monica
1
@ Glen_b Maksud Anda sebenarnya adalah alasan mengapa saya datang ke sini untuk bertanya. Saya tidak bisa mengerti bagaimana R bisa salah dalam perhitungan dasar seperti itu, tetapi saya tidak tahu mengapa mereka berbeda. Saya acara mengintip di sekitar kode sumber. Tetapi pada akhirnya, kesalahannya adalah di tempat terakhir yang saya pikir saya perhatikan, sebagian besar karena bagian kalkulus berada pada batas pengetahuan saya. Saya belajar banyak dari jawabannya!
post-hoc
2
Ya, penting untuk mencoba mencari tahu mengapa mereka berbeda; masuk akal untuk bertanya di sini jika Anda tidak bisa menyelesaikannya. Saya mencoba menyarankan mengapa tempat terakhir yang Anda pertimbangkan mungkin menjadi salah satu tempat pertama yang terlihat. Saya terperangkap dengan membuat perubahan 'penyederhanaan' menit-menit terakhir pada contoh pada satu atau dua kesempatan sendiri.
Glen_b -Reinstate Monica

Jawaban:

25

Sepertinya penulis membuat kesalahan matematis di suatu tempat.

Jika Anda memperluas jumlah penjumlahan kuadrat

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

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

Sbm

dS/db=6b+14m-30→3b+7m-15=0
dS/dm=14b+42m-94→7b+21m-47=0

Memecahkan

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

R mengatakan ini memang 2.571429 ...

Berdasarkan tautan ini, sepertinya ini berasal dari kursus Coursera ...? Mungkin ada salah-transkripsi data di suatu tempat?

∑(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

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

Ben Bolker
sumber
2
Wow. Ya kamu benar. Ini dari kursus Coursera dan dari video, bukan transkripsi. Jadi saya kira dia menyederhanakannya untuk membuat perhitungan lebih sederhana untuk video dan tidak mengharapkan siapa pun untuk mencoba dan mengulanginya. Kebetulan itu adalah video pertama yang saya lihat jadi saya mencoba untuk mengikuti. Sudah jelas bahwa saya perlu meningkatkan keterampilan dalam soal matematika. Saya pikir menemukan kesalahannya. Istilah konstan, yang Anda katakan tidak masalah, mungkin nilai yang benar yang melalui perhitungannya. Saya akan memeriksa jawaban Anda beberapa kali lagi untuk belajar sendiri. Saya sangat menghargai itu!
post-hoc
Saya tidak berpikir istilah konstan akan membuang perhitungan. Itu tidak akan mempengaruhi estimasi kemiringan dan intersep (menghilang ketika kita mengambil turunannya), hanya perkiraan sisa SSQ / standar deviasi.
Ben Bolker