R语言总结（四）

二十六：多元线性回归分析

 

 #####多元线性回归
#多元线性回归的函数是lm()
#lm(因变量~自变量1+自变量2+...+自变量n,data = 数据框)
res <- lm(mydata$工作满意度~mydata$环境满意度+
     mydata$薪资满意度+mydata$工作内容满意度+mydata$工作能力,
   data = mydata)
full_res <- summary(res)
tab <- round(full_res$coefficients,3)
#模型的诊断

从以上的分析结果可以得出：

1、本次的多元线性回归模型的拟合度较好，调整R方=0.523>0.5。拟合度较好，意味着自变量能够比较好地反映出因变量的变化情况。换句话说，就是因变量的变化主要来源于自变量的影响。

2、通过F检验，得出回归模型显著（F（4,58）=17.99，P<0.001），意味着纳入本次回归模型的4个自变量中至少有一个自变量可以显著影响因变量。那么，究竟有几个自变量可以显著影响因变量？并且是哪几个？因此，还需要继续对每一个自变量的回归系数进行检验。

3、通过回归系数的T检验可以得出：

• 环境满意度不能显著影响工作满意度（T=1.453，P=0.152>0.05）。
• 薪资满意度可以显著影响工作满意度（T=2.584，P=0.012<0.05）；更进一步，薪资满意度的回归系数为0.489>0,因此，薪资满意度可以显著正向影响工作满意度，二者之间具体的定量关系是：员工的薪资满意度每提升1，工作满意度则随之提升0.489。
• 工作内容满意度可以显著影响工作满意度（T=-2.728，P=0.008<0.05）；更进一步，工作内容满意度的回归系数为-0.528<0,因此，工作内容满意度可以显著负向影响工作满意度，二者之间具体的定量关系是：员工的工作内容满意度每提升1，工作满意度则随之降低0.528。
• 工作能力可以显著影响工作满意度（T=5.751，P<0.05）；更进一步，工作能力的回归系数为0.98>0,因此，工作能力可以显著正向影响工作满意度，二者之间具体的定量关系是：员工的工作能力每提升1，工作满意度则随之升高0.98。

最后，基于以上的分析结果可以得出以下的回归方程：

工作满意度=0.489*薪资满意度-0.528*工作内容满意度+0.98*工作能力

回归模型的诊断：

1. 线性回归模型的数据诊断
   1. 线性回归模型要求样本之间不存在序列相关（独立填写）
   #1.1 DW 序列相关
   #计算DW需要使用lmtest程序包
   library(lmtest)
   dwtest(res)
2. 

考察样本是否存在序列相关的统计量是DW，本次回归模型的DW值为1.13，DW只有2附近（1.9-2.1），才认为样本之间不存在序列相关。本次的运算结果意味着模型存在轻微序列相关。

1. 线性回归模型要求自变量之间不存在共线性（两个自变量反应同一个内容）

线性回归可以通过VIF值考察自变量之间的共线性。本次回归模型的4个自变量的VIF值全部小于5，意味着不存在多重共线性。

#1.2 VIF值
#VIF值的计算会用到car包中的vif()函数
#VIF值小于5则认为自变量之间不存在共线性
library(car)
vif <- round(vif(res),3)
vif <- c('NA', vif)
vif
tab <- data.frame(tab, vif)
write.csv(tab, 'tab.csv')

• 线性回归模型优良状况

线性回归模型的优良状况全部由残差来反应。预测值和实际值之间的差异，就是残差。

2.1 残差的数值不能太大

2.2 没有异方差（残差不能存在趋势）

2.3 残差要服从正态分布

本次的多元线性回归模型的残差数值不大，不存在异方差，且基本服从正态分布。具体情况见下图：

二十七：二元logistic回归分析

 ######二元logistic回归分析###
#1、导入数据
library(Hmisc)   #导入SPSS格式的数据需要用到的程序包
mydata <- spss.get('NO.31 配套数据.sav',
                   use.value.labels = TRUE)
str(mydata)
mydata$年龄 <- as.numeric(mydata$年龄)
#确保自变量的类型是准确的之后，开始进行二元logistic回归分析

#2、二元logistic回归分析
#使用的函数是glm()广义线性模型
#general linear model  函数glm()
#linear model   函数是lm()
#glm()的代码格式
#glm(方程, family = XXXXX, data = XXXX)
#family = XXXXX表示用哪一种模型 binomial二分类
model <- glm(mydata$糖尿病=='有'~mydata$BMI+mydata$性别+
               mydata$吸烟+mydata$饮酒+mydata$家族史+
               mydata$年龄+mydata$腰围, 
             family = binomial, 
             data = mydata)   #建立的二元logistic回归模型
View(model)

res <- summary(model)  #查看完整的运算结果

res$coefficients
tab <- round(res$coefficients,3)
model$coefficients    #仅仅只有回归系数
#summary(model)可以看到完整的回归系数和回归系数的显著性检验
#model仅仅只可以看到回归系数，看不到回归系数的显著性检验
#计算OR值只需要回归系数。
OR <- round(exp(model$coefficients),3)
OR

#计算OR的95%置信区间
model$coefficients
confint(model)
ORCI <- round(exp(confint(model)),3)
#OR（下~上）
OR的置信区间 <-paste0(OR,'(',ORCI[,1],'~', ORCI[,2],')')
fianl <- data.frame(tab,OR的置信区间)
#导出结果
write.csv(fianl,'fianl.csv')

从以上的分析结果可以清晰得出：

1. 性别、吸烟、饮酒、腰围、BMI这几个因素的回归系数均没有通过显著性水平为0.05的检验，P值全部大于0.05，意味着这几个因素不是糖尿病患病的影响因素。
2. 家族史是糖尿病患病的影响因素，P<0.05；有家族史患糖尿病的可能性是无家族史的6.361倍（OR值，计算exp(回归系数）），置信区间是3.065~14.026

年龄是糖尿病患病的影响因素，P=0.034<0.05；年龄大1岁出现糖尿病的可能性是年龄小1岁的1.052倍（年龄每增加1岁，患病的可能性就增加5.2%）。

接下来，对于以上的回归模型，还需要诊断模型的优良度。

#二元logistic回归模型的显著性检验
#也是用到lmtest
library(lmtest)
modelnull <- glm(mydata$糖尿病=='有'~1, 
                 family = binomial,
                 data = mydata)
#空模型和全模型的差异对比函数是lrtest
lrtest(model, modelnull)

    回归模型的显著性（自变量是否影响因变量）

   通过对比空模型和全模型之间的差异，得出P<0.05,意味着含有自变量的模型与不含有自变量的模型之间存在显著差异，即自变量中至少有一个可以显著影响因变量，即回归模型吧显著。

#伪R方
#伪R方的计算需要用到pscl程序包中的pR2()函数
library(pscl)
pR2(model) 

    回归模型的拟合度

   本回归模型的伪R2为0.288，拟合一般，但是对于二元logistic回归模型而言，伪R2对拟合度的反映不太准确。因此还需要继续考察模型预测的准确率。

#预测准确率
prob <- predict(model, type = 'response')
#割点是0.5，即>=0.5认为有病，<0.5无病。
type <- ifelse(prob >= 0.5, '有1', '无1')

mydata <- data.frame(mydata, type)
aaa <- table(mydata$糖尿病, mydata$type.1)
aaa <- data.frame(aaa)
(83+61)/(83+17+39+61)
#得出模型的预测准确率为72%
write.csv(aaa,'aaa.csv')

   预测准确率

  模型预测的准确率为72%，即模型的预测能力良好。见下表：