library(faraway)
data(savings)
#建立複線性迴歸模型
mdl<-lm(sr~pop15+pop75+dpi+ddpi,savings)
#使用指令coefficients()得到整體迴歸模型迴歸係數beta，以savings資料為例:
#sr=28.56609+(-0.46119)*pop15+(-1.69150)*pop75+(-0.00034)*dpi+(0.40969)*ddpi
round(coefficients(mdl),5)
 
#檢定各變數的共線性
vif(mdl)
#這邊有發現兩個變數超過5 分別是pop75、pop15 
#這兩個變數要忽略? 還是先偏F檢定再對偏F檢定後的模型做指令vif()?
 
#使用指令summary()可以彙總物件mdl的內容，了解上述模型是否可簡化，及探討模型的解釋能力。
summary(mdl)
#以savings資料為例，F-statistic: 5.756 on 4 and 45 DF，p-value: 0.0007904知道引入解釋變數是有統計意義。
這個模型mdl的能解釋y的變異達 33.85%，後面須看看配適殘差圖呈現的型態及觀察離群值是否有影響。
再由beta係數的彙總資料發現pop75的p-value=0.125530 > 0.05和pop75的p-value=0.719173 > 0.05 皆不拒絕Ho可以考慮把pop75或dpi拿掉。
#以上的結果要成立必須是，residuals必須是滿足常態假設，且變異數為常數。
#先利用偏F檢定判斷上面兩個不拒絕Ho之變數是否該忽略。
#設一個忽略掉pop75、dpi這個變數的新的模型為mdl1。
 
mdl1<-lm(sr~pop15+ddpi,savings)
#再使用anova()指令對mdl和mdl1兩個模型做變異數分析，檢定pop75和dpi可以同時忽略(Ho : pop75=dpi=0)。
anova(mdl1,mdl)
#由變異數分析資料發現的p-value=0.19> 0.05 所以不拒絕Ho，考慮用Model1。
#這邊我省略了model selection是因為我沒打算要加其他變數 想說用偏F檢定篩選掉可忽略的變數即可
 
#殘差VS配適值關係圖
plot(fitted(mdl1),residuals(mdl1))
mean(fitted(mdl1))
#使用var.test()指令檢定該模型mdl1的變異數是否為常數，配適值範圍可用模型之配適值之平均值。
var.test(residuals(mdl1)[fitted(mdl1)>9.671],residuals(mdl1)[fitted(mdl1)<9.671])
#使用qqnorm()指令跑出QQ-norm圖。
qqnorm(residuals(mdl1))
qqline(residuals(mdl1))
#使用shapiro.test()指令檢定此模型mdl1之常態假設是否成立。
shapiro.test(residuals(mdl1))
#常態假設以及變異數檢定皆成立，最後再檢查是否有離群值存在，
 
#載入套件car，使用outlierTest()指令找離群值。
library(car)
outlierTest(mdl1)
#可以找到Zambia這個離群值
 
library(car)
influencePlot(mdl1)
#但使用influencePlot()會找到兩個影響點Zambia和Libya
#這邊是要把Zambia刪除後的模型再做檢定 還是兩個影響點都須刪除?
 
#另外有沒有哪些重要的步驟我卻忽略了@@? 目標是找到最適模型即可! 有的話麻煩c大補充一下><
 

bGlicmFyeShmYXJhd2F5KQpkYXRhKHNhdmluZ3MpCiPlu7rnq4vopIfnt5rmgKfov7TmrbjmqKHlnosKbWRsPC1sbShzcn5wb3AxNStwb3A3NStkcGkrZGRwaSxzYXZpbmdzKQoj5L2/55So5oyH5LukY29lZmZpY2llbnRzKCnlvpfliLDmlbTpq5Tov7TmrbjmqKHlnovov7Tmrbjkv4LmlbhiZXRh77yM5Lulc2F2aW5nc+izh+aWmeeCuuS+izoKI3NyPTI4LjU2NjA5KygtMC40NjExOSkqcG9wMTUrKC0xLjY5MTUwKSpwb3A3NSsoLTAuMDAwMzQpKmRwaSsoMC40MDk2OSkqZGRwaQpyb3VuZChjb2VmZmljaWVudHMobWRsKSw1KQoKI+aqouWumuWQhOiuiuaVuOeahOWFsee3muaApwp2aWYobWRsKQoj6YCZ6YKK5pyJ55m854++5YWp5YCL6K6K5pW46LaF6YGONSDliIbliKXmmK9wb3A3NeOAgXBvcDE1IAoj6YCZ5YWp5YCL6K6K5pW46KaB5b+955WlPyDpgoTmmK/lhYjlgY9G5qqi5a6a5YaN5bCN5YGPRuaqouWumuW+jOeahOaooeWei+WBmuaMh+S7pHZpZigpPwoKI+S9v+eUqOaMh+S7pHN1bW1hcnkoKeWPr+S7peW9mee4veeJqeS7tm1kbOeahOWFp+Wuue+8jOS6huino+S4iui/sOaooeWei+aYr+WQpuWPr+ewoeWMlu+8jOWPiuaOouiojuaooeWei+eahOino+mHi+iDveWKm+OAggpzdW1tYXJ5KG1kbCkKI+S7pXNhdmluZ3Pos4fmlpnngrrkvovvvIxGLXN0YXRpc3RpYzogNS43NTYgb24gNCBhbmQgNDUgREbvvIxwLXZhbHVlOiAwLjAwMDc5MDTnn6XpgZPlvJXlhaXop6Pph4vorormlbjmmK/mnInntbHoqIjmhI/nvqnjgIIK6YCZ5YCL5qih5Z6LbWRs55qE6IO96Kej6YeLeeeahOiuiueVsOmBlCAzMy44NSXvvIzlvozpnaLpoIjnnIvnnIvphY3pganmrpjlt67lnJblkYjnj77nmoTlnovmhYvlj4rop4Dlr5/pm6LnvqTlgLzmmK/lkKbmnInlvbHpn7/jgIIK5YaN55SxYmV0YeS/guaVuOeahOW9mee4veizh+aWmeeZvOePvnBvcDc155qEcC12YWx1ZT0wLjEyNTUzMCA+IDAuMDXlkoxwb3A3NeeahHAtdmFsdWU9MC43MTkxNzMgPiAwLjA1IOeahuS4jeaLkue1lUhv5Y+v5Lul6ICD5oWu5oqKcG9wNzXmiJZkcGnmi7/mjonjgIIKI+S7peS4iueahOe1kOaenOimgeaIkOeri+W/hemgiOaYr++8jHJlc2lkdWFsc+W/hemgiOaYr+a7v+i2s+W4uOaFi+WBh+iore+8jOS4lOiuiueVsOaVuOeCuuW4uOaVuOOAggoj5YWI5Yip55So5YGPRuaqouWumuWIpOaWt+S4iumdouWFqeWAi+S4jeaLkue1lUhv5LmL6K6K5pW45piv5ZCm6Kmy5b+955Wl44CCCiPoqK3kuIDlgIvlv73nlaXmjolwb3A3NeOAgWRwaemAmeWAi+iuiuaVuOeahOaWsOeahOaooeWei+eCum1kbDHjgIIKCm1kbDE8LWxtKHNyfnBvcDE1K2RkcGksc2F2aW5ncykKI+WGjeS9v+eUqGFub3ZhKCnmjIfku6TlsI1tZGzlkoxtZGwx5YWp5YCL5qih5Z6L5YGa6K6K55Ww5pW45YiG5p6Q77yM5qqi5a6acG9wNzXlkoxkcGnlj6/ku6XlkIzmmYLlv73nlaUoSG8gOiBwb3A3NT1kcGk9MCnjgIIKYW5vdmEobWRsMSxtZGwpCiPnlLHorornlbDmlbjliIbmnpDos4fmlpnnmbznj77nmoRwLXZhbHVlPTAuMTk+IDAuMDUg5omA5Lul5LiN5ouS57WVSG/vvIzogIPmha7nlKhNb2RlbDHjgIIKI+mAmemCiuaIkeecgeeVpeS6hm1vZGVsIHNlbGVjdGlvbuaYr+WboOeCuuaIkeaykuaJk+eul+imgeWKoOWFtuS7luiuiuaVuCDmg7PoqqrnlKjlgY9G5qqi5a6a56+p6YG45o6J5Y+v5b+955Wl55qE6K6K5pW45Y2z5Y+vCgoj5q6Y5beuVlPphY3pganlgLzpl5zkv4LlnJYKcGxvdChmaXR0ZWQobWRsMSkscmVzaWR1YWxzKG1kbDEpKQptZWFuKGZpdHRlZChtZGwxKSkKI+S9v+eUqHZhci50ZXN0KCnmjIfku6TmqqLlrproqbLmqKHlnottZGwx55qE6K6K55Ww5pW45piv5ZCm54K65bi45pW477yM6YWN6YGp5YC856+E5ZyN5Y+v55So5qih5Z6L5LmL6YWN6YGp5YC85LmL5bmz5Z2H5YC844CCCnZhci50ZXN0KHJlc2lkdWFscyhtZGwxKVtmaXR0ZWQobWRsMSk+OS42NzFdLHJlc2lkdWFscyhtZGwxKVtmaXR0ZWQobWRsMSk8OS42NzFdKQoj5L2/55SocXFub3JtKCnmjIfku6Tot5Hlh7pRUS1ub3Jt5ZyW44CCCnFxbm9ybShyZXNpZHVhbHMobWRsMSkpCnFxbGluZShyZXNpZHVhbHMobWRsMSkpCiPkvb/nlKhzaGFwaXJvLnRlc3QoKeaMh+S7pOaqouWumuatpOaooeWei21kbDHkuYvluLjmhYvlgYfoqK3mmK/lkKbmiJDnq4vjgIIKc2hhcGlyby50ZXN0KHJlc2lkdWFscyhtZGwxKSkKI+W4uOaFi+WBh+ioreS7peWPiuiuiueVsOaVuOaqouWumueahuaIkOeri++8jOacgOW+jOWGjeaqouafpeaYr+WQpuaciemboue+pOWAvOWtmOWcqO+8jAoKI+i8ieWFpeWll+S7tmNhcu+8jOS9v+eUqG91dGxpZXJUZXN0KCnmjIfku6Tmib7pm6LnvqTlgLzjgIIKbGlicmFyeShjYXIpCm91dGxpZXJUZXN0KG1kbDEpCiPlj6/ku6Xmib7liLBaYW1iaWHpgJnlgIvpm6LnvqTlgLwKCmxpYnJhcnkoY2FyKQppbmZsdWVuY2VQbG90KG1kbDEpCiPkvYbkvb/nlKhpbmZsdWVuY2VQbG90KCnmnIPmib7liLDlhanlgIvlvbHpn7/pu55aYW1iaWHlkoxMaWJ5YQoj6YCZ6YKK5piv6KaB5oqKWmFtYmlh5Yiq6Zmk5b6M55qE5qih5Z6L5YaN5YGa5qqi5a6aIOmChOaYr+WFqeWAi+W9semfv+m7numDvemgiOWIqumZpD8KCiPlj6blpJbmnInmspLmnInlk6rkupvph43opoHnmoTmraXpqZ/miJHljbvlv73nlaXkuoZAQD8g55uu5qiZ5piv5om+5Yiw5pyA6YGp5qih5Z6L5Y2z5Y+vISDmnInnmoToqbHpurvnhalj5aSn6KOc5YWF5LiA5LiLPjwK