各类统计方法R语言实现（六）
2020-05-06 本文已影响0人
拾光_2020
今天是各类统计方法R语言实现的第六期，我们主要介绍多元线性回归、回归诊断。
多元线性回归
多元线性回归指的是用多个自变量预测一个因变量，且自变量与因变量之间为线性关系，在分析过程中要考虑交互项的问题。
数据输入
使用的数据自变量为第一次推文各类统计方法R语言实现（一）中所用的mtcar数据集，此数据集中mpg表示每加仑油行驶英里数，hp马力，wt车重。此次分析的因变量是mpg，自变量是hp与wt。
data(mtcars)
head(mtcars)
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
#计算描述统计分析
mtcars$cyl<-as.factor(mtcars$cyl)
mtcars$vs<-as.factor(mtcars$vs)
mtcars$am<-as.factor(mtcars$am)
mtcars$gear<-as.factor(mtcars$gear)
mtcars$carb<-as.factor(mtcars$carb)
#计算描述统计分析
summary(mtcars)
## mpg cyl disp hp drat
## Min. :10.40 4:11 Min. : 71.1 Min. : 52.0 Min. :2.760
## 1st Qu.:15.43 6: 7 1st Qu.:120.8 1st Qu.: 96.5 1st Qu.:3.080
## Median :19.20 8:14 Median :196.3 Median :123.0 Median :3.695
## Mean :20.09 Mean :230.7 Mean :146.7 Mean :3.597
## 3rd Qu.:22.80 3rd Qu.:326.0 3rd Qu.:180.0 3rd Qu.:3.920
## Max. :33.90 Max. :472.0 Max. :335.0 Max. :4.930
## wt qsec vs am gear carb
## Min. :1.513 Min. :14.50 0:18 0:19 3:15 1: 7
## 1st Qu.:2.581 1st Qu.:16.89 1:14 1:13 4:12 2:10
## Median :3.325 Median :17.71 5: 5 3: 3
## Mean :3.217 Mean :17.85 4:10
## 3rd Qu.:3.610 3rd Qu.:18.90 6: 1
## Max. :5.424 Max. :22.90 8: 1
代码展示，首先不考虑交互项
#模型拟合
fit<-lm(mpg~hp+wt,data=mtcars)
##展示模型
summary(fit)
##
## Call:
## lm(formula = mpg ~ hp + wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.941 -1.600 -0.182 1.050 5.854
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.22727 1.59879 23.285 < 2e-16 ***
## hp -0.03177 0.00903 -3.519 0.00145 **
## wt -3.87783 0.63273 -6.129 1.12e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.593 on 29 degrees of freedom
## Multiple R-squared: 0.8268, Adjusted R-squared: 0.8148
## F-statistic: 69.21 on 2 and 29 DF, p-value: 9.109e-12
##利用模型计算预测值
fitted(fit)
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## 23.572329 22.583483 25.275819 21.265020
## Hornet Sportabout Valiant Duster 360 Merc 240D
## 18.327267 20.473816 15.599042 22.887067
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 21.993673 19.979460 19.979460 15.725369
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## 17.043831 16.849939 10.355205 9.362733
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 9.192487 26.599028 29.312380 28.046209
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## 24.586441 18.811364 19.140979 14.552028
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 16.756745 27.626653 26.037374 27.769769
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## 16.546489 20.925413 12.739477 22.983649
#计算每一个预测值与实际值之间的差值，即残差
residuals(fit)
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## -2.57232940 -1.58348256 -2.47581872 0.13497989
## Hornet Sportabout Valiant Duster 360 Merc 240D
## 0.37273336 -2.37381631 -1.29904236 1.51293266
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 0.80632669 -0.77945988 -2.17945988 0.67463146
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## 0.25616901 -1.64993945 0.04479541 1.03726743
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 5.50751301 5.80097202 1.08761978 5.85379085
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## -3.08644148 -3.31136386 -3.94097947 -1.25202805
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 2.44325481 -0.32665313 -0.03737415 2.63023081
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## -0.74648866 -1.22541324 2.26052287 -1.58364943
结果解读：
回归系数Estimate 的含义是：一个预测变量每增加一个单位，其他预测变量保持不变时，因变量要增加的数量。本例中hp与wt的p值均小于0.05，表明两个指标均可纳入模型，保持马力不变，wt每增加一个单位（1000磅），mpg减少3.87783英里/加仑(US)
所有预测变量解释了82.68%的方差。
代码展示，考虑交互项
#模型拟合
fit<-lm(mpg~hp+wt+hp:wt,data=mtcars)
##展示模型
summary(fit)
##
## Call:
## lm(formula = mpg ~ hp + wt + hp:wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0632 -1.6491 -0.7362 1.4211 4.5513
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.80842 3.60516 13.816 5.01e-14 ***
## hp -0.12010 0.02470 -4.863 4.04e-05 ***
## wt -8.21662 1.26971 -6.471 5.20e-07 ***
## hp:wt 0.02785 0.00742 3.753 0.000811 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.153 on 28 degrees of freedom
## Multiple R-squared: 0.8848, Adjusted R-squared: 0.8724
## F-statistic: 71.66 on 3 and 28 DF, p-value: 2.981e-13
##利用模型计算预测值
fitted(fit)
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## 23.09547 21.78138 25.58488 20.02924
## Hornet Sportabout Valiant Duster 360 Merc 240D
## 17.28996 18.88542 15.40745 21.65887
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 20.84992 18.55379 18.55379 15.14994
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## 16.23929 16.07909 12.02179 11.89490
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 12.50221 27.84866 32.63195 30.24587
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## 24.56317 17.57441 17.91776 15.03111
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 15.93596 29.53900 26.71871 28.56630
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## 15.36033 19.52990 13.54587 22.31363
#计算每一个预测值与实际值之间的差值，即残差
residuals(fit)
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## -2.0954741 -0.7813755 -2.7848771 1.3707560
## Hornet Sportabout Valiant Duster 360 Merc 240D
## 1.4100448 -0.7854161 -1.1074453 2.7411309
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 1.9500834 0.6462128 -0.7537872 1.2500604
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## 1.0607148 -0.8790873 -1.6217868 -1.4949003
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 2.1977932 4.5513369 -2.2319540 3.6541302
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## -3.0631732 -2.0744146 -2.7177637 -1.7311118
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 3.2640401 -2.2390044 -0.7187056 1.8336953
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## 0.4396692 0.1701019 1.4541328 -0.9136259
可以看到交互项hp:wt也是显著的表明因变量与其中一个自变量的关系依赖于另一个预测变量的水平，在本例中，每加仑汽油行驶英里数与汽车马力的关系依赖于车重不同而不同。
图形展示交互项的结果
我们也可以通过图形来展示交互项的结果
library(effects)
## Warning: package 'effects' was built under R version 3.6.3
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
plot(effect("hp:wt", fit, xlevels=list(wt=c(2.2,3.2,4.2))), multiline=TRUE)
image.png
结果解释： 1、随着车重的增加，马力与每加仑汽油行驶英里数的关系会减弱。 2、当wt=4.2时，直线几乎水平，表明此时随着hp的增加，mpg不会发生改变。 通过该可视化结果，我们更易理解交互项的意义。
回归诊断
当构建好模型之后，如何判断模型是否合适呢？
主要是看模型在多大程度上满足最小二乘（OLS）模型统计假设，这在上次推文中介绍过，主要有以下四项：（括号中为假设相应的统计学指标）
①正态性：对于固定自变量，因变量成正态分布。（表明残差服从正态分布）
②独立性：因变量取值相互独立。（依赖于数据来源，若数据来自随机取样，则一般为独立；若来自某一家庭之类的互项有关联的团体，则不一定独立。）
③线性：自变量和因变量之间为线性相关。（表明残差与预测值没有关联）
④同方差性：因变量方差不随自变量的水平不同而变化。（位置尺度图（Scale-Location Graph）水平线周围的点随机分布）
对于简单线性回归的诊断
#简单线性模型
fit<-lm(weight~height,data=women)
par(mfrow=c(2,2))
plot(fit)
image.png
结果解释：
①左上角代表的是线性，横轴为预测值，纵轴为残差，可以发现二者具有曲线关系，表明可能需要尝试多项式回归。
②右上角为Q-Q图，反映线性假设，与之前的正态性检验是类似的，此处纵坐标为标准化残差。在回归分析中，测定值与按回归方程预测的值之差称为残差，以δ表示。（δ-残差的均值）/残差的标准差，称为标准化残差，以δ*表示，即z-标准化之后的残差。这张图看起来勉勉强强，如果需要确认是否服从正态分布，还可做正态分布检验（见下面代码）。
③左下角反映同方差性，水平线周围的点随机分布，该图大致满足这一条件。
④右下角为残差与杠杆图，可用于鉴别离群点、高杠杆点和强影响点，但是该图可读性不强，下一次推文会介绍更好的呈现方法。
以下为该三类点的定义：
a.离群点：拟合回归模型对其预测效果不佳（即残差的绝对值较大）。
b.有高杠杆值的变量表明它是一个异常的自变量组合。
c.强影响点表明他对模型参数的估计产生的影响过大。
#计算残差
residuals(fit)
## 1 2 3 4 5 6
## 2.41666667 0.96666667 0.51666667 0.06666667 -0.38333333 -0.83333333
## 7 8 9 10 11 12
## -1.28333333 -1.73333333 -1.18333333 -1.63333333 -1.08333333 -0.53333333
## 13 14 15
## 0.01666667 1.56666667 3.11666667
shapiro.test(residuals(fit))
##
## Shapiro-Wilk normality test
##
## data: residuals(fit)
## W = 0.91909, p-value = 0.1866
#z-score标准化
(residuals(fit) - mean(residuals(fit)))/sd(residuals(fit))
## 1 2 3 4 5 6
## 1.64451468 0.65780587 0.35158590 0.04536592 -0.26085405 -0.56707403
## 7 8 9 10 11 12
## -0.87329400 -1.17951397 -0.80524512 -1.11146509 -0.73719623 -0.36292738
## 13 14 15
## 0.01134148 1.06609917 2.12085686
shapiro.test((residuals(fit) - mean(residuals(fit)))/sd(residuals(fit)))
##
## Shapiro-Wilk normality test
##
## data: (residuals(fit) - mean(residuals(fit)))/sd(residuals(fit))
## W = 0.91909, p-value = 0.1866
#z-score标准化简化代码
scale(residuals(fit))
## [,1]
## 1 1.64451468
## 2 0.65780587
## 3 0.35158590
## 4 0.04536592
## 5 -0.26085405
## 6 -0.56707403
## 7 -0.87329400
## 8 -1.17951397
## 9 -0.80524512
## 10 -1.11146509
## 11 -0.73719623
## 12 -0.36292738
## 13 0.01134148
## 14 1.06609917
## 15 2.12085686
## attr(,"scaled:center")
## [1] 4.810966e-17
## attr(,"scaled:scale")
## [1] 1.469532
shapiro.test(scale(residuals(fit)))
##
## Shapiro-Wilk normality test
##
## data: scale(residuals(fit))
## W = 0.91909, p-value = 0.1866
可以看到，上面展示了两种z-score标准化方法，z-score标准化前后Shapiro-Wilk正态性检验结果不变。 该残差数据p>0.05，服从正态分布。
对于多项式回归的诊断
#多项式回归模型
fit2<-lm(weight ~ height + I(height^2),data = women)
summary(fit2)
##
## Call:
## lm(formula = weight ~ height + I(height^2), data = women)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50941 -0.29611 -0.00941 0.28615 0.59706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 261.87818 25.19677 10.393 2.36e-07 ***
## height -7.34832 0.77769 -9.449 6.58e-07 ***
## I(height^2) 0.08306 0.00598 13.891 9.32e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3841 on 12 degrees of freedom
## Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994
## F-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(fit2)
image.png
可以看出基本满足线性、残差正态性、同方差性。点13影响了线性假设，点15的cook距离较大，类似于强影响点。
newfit <- lm(weight~ height + I(height^2), data=women[-c(13,15),])
summary(newfit)
##
## Call:
## lm(formula = weight ~ height + I(height^2), data = women[-c(13,
## 15), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.43326 -0.24742 0.01186 0.23997 0.39241
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 235.286916 25.546802 9.210 3.36e-06 ***
## height -6.509449 0.796679 -8.171 9.78e-06 ***
## I(height^2) 0.076471 0.006194 12.347 2.23e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3198 on 10 degrees of freedom
## Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994
## F-statistic: 1.06e+04 on 2 and 10 DF, p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(newfit)
image.png
注意此处删除数据需要非常谨慎，因为本应是用模型匹配数据，而非数据匹配模型。
线性模型假设的综合验证
使用gvlma包中的gvlma()函数。
fit2<-lm(weight ~ height+ I(height^2),data = women)
summary(fit2)
##
## Call:
## lm(formula = weight ~ height + I(height^2), data = women)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50941 -0.29611 -0.00941 0.28615 0.59706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 261.87818 25.19677 10.393 2.36e-07 ***
## height -7.34832 0.77769 -9.449 6.58e-07 ***
## I(height^2) 0.08306 0.00598 13.891 9.32e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3841 on 12 degrees of freedom
## Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994
## F-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16
library(gvlma)
gvmodel<-gvlma(fit2)
summary(gvmodel)
##
## Call:
## lm(formula = weight ~ height + I(height^2), data = women)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50941 -0.29611 -0.00941 0.28615 0.59706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 261.87818 25.19677 10.393 2.36e-07 ***
## height -7.34832 0.77769 -9.449 6.58e-07 ***
## I(height^2) 0.08306 0.00598 13.891 9.32e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3841 on 12 degrees of freedom
## Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994
## F-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = fit2)
##
## Value p-value Decision
## Global Stat 10.858254 0.028204 Assumptions NOT satisfied!
## Skewness 0.005991 0.938306 Assumptions acceptable.
## Kurtosis 1.027270 0.310801 Assumptions acceptable.
## Link Function 9.128078 0.002517 Assumptions NOT satisfied!
## Heteroscedasticity 0.696916 0.403822 Assumptions acceptable.
可以看出Global Stat不满足条件，表明不符合某项或某些假设条件，可以用之前的方法来确定不符合哪项或哪些假设条件。
好了，今天的R语言实现统计方法系列推文暂时告一段落，我们下次再见吧！ 小伙伴们如果有什么统计上的问题，或者如果想要学习什么方面的生物信息内容，可以在微信群或者知识星球提问，没准哪天的推文就是专门解答你的问题哦！
