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)
#使用指令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。
#殘差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)
#文字的離群值可以照C大的指令刪除了 是我不小心指令有打錯= =+
#但有問題是outlierTest()指令只會跑出一個離群值 而且不一定會是影響點 應該是不能亂動
#有個想法是說利用cook's distance plot找出影響點 再將影響點從資料中刪除 再重新建模再做檢定 不知到是否可行?
cutoff <- 4/((nrow(savings)-length(mdl1$coefficients)-2)) 
plot(mdl1, which=4, cook.levels=cutoff)

bGlicmFyeShmYXJhd2F5KQpkYXRhKHNhdmluZ3MpCiPlu7rnq4vopIfnt5rmgKfov7TmrbjmqKHlnosKbWRsPC1sbShzcn5wb3AxNStwb3A3NStkcGkrZGRwaSxzYXZpbmdzKQoj5L2/55So5oyH5LukY29lZmZpY2llbnRzKCnlvpfliLDmlbTpq5Tov7TmrbjmqKHlnovov7Tmrbjkv4LmlbhiZXRh77yM5Lulc2F2aW5nc+izh+aWmeeCuuS+izoKI3NyPTI4LjU2NjA5KygtMC40NjExOSkqcG9wMTUrKC0xLjY5MTUwKSpwb3A3NSsoLTAuMDAwMzQpKmRwaSsoMC40MDk2OSkqZGRwaQpyb3VuZChjb2VmZmljaWVudHMobWRsKSw1KQoj5L2/55So5oyH5Lukc3VtbWFyeSgp5Y+v5Lul5b2Z57i954mp5Lu2bWRs55qE5YWn5a6577yM5LqG6Kej5LiK6L+w5qih5Z6L5piv5ZCm5Y+v57Ch5YyW77yM5Y+K5o6i6KiO5qih5Z6L55qE6Kej6YeL6IO95Yqb44CCCnN1bW1hcnkobWRsKQoj5Lulc2F2aW5nc+izh+aWmeeCuuS+i++8jEYtc3RhdGlzdGljOiA1Ljc1NiBvbiA0IGFuZCA0NSBERu+8jHAtdmFsdWU6IDAuMDAwNzkwNOefpemBk+W8leWFpeino+mHi+iuiuaVuOaYr+aciee1seioiOaEj+e+qeOAggrpgJnlgIvmqKHlnottZGznmoTog73op6Pph4t555qE6K6K55Ww6YGUIDMzLjg1Je+8jOW+jOmdoumgiOeci+eci+mFjemBqeaumOW3ruWcluWRiOePvueahOWei+aFi+WPiuingOWvn+mboue+pOWAvOaYr+WQpuacieW9semfv+OAggrlho3nlLFiZXRh5L+C5pW455qE5b2Z57i96LOH5paZ55m854++cG9wNzXnmoRwLXZhbHVlPTAuMTI1NTMwID4gMC4wNeWSjHBvcDc155qEcC12YWx1ZT0wLjcxOTE3MyA+IDAuMDUg55qG5LiN5ouS57WVSG/lj6/ku6XogIPmha7miopwb3A3NeaIlmRwaeaLv+aOieOAggoj5Lul5LiK55qE57WQ5p6c6KaB5oiQ56uL5b+F6aCI5piv77yMcmVzaWR1YWxz5b+F6aCI5piv5ru/6Laz5bi45oWL5YGH6Kit77yM5LiU6K6K55Ww5pW454K65bi45pW444CCCiPlhYjliKnnlKjlgY9G5qqi5a6a5Yik5pa35LiK6Z2i5YWp5YCL5LiN5ouS57WVSG/kuYvorormlbjmmK/lkKboqbLlv73nlaXjgIIKI+ioreS4gOWAi+W/veeVpeaOiXBvcDc144CBZHBp6YCZ5YCL6K6K5pW455qE5paw55qE5qih5Z6L54K6bWRsMeOAggptZGwxPC1sbShzcn5wb3AxNStkZHBpLHNhdmluZ3MpCiPlho3kvb/nlKhhbm92YSgp5oyH5Luk5bCNbWRs5ZKMbWRsMeWFqeWAi+aooeWei+WBmuiuiueVsOaVuOWIhuaekO+8jOaqouWumnBvcDc15ZKMZHBp5Y+v5Lul5ZCM5pmC5b+955WlKEhvIDogcG9wNzU9ZHBpPTAp44CCCmFub3ZhKG1kbDEsbWRsKQoj55Sx6K6K55Ww5pW45YiG5p6Q6LOH5paZ55m854++55qEcC12YWx1ZT0wLjE5PiAwLjA1IOaJgOS7peS4jeaLkue1lUhv77yM6ICD5oWu55SoTW9kZWwx44CCCiPmrpjlt65WU+mFjemBqeWAvOmXnOS/guWclgpwbG90KGZpdHRlZChtZGwxKSxyZXNpZHVhbHMobWRsMSkpCm1lYW4oZml0dGVkKG1kbDEpKQoj5L2/55SodmFyLnRlc3QoKeaMh+S7pOaqouWumuipsuaooeWei21kbDHnmoTorornlbDmlbjmmK/lkKbngrrluLjmlbjvvIzphY3pganlgLznr4TlnI3lj6/nlKjmqKHlnovkuYvphY3pganlgLzkuYvlubPlnYflgLzjgIIKdmFyLnRlc3QocmVzaWR1YWxzKG1kbDEpW2ZpdHRlZChtZGwxKT45LjY3MV0scmVzaWR1YWxzKG1kbDEpW2ZpdHRlZChtZGwxKTw5LjY3MV0pCiPkvb/nlKhxcW5vcm0oKeaMh+S7pOi3keWHulFRLW5vcm3lnJbjgIIKcXFub3JtKHJlc2lkdWFscyhtZGwxKSkKcXFsaW5lKHJlc2lkdWFscyhtZGwxKSkKI+S9v+eUqHNoYXBpcm8udGVzdCgp5oyH5Luk5qqi5a6a5q2k5qih5Z6LbWRsMeS5i+W4uOaFi+WBh+ioreaYr+WQpuaIkOeri+OAggpzaGFwaXJvLnRlc3QocmVzaWR1YWxzKG1kbDEpKQoj5bi45oWL5YGH6Kit5Lul5Y+K6K6K55Ww5pW45qqi5a6a55qG5oiQ56uL77yM5pyA5b6M5YaN5qqi5p+l5piv5ZCm5pyJ6Zui576k5YC85a2Y5Zyo77yMCiPovInlhaXlpZfku7ZjYXLvvIzkvb/nlKhvdXRsaWVyVGVzdCgp5oyH5Luk5om+6Zui576k5YC844CCCmxpYnJhcnkoY2FyKQpvdXRsaWVyVGVzdChtZGwxKQoj5paH5a2X55qE6Zui576k5YC85Y+v5Lul54WnQ+Wkp+eahOaMh+S7pOWIqumZpOS6hiDmmK/miJHkuI3lsI/lv4PmjIfku6TmnInmiZPpjK89ID0rCiPkvYbmnInllY/poYzmmK9vdXRsaWVyVGVzdCgp5oyH5Luk5Y+q5pyD6LeR5Ye65LiA5YCL6Zui576k5YC8IOiAjOS4lOS4jeS4gOWumuacg+aYr+W9semfv+m7niDmh4noqbLmmK/kuI3og73kuoLli5UKI+acieWAi+aDs+azleaYr+iqquWIqeeUqGNvb2sncyBkaXN0YW5jZSBwbG905om+5Ye65b2x6Z+/6bueIOWGjeWwh+W9semfv+m7nuW+nuizh+aWmeS4reWIqumZpCDlho3ph43mlrDlu7rmqKHlho3lgZrmqqLlrpog5LiN55+l5Yiw5piv5ZCm5Y+v6KGMPwpjdXRvZmYgPC0gNC8oKG5yb3coc2F2aW5ncyktbGVuZ3RoKG1kbDEkY29lZmZpY2llbnRzKS0yKSkgCnBsb3QobWRsMSwgd2hpY2g9NCwgY29vay5sZXZlbHM9Y3V0b2ZmKQ==