題:
為什麼R的lm()返回的係數估計與我的教科書不同?
post-hoc
2014-06-25 03:44:31 UTC
view on stackexchange narkive permalink

背景

我正在嘗試了解關於擬合模型的課程中的 first 示例(因此,這似乎很簡單)。我已經手工完成了計算,並且它們與示例匹配,但是當我在R中重複計算時,模型係數不可用。我以為差異可能是由於教科書使用總體方差($ \ sigma ^ 2 $)而R可能是使用樣本方差($ S ^ 2 $),但我看不到這些在計算中的位置。例如,如果 lm()在某處使用 var(),則 var()上的幫助部分將指出:

使用分母n-1給出iid(協)方差的無偏估計觀察。

我已經查看了 lm() lm.fit()的代碼,但都沒有使用 var(),但 lm.fit()會將數據傳遞到編譯的Ccode( z <- .Call(C_Cdqrls,x,y,tol,FALSE)代碼>)我無法訪問。

問題

有人可以解釋為什麼R給出不同的結果嗎?即使使用樣本方差與總體方差有差異,為什麼係數估計也不同?

數據

適合一條線以根據學校年級預測鞋子的大小。

#模型datamod.dat <- read.table(text ='grade shoe 1 1 2 5 4 9',header = T);#meanmod.mu <-均值(mod.dat $ shoe);#變異性mod.var <- sum((mod.dat $ shoe-mod.mu)^ 2)#來自教科書的模型係數mod.m <-8 / 3; mod.b <- -1;#預測值(1.666667 4.333333 9.666667)mod.man.pred <- mod.dat $ grade * mod.m + mod.b;#殘差(-0.6666667 0.6666667 -0.6666667)mod.man.resid <-(mod .dat $ shoe-mod.man.pred)#剩餘方差(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;#但lm()給出不同的結果:摘要(lm(鞋〜等級,數據= mod .dat))
  Call:lm(公式=鞋子〜成績,數據= mod.dat)殘差:1 2 3- 0.5714 0.8571 -0.2857係數:估計標準誤差t值Pr(> | t |)(攔截)-1.0000 1.3093 -0.764 0.585級2.5714 0.4949 5.196 0.121殘餘標準誤差:1個自由度上為1.069多個R平方:0.9643,調整後R平方:0.9286 F統計:在1對1和1對DF上,p值:0.121  

Edit

Ben Bolker所示,看起來老師在犯錯誤有時。看來R的計算是正確的。這個故事的寓意:不要僅僅因為老師說的是真實而相信某事。親自驗證!

您可以編輯問題以同時包含您手工進行的計算以匹配教科書的結果嗎?
@JuhoKokkala他們已經在註釋中。我試圖節省空間,但可能並不明顯。
仔細檢查`mod.m = 8 / 3`。因為如果您設置`mod.m = 2.5714`,那麼它們似乎是相同的。
據我所知,註釋中的任何地方都沒有計算係數mod.m = 8/3和mod.b = -1,所以這並不明顯。正如上面的@Stat所評論的,該錯誤似乎出在計算mod.m中。
@Stat在文本中列出了值'8/3'和`-1'作為通過解決以下等式而產生最小誤差的值:$ 21m ^ 2 + 14mb + 3b ^ 2-94m-30b + 81 $。老實說,我不完全理解這一部分,但也許您是在說這些估計值是出於教學目的而四捨五入的?
@JuhoKokkala您是正確的。我只是接受了教科書中的那個價值,因為我真的不明白那部分。係數應被認為是最小化線路誤差的最佳解決方案,但也許將其舍入以使其他計算更容易?他沒有解釋他如何到達“ 8/3”和“ -1”,老實說,我也不知道一個人怎麼到達那裡!我應該相信trust的結果嗎?我問這個問題是因為我試圖重新分析其他人的結果,並且我遇到了同樣的問題,即ℝ沒有給出與論文作者相同的估計。
重要的是要記住,“任何人”都可能會犯錯誤-您的老師,您,此處的回答者,R程序員-任何人。因此,當試圖找出當事情分歧時錯誤可能出在哪裡時,請考慮有多少其他人正在檢查每件事。在R中的lm函數的情況下,實際上有成千上萬的人通過將結果與其他事物進行比較來檢查結果,並且每次代碼中的任何更改都對照已知示例檢查lm的輸出。在這裡有了答案,至少有幾個人可能會檢查(您的問題已被查看過29次)。
@Glen_b您的意思實際上就是我來這裡詢問的原因。我不明白在這樣的基本計算中R怎麼可能是錯誤的,但我不知道為什麼它們不一樣。我監聽了源代碼。但是最後,錯誤出現在我想看的最後一個地方,主要是因為演算部分處於我所知的範圍之內。我從答案中學到了很多東西!
是的,重要的是設法弄清楚它們為何不同。在這裡問是否無法解決是很有意義的。我試圖提出一個建議,為什麼您考慮的最後一個地方可能卻成為了第一個出現的地方。我自己一次或兩次對示例進行最後的“簡化”更改而被吸引住了。
一 回答:
Ben Bolker
2014-06-25 05:20:07 UTC
view on stackexchange narkive permalink

看來作者在某個地方犯了數學錯誤。

如果擴大平方和偏差

$$ S =((b + m)-1 )^ 2 +(((b + 2m)-5)^ 2 +((b + 4m)-9)^ 2 $$您將獲得$$ \ begin {split} S = & b ^ 2 + 2 b m + m ^ 2 + 1-2 b-2 m \\ + & b ^ 2 + 25-10 b -20 m \\ + & b ^ 2 + 8 b m + 16 m ^ 2 + 81 -18 b -72 m \ end {split} $$

減少到$$ 3 b ^ 2 + 14 bm + 21 m ^ 2 + 107-30 b-94 m $$

現在,我們需要通過針對$ b $和$ m設置$ S $的導數來嘗試最小化此常量。 $到零並求解系統。$$ dS / db = 6 b + 14 m -30 \到3 b +7 m-15 = 0 $$$$ ds / dm = 14 b +42 m -94 \到7 b + 21 m -47 = 0 $$

解決

$$ \開始{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說這確實是2.571429 ...

基於此鏈接,這似乎來自Coursera課程...?也許某處的數據有錯誤的轉錄?

進行此計算的另一種獨立方法是,知道估計的回歸斜率等於叉積($ \ sum(y -\ bar y)(x- \ bar x)$)除以平方和($ \ sum(x- \ bar x)^ 2 $)。

  g <- c (1,2,4)g0 <- g-平均值(g)s <- c(1,5,9)s0 <- s-平均值sum(g0 * s0)/(sum(g0 ^ 2) )## [1] 2.571429  

如果認為鞋子的尺碼是$ \ {1,11 / 3,9 \} $而不是$ \ {1,5,9 \} $然後斜率變為8/3 ...

哇。你是對的。它來自Coursera課程,來自視頻,而不是轉錄。因此,我猜他是簡化了視頻的計算,卻沒想到有人會重複嘗試。它恰好是我看到的第一個視頻,所以我嘗試跟進。顯然,我需要提高數學水平。我認為雖然發現了錯誤。您說的無關緊要的常數項可能是通過他的計算得出的正確值。我將再幾次瀏覽您的答案以自學。我真的很感激!
我認為常數項不會影響計算。它不會影響斜率和截距的估計(當我們採用導數時,它會消失),只會影響殘餘SSQ /標準偏差的估計。


該問答將自動從英語翻譯而來。原始內容可在stackexchange上找到,我們感謝它分發的cc by-sa 3.0許可。
Loading...