例：swiss數據集
library(datasets)data(swiss)head(swiss)head(swiss) Fertility Agriculture Examination Education Catholic Infant.MortalityCourtelary 80.2 17.0 15 12 9.96 22.2Delemont 83.1 45.1 6 9 84.84 22.2Franches-Mnt 92.5 39.7 5 5 93.40 20.2Moutier 85.8 36.5 12 7 33.77 20.3Neuveville 76.9 43.5 17 15 5.16 20.6Porrentruy 76.1 35.3 9 7 90.57 26.6
hist(swiss$Catholic)圖1.關於天主教徒比例的各省數量的分布
這個分布是雙峰的，因為大多數省份要麼是天主教徒佔絕大多數，要麼是新教徒佔絕大多數。
library(dplyr)swiss = mutate(swiss, CatholicBin = 1 * (Catholic > 50))head(swiss) Fertility Agriculture Examination Education Catholic Infant.Mortality CatholicBin1 80.2 17.0 15 12 9.96 22.2 02 83.1 45.1 6 9 84.84 22.2 13 92.5 39.7 5 5 93.40 20.2 14 85.8 36.5 12 7 33.77 20.3 05 76.9 43.5 17 15 5.16 20.6 06 76.1 35.3 9 7 90.57 26.6 1使用dplyr包創建一個天主教的二值變量CatholicBin，如果該省多數是天主教徒，則為1；如果多數是新教徒，則為0。
g <- ggplot(swiss, aes(x = Agriculture, y = Fertility, color = factor(CatholicBin)))g <- g + geom_point(size = 6, color = "black") + geom_point(size = 4)g <- g + xlab("% in Agriculture") + ylab("Fertility")g圖2.不同宗教信仰下變量農業與生育率關係的散點圖
fit <- lm(Fertility ~ Agriculture, data = swiss)g1 <- gg1 <- g1 + geom_abline(intercept = coef(fit)[1], slope = coef(fit)[2], size = 2)g1summary(fit)$coef Estimate Std. Error t value Pr(>|t|)(Intercept) 60.3044 4.25126 14.185 3.216e-18Agriculture 0.1942 0.07671 2.532 1.492e-02圖3.不考慮宗教繪製農業與生育率的回歸線
不考慮宗教信仰，擬合生育率與從事農業比例之間的線性回歸模型，截距約為60，斜率約為0.19。
fit <- lm(Fertility ~ Agriculture + factor(CatholicBin), data = swiss)summary(fit)$coef Estimate Std. Error t value Pr(>|t|)(Intercept) 60.8322 4.1059 14.816 1.032e-18Agriculture 0.1242 0.0811 1.531 1.329e-01factor(CatholicBin)1 7.8843 3.7484 2.103 4.118e-02
g1 <- gg1 <- g1 + geom_abline(intercept = coef(fit)[1], slope = coef(fit)[2], size = 2)g1 <- g1 + geom_abline(intercept = coef(fit)[1] + coef(fit)[3], slope = coef(fit)[2], size = 2)g1圖4.宗教作為因子變量繪製農業與生育率的回歸線
將變量CatholicBin作為因子變量，且這個變量與從事農業比例沒有交互作用，擬合線性回歸模型，可以得到兩條斜率相同，截距不同的回歸線。
factor()將CatholicBin變為因子變量，否則R會將其視為連續回歸變量。係數(Intercept)與宗教信仰無關，係數factor(CatholicBin)1是天主教徒省份的截距變化。
fit <- lm(Fertility ~ Agriculture * factor(CatholicBin), data = swiss)summary(fit)$coef Estimate Std. Error t value Pr(>|t|)(Intercept) 62.04993 4.78916 12.9563 1.919e-16Agriculture 0.09612 0.09881 0.9727 3.361e-01factor(CatholicBin)1 2.85770 10.62644 0.2689 7.893e-01Agriculture:factor(CatholicBin)1 0.08914 0.17611 0.5061 6.153e-01
g1 <- gg1 <- g1 + geom_abline(intercept = coef(fit)[1], slope = coef(fit)[2], size = 2)g1 <- g1 + geom_abline(intercept = coef(fit)[1] + coef(fit)[3], slope = coef(fit)[2] + coef(fit)[4], size = 2)g1圖5.考慮宗教與農業交互作用擬合線性回歸
・對於新教徒佔大多數的省份，回歸模型截距為第一個係數(Intercept) ，約為62，斜率為第二個係數Agriculture，約為0.096。
・對於天主教徒佔大多數的省份，回歸模型截距為第一個係數(Intercept) 加上第三個係數factor(CatholicBin)1 ，約為64.8，斜率為第二個係數Agriculture加上第四個係數Agriculture:factor(CatholicBin)1，約為0.185。
由於天主教徒作用的係數都為正，故截距較高的回歸線為天主教徒回歸線。
殘差和回歸診斷例：swiss數據集
R中提供了大量檢驗回歸分析中統計假設的方法。最常見的方法就是對lm()函數返回的對象使用plot()函數，可以生成評價模型擬合情況的四幅圖形。
data(swiss)par(mfrow=c(2,2)) fit<-lm(Fertility~.,data=swiss)plot(fit)圖6.生育率與其他變量之間關係的回歸診斷圖
➢「殘差圖與擬合圖」（Residuals vs Fitted，左上），若因變量與自變量線性相關，那麼殘差值與預測（擬合）值就沒有任何系統關聯。——線性
➢「正態Q-Q圖」（Normal Q-Q，右上）是在正態分布對應的值下，標準化殘差的概率圖。——正態性。當預測變量值固定時，因變量成正態分布，則殘差值也應該是一個均值為0的正態分布。若滿足正態假設，圖上的點應該落在呈45度角的直線上；否則就違反了正態性的假設。
➢「位置尺度圖」（Scale-Location Graph，左下），若滿足不變方差假設，水平線周圍的點應該隨機分布。——同方差性
➢「殘差與槓桿圖」（Residuals vs Leverage，右下）可以鑑別出離群點、高槓桿值點和強影響點。
圖7.異常觀測點
離群值可能具有不同程度的影響。離群值可以符合回歸關係（如右上角點）。槓桿的概念 僅僅是數據點距自變量軸中心有多遠的概念，因變量值不參與計算一個觀測點的槓桿值。影響力在於該點對模型參數的估計發揮作用的大小。
◎左上角點是離群點具有低槓桿值，低影響力，以不符合回歸關係的方式存在。
◎左下角點不是離群點，槓桿值低，影響力低。
◎右上角點是離群點，也是高槓桿值點，具有較高的槓桿作用，但服從其他觀測點的回歸關係，具有較低的實際影響力。
◎右下角是離群點，也是強影響點，具有很高的槓桿作用，在擬合中則會對模型的參數估計發揮作用。
離群點可能是真實過程也可能是虛假過程產生的。擬合模型時希望從模型中刪除虛假過程產生的離群點同時保留真實過程產生的離群點。influence.measures()函數可以提供影響力的衡量指標?influence.measures。➢rstandard：標準化殘差，殘差除以其標準差。➢rstudent：標準化殘差，即殘差除以其標準差，在計算標準差的時刪除了第i個數據點使標準化殘差服從t分布 ➢hatvalues：衡量槓桿值，可用於判斷數據的輸入錯誤➢dffits：當刪除第i個點擬合模型時，預測值的變化➢dfbetas：當刪除第i個點擬合模型時，各個係數的變化➢resid(fit) / (1 - hatvalues(fit))：fit為線性模型擬合時返回PRESS殘差，如交叉驗證中留一驗證(leave-one-out cross-validation, LOOCV）的殘差，即數據i點的值與刪除該點擬合的模型在該點預測值的差值。PRESS殘差在檢測異常值和確定影響方面不是特別有用，而是用於確定模型擬合度。例：回歸診斷示例
n<-100x<-c(10,rnorm(n))y<-c(10,c(rnorm(n))) plot(x,y,frame=FALSE,cex=2,pch=21,bg="lightblue",col="black") abline(lm(y~x))圖8.生成有離群點的隨機數據集
隨機生成成對的獨立隨機數據，添加了離群點(10,10)。
fit<-lm(y~x)head(dfbetas(fit)) (Intercept) x1 -0.160472495 6.36542458932 0.044100714 -0.00069552743 0.054595229 -0.02281624184 -0.014275902 -0.01122142085 -0.005473504 -0.00122432526 0.059918030 -0.0025620951
round(dfbetas(fit)[1:10,2],3) #round()函數保留小數點後3位 1 2 3 4 5 6 7 8 9 10 6.365 -0.001 -0.023 -0.011 -0.001 -0.003 -0.004 -0.021 0.029 -0.059
round(hatvalues(fit)[1:10],3) 1 2 3 4 5 6 7 8 9 10 0.503 0.010 0.012 0.015 0.010 0.010 0.010 0.021 0.010 0.017可以看到第1個數據點(10,10)的dfbetas和hatvalues比其餘點的大的多，由於該點的存在，數據被認為有很強的相關性，否則數據相關係數估計為0，影響較大。hatvalues取值必須在0-1之間。
n<-100x<-c(10,rnorm(n))y<-c(10,x[2:n]+rnorm(n,sd=0.2))plot(x,y,frame=FALSE,cex=2,pch=21,bg="lightblue",col="black") abline(lm(y[2:n]~x[2:n]))圖9.添加位於回歸線上的離群點
沿著存在回歸關係的數據的回歸線生成數據，並添加位於回歸線上的離群點(10,10)。
fit2<-lm(y~x)round(dfbetas(fit2)[1:10,2],3) 1 2 3 4 5 6 7 8 9 10 0.585 -0.005 0.176 -0.001 -0.010 -0.013 -0.063 0.113 0.080 -0.021 round(hatvalues(fit2)[1:10],3) 1 2 3 4 5 6 7 8 9 10 0.509 0.010 0.019 0.010 0.011 0.011 0.015 0.028 0.015 0.011可以看到離群點(10,10)的dfbetas仍然很大，但沒有之前情況那麼大，對擬合的影響沒有那麼大；hatvalues依舊比其他大很多。該離群點超出了X值的範圍，但滿足其他觀測點的線性回歸關係，具有較大的槓桿值，但影響力並不大。
例：Stefanski和合作者共同創建的一個數據集，用來說明殘差圖的作用
dat<-read.table('http://www4.stat.ncsu.edu/~stefanski/NSF_Supported/Hidden_Images/orly_owl_files/orly_owl_Lin_4p_5_flat.txt')pairs(dat)圖10.各變量關係的散點圖矩陣 圖中有許多重複繪製的點。
summary(lm(V1~.-1,data=dat))$coef Estimate Std. Error t value Pr(>|t|)V2 0.9856157 0.12798121 7.701253 1.989126e-14V3 0.9714707 0.12663829 7.671225 2.500259e-14V4 0.8606368 0.11958267 7.197003 8.301184e-13V5 0.9266981 0.08328434 11.126919 4.778110e-28以V1為因變量，擬合多元線性模型。
fit<-lm(V1~.-1,data=dat)plot(predict(fit),resid(fit),pch='.')圖11.繪製殘差圖 殘差圖可以放大線性模型擬合不佳的方面，指出擬合模型遺漏的趨勢。