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

[celestialgod (攸藍)]

※ 引述《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
推 diplazium:推c大的用心解說(原來可以用lattice畫Q-Q plot啊~)
→ celestialgod:為了弄出漂亮的東西 還有他比較方便 哈哈
→ celestialgod:比較方便於畫其他分配...qqplot這個function不愛用
推 Wush978:推
→ Wush978:Btw, 像這類文章，如果用knitr套件的Rmd寫，排版會好很多
去用了一下，不過似乎無法編譯中文，不論我的input是utf8 or chi-big5...

→ celestialgod:謝謝版大建議QQ
推 sseeaann:感謝c大的解說 對小弟非常有用 另外圖上標記點是離群值?
→ sseeaann:因為我做的圖適用plot()指令和qqnorm()指令來做 不會特別
→ sseeaann:標記離群值 得自己去找 能再請教c大模型的處理共線性指令
處理共線性最快的方法就是用 vif的指令 不知道在MASS還是car裡面就有
※ 編輯: celestialgod    來自: 218.164.79.227       (09/05 22:14)
推 sseeaann:好的 我會試試看^^
推 sseeaann:另外請教一下c大 預測和信賴區間的指令您知道怎麼表示嗎?
# 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
→ celestialgod:簡單說明第二塊，就是我創了兩個地區的人口等變項
→ celestialgod:predict要求新資料必須有一樣的變數名稱
→ celestialgod:最後調整interval即是預測區間
# 此外，迴歸係數的信賴區間可由下列code獲得
confint(lm.fit)
推 sseeaann:請問一下 我發現var.test()指令好像也能得到95%信賴區間
→ sseeaann:另外刪除離群值的模型還是有可能常態假設不成立或
→ sseeaann:等變異不相等? 假如刪除前的模型常態假設和等變異相等
# 刪除前後的不同，這我可能無法回答你
# 但是就我所學來說，刪除離群值這件事真的很少做...
# 我是念純統計的，我們老師通常都說資料最大
# 如果不是記錄錯誤的話，通常都會納入考量...
# 或是influence point才會考慮刪除。
→ sseeaann:另外真的感謝c大 幫了小弟很大的忙^^
推 Wush978:encoding issue 滿煩的，我自己是用utf8編譯都沒問題
推 Wush978:http://rpubs.com/wush978
→ celestialgod:windows下這個議題更麻煩XDDD
推 sseeaann:c大有用過var.test()這指令求等變異嗎? 還像也會跑出
→ sseeaann:95%的信賴區間 請問這是哪個值的信賴區間@@?
variance ratio的信賴區間
我看了一下，那個只能用在category的covariate...
連續型的covariate不能用那個檢驗等變異....
※ 編輯: celestialgod    來自: 218.164.79.227       (09/06 01:15)
推 sseeaann:了解 原來是我搞錯了><