PTT數位生活區 / R_Language

Re: [問題] 用R做複線性迴歸方程分析

看板R_Language作者 (攸藍)時間11年前 (2013/09/05 17:48), 編輯推噓10(10推 0噓 14→)
留言24則, 4人參與, 最新討論串2/3 (看更多)
※ 引述《sseeaann (屌哥)》之銘言： : 資料檔我是載入套件faraway 用data(savings)資料檔去分析 : 指令大部分都了解 但可能我迴歸觀念不太好 有些步驟會搞混 : ex:像是離群值要常態假設前刪除還是確定為常態假設再刪除 : 或是刪除變數該用stepwise或是用full model和reduced model的F分配來檢定 : 這些步驟有點不知道用的時機點在哪裡 : 有沒有比較懂迴歸分析的大大可以稍微詳述一下整套迴歸的步驟(利用此資料檔)? : 另外小弟我晚一點會在把我打的程式碼PO上來 目前正在打當中^^ : → sseeaann:我在統計版也有發文 不過都沒人回-.- : → sseeaann:我現在比較想知道的是說 當整體模型檢定出現兩個以上變數 : → sseeaann:T值>0.05 不拒絕Ho 那是否可以直接剔除掉? : → sseeaann:如果直接剔除 是要用哪個方法 偏F檢定還是用stepwise Reg 你想要了解的部分就是regression analysis中 model checking 跟 model selection的部分 我就分這兩個part來回答 一邊附上我的程式做解說 library(faraway) data(savings) # Part I. Model checking # 一般來說，model checking，主要就是error term的常態假設以及等變異假設的確認 # 另外，還有利用leverage or Cook's distance查看有沒有outliers # 在R裡面有很方面的方式去做這三件事，程式如下： lm.fit = lm(sr~pop75+pop15+dpi+ddpi, data=savings) summary(lm.fit) par(mfrow = c(2, 2), oma = c(0, 0, 2, 0)) plot(lm.fit) # a. Normality Checking # 右上角就是確認你的殘差，通常就是直接看圖片大概跟QQ-line (圖中的broken line) # 接近，大概我們都會認為error term符合常態假設，至於檢定，就我的經驗來說， # 我實在不認為Goodness of fit檢定有良好的Power..... # 常見的檢定像是：Shapiro-Wilk test of normality # R當然有相對的程式如下： shapiro.test(resid(lm.fit)) # H_0: Normality # 或是你要自己畫QQ plot，我這提供一個範本 library(lattice) qqmath(~ res, data = data.frame(res=resid(lm.fit)), distribution = qnorm, aspect="xy",prepanel = prepanel.qqmathline, panel = function(x, ...) { panel.qqmathline(x, ...) panel.qqmath(x, ...) }) # b. Equal Variance # 左上角的圖就是用來check這個假設，一般來說，不要有特別的pattern # 像是漸小，漸大，先小後大等等，或是有二次曲線出現等 # 通常就認為假設是OK的，我一樣不喜歡等變異假設 # 一樣提供R的等變異假設檢定，在lmtest這個package裡面有 library(lmtest) gqtest(sr~pop75+pop15+dpi+ddpi, data=savings) # 未測試，直接google的 # equal variance圖的其他畫法 xyplot(res~fit, data=data.frame(res=resid(lm.fit), fit=fitted(lm.fit))) # c. Outliers checking # 通常就看右下角跟左下角兩張圖 # 左下角那張 通常不要超過sqrt(3)=1.732就可以了 # 右下角那張看leverage跟Cook's distance，但是至今我還是不太會看... # 我記得我之前學到的是通常不要超過average leverage的0.7倍吧 # 不過R會自動幫你註記他認為有問題的點，機制我就不知道了 # 其他畫法 xyplot(stdres~fit,data=data.frame(fit=fitted(lm.fit),stdres=rstandard(lm.fit))) # d. others # 獨立性檢定就很少再用，因為通常都只能檢驗特定期數，很少用 # R code如下： library(car) durbinWatsonTest(lm.fit) # Part II. Model selection # 這個part就很複雜了，因為還要牽扯到每一個應用層面的關係 # 通常不同的應用層面，有不同的選擇標準以及方式 # 常見的有stepwise, backward, forward # 標準有R^2, AIC, BIC, p-value等 # 通常先進行這個Part，才進行model checking的部分 # 我這裡是因為它比較複雜先放在後面說 # 常見的code像是： step(lm.fit) step(lm.fit, direction="backward") step(lm.fit, direction="forward") library(MASS) stepAIC(lm.fit) # 共線性的檢驗 vif(lm.fit) # 如果大於5，就是有共線性的存在 # 便要去確定covariate的correlation matrix # 通常不建議有超過0.7相關係數的變數存在 # 其他解決方案則有 lasso 跟 ridge等方式 # 至於留存問題，通常根據你對資料的了解而決定 # 但是這個部分，真的不是我在這裡可以說完的，我就不獻醜了.... # 有錯煩請指教，獻醜了，謝謝 剛有google到下面兩個網址，你可以去瞧瞧 http://www.statmethods.net/stats/regression.html http://www.statmethods.net/stats/rdiagnostics.html -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 218.164.79.227
推
09/05 18:00, , 1F
推c大的用心解說(原來可以用lattice畫Q-Q plot啊~)
09/05 18:00, 1F
→
09/05 18:04, , 2F
為了弄出漂亮的東西 還有他比較方便 哈哈
09/05 18:04, 2F
→
09/05 18:09, , 3F
比較方便於畫其他分配...qqplot這個function不愛用
09/05 18:09, 3F
推
09/05 19:50, , 4F
推
09/05 19:50, 4F
→
09/05 19:50, , 5F
Btw, 像這類文章，如果用knitr套件的Rmd寫，排版會好很多
09/05 19:50, 5F
去用了一下，不過似乎無法編譯中文，不論我的input是utf8 or chi-big5...
→
09/05 21:48, , 6F
謝謝版大建議QQ
09/05 21:48, 6F
推
09/05 21:55, , 7F
感謝c大的解說 對小弟非常有用 另外圖上標記點是離群值?
09/05 21:55, 7F
→
09/05 21:56, , 8F
因為我做的圖適用plot()指令和qqnorm()指令來做 不會特別
09/05 21:56, 8F
→
09/05 21:58, , 9F
標記離群值 得自己去找 能再請教c大模型的處理共線性指令
09/05 21:58, 9F
處理共線性最快的方法就是用 vif的指令 不知道在MASS還是car裡面就有 ※ 編輯: celestialgod 來自: 218.164.79.227 (09/05 22:14)
推
09/05 22:23, , 10F
好的 我會試試看^^
09/05 22:23, 10F
推
09/05 23:40, , 11F
另外請教一下c大 預測和信賴區間的指令您知道怎麼表示嗎?
09/05 23:40, 11F
# confidence interval for fitted value predict(lm.fit, interval="confidence", level=0.95) # confidence interval for prediction pre.data = matrix(c(30,10,700,4,40,5,1000,2),2,,TRUE) rownames(pre.data) = c("A_region", "B_region") newdata=data.frame(pre.data) names(newdata)=names(savings)[2:5] predict(lm.fit, newdata=newdata ,interval="prediction", level=0.95) # To ger more information by R>?predict.lm
→
09/05 23:49, , 12F
簡單說明第二塊，就是我創了兩個地區的人口等變項
09/05 23:49, 12F
→
09/05 23:49, , 13F
predict要求新資料必須有一樣的變數名稱
09/05 23:49, 13F
→
09/05 23:50, , 14F
最後調整interval即是預測區間
09/05 23:50, 14F
# 此外，迴歸係數的信賴區間可由下列code獲得 confint(lm.fit)
推
09/06 00:00, , 15F
請問一下 我發現var.test()指令好像也能得到95%信賴區間
09/06 00:00, 15F
→
09/06 00:03, , 16F
另外刪除離群值的模型還是有可能常態假設不成立或
09/06 00:03, 16F
→
09/06 00:03, , 17F
等變異不相等? 假如刪除前的模型常態假設和等變異相等
09/06 00:03, 17F
# 刪除前後的不同，這我可能無法回答你 # 但是就我所學來說，刪除離群值這件事真的很少做... # 我是念純統計的，我們老師通常都說資料最大 # 如果不是記錄錯誤的話，通常都會納入考量... # 或是influence point才會考慮刪除。
→
09/06 00:04, , 18F
另外真的感謝c大 幫了小弟很大的忙^^
09/06 00:04, 18F
推
09/06 00:25, , 19F
encoding issue 滿煩的，我自己是用utf8編譯都沒問題
09/06 00:25, 19F
推
09/06 00:26, , 20F
09/06 00:26, 20F
→
09/06 00:37, , 21F
windows下這個議題更麻煩XDDD
09/06 00:37, 21F
推
09/06 00:50, , 22F
c大有用過var.test()這指令求等變異嗎? 還像也會跑出
09/06 00:50, 22F
→
09/06 00:51, , 23F
95%的信賴區間 請問這是哪個值的信賴區間@@?
09/06 00:51, 23F
variance ratio的信賴區間 我看了一下，那個只能用在category的covariate... 連續型的covariate不能用那個檢驗等變異.... ※ 編輯: celestialgod 來自: 218.164.79.227 (09/06 01:15)
推
09/06 02:27, , 24F
了解 原來是我搞錯了><
09/06 02:27, 24F
‣ 更多 celestialgod 的文章
文章代碼(AID): #1IA5Dorj (R_Language)
R_Language 近期熱門文章
更多 近期熱門文章 >>
PTT數位生活區 即時熱門文章
更多 即時熱門文章 >>
‣ 更多 celestialgod 的文章
文章代碼(AID): #1IA5Dorj (R_Language)