R语言——Logistic回归

1.数据

The birthwt data frame has 189 rows and 10 columns. The data were collected at Baystate Medical Center, Springfield, Mass during 1986.

variable	label
low	indicator of birth weight less than 2.5 kg.
age	mother's age in years.
lwt	mother's weight in pounds at last menstrual period.
race	mother's race (1 = white, 2 = black, 3 = other).
smoke	smoking status during pregnancy.
ptl	number of previous premature labours.
ht	history of hypertension.
ui	presence of uterine irritability.
ftv	number of physician visits during the first trimester.
bwt	birth weight in grams.

加载MASS包中的birthwt数据

#加载MASS包中的birthwt数据
data(birthwt,package = 'MASS')
#查看数据的结构
str(birthwt)


根据数据的结构进行数据处理，原始变量的race、smoke、ht、ui都是数值型变量，在分析前需要转化成因子。另外ptl表示先前早产次数、变量ftv表示怀孕前三个月探访医生次数，可以先查看两个变量的频数分布。

#将数据框添加到搜索路径
attach(birthwt)
#加载epiDisplay包
library(epiDisplay)
#查看ptl和ftv的频数分布
tab1(ptl) 
tab1(ftv) 

结果发现，84.1%的孕妇没有早产史，早产次数超过一次的孕妇也很少，所以有必要将ptl转化为二分类的因子（0：没有早产史，1 ：有早产史）。同理我们也可以将ftv转化为一个三分类的因子（0：未探访医生，1探访医生1次，2 ：探访医生大于1次）。具体代码操作如下,并将整理好的数据命名为birthweight。

#加载dplyr包
library(dplyr)
birthweight <- birthwt %>%
  mutate(race=factor(race,labels = c('white','black','other')),
         smoke=factor(smoke,labels = c('no','yes')),
         ptl=ifelse(ptl>0,'1 ',ptl),
         ptl=factor(ptl,labels = c('no','yes')),
         ht=factor(ht,labels = c('no','yes')),
         ui=factor(ui,labels = c('no','yes')),
         ftv=ifelse(ftv>1,'2 ',ftv),
         ftv=factor(ftv,labels = c('no','1次','>1次')))
#查看数据结构
str(birthweight)


2.模型的建立

实际上变量low变量是由bwt生成的，我们将其作为结果变量（取0和1），建立二元logistic回归模型探索新生儿低体重的影响因素。在R语言中，函数glm()用于拟合包括Logistic回归在内的广义线性模型。

#建立二元Logistic回归模型
glm1 <-glm(low~age lwt race smoke ptl ht ui ftv,
           family = binomial,data = birthweight) 
#使用summary()函数提取模型的汇总结果
summary(glm1)

上面输出第一部分展示了模型调用的公式。第二部分给出了残差的分布。第三部分给出了回归系数的估计、标准误和显著性检验结果。我们通常对第三部分最感兴趣。接着给出的是两个残差平方和，Null deviance表示只包含常数项的模型的残差平方和，Residual deviance表示现模型的残差平方和。我们感兴趣的是两个残差平方和的差值，这里是234.67-195.48=39.19，这个值可以用来检验模型的显著性。残差平方和渐进服从卡方分布，我们可以用pchisq()函数得到。结果的最后部分还包括了模型的AIC值和模型的拟合过程中迭代次数。AIC值将模型中参数的个数考虑进来，度量模型拟合的好坏，可以用于比较模型。迭代次数是算法指标，一般不关注，但是这个值如果太大，比如25，表示对现有数据而言模型太复杂。

3.自变量的筛选

与多重线性回归模型类似，当自变量的个数较多时，为了使建立的Logistic回归模型比较稳定和便于解释，应尽量将回归效果显著的自变量筛选入模型，将效果不显著的变量剔除在外。我们可以使用drop1()函数得到剔除一个自变量后的模型的AIC值。

#drop1()函数得到剔除一个自变量后的模型的AIC值。
drop1(glm1)


从上面数据可以看出，第一步需要在模型中剔除变量ftv，因为剔除该变量后，模型的AIC值最小。第二步可以对剔除ftv之后的模型使用drop1()函数，以此类推。
可以使用step()函数一次性自动完成这个过程。

#逐步回归,direction可以定义是向前向后、前进法、后退法，trace=T表示显示筛选详细过程
glm2 <- step(glm1,direction='both',trace=T)
summary(glm2)

4.模型的比较

#将上述两个模型进行比较
anova(glm1,glm2,test = 'Chisq')

anova()函数中的test参数可以设为LRT,得到的结果与上面相同。结果显示P=0.4981>0.05,表明glm1和glm2两个模型拟合的一样好。也就是说年龄和孕早期探访医生的次数不会显著提高模型的预测精度，我们更倾向于使用更简单的模型。函数AIC也可以比较模型。

#使用AIC函数比较模型
AIC(glm1,glm2)

5.回归系数的解释

#查看回归系数
coef(glm2)


在logistic回归中，模型拟合的是响应变量y=1时对数优势比（log(odds)）。因此，回归系数的含义是当其他变量保持不变时，一个单位预测变量的变化可引起的响应变量对数优势比的变化，为了便于解释，可以将结果指数化转换成优势比。

exp(coef(glm2))

结果表明，lwt（母亲怀孕期体重）每增加1磅，新生儿体重的优势将乘以0.9842（保持其他变量不变）。对于连续变量，一个单位的变化可能并不太好解释。如果lwt增加10磅，优势将乘以0.9842^10，即约为0.817。
黑人分娩低体重儿的优势约为白人孕妇的3.67倍，其他种族分娩低体重儿的优势约为白人孕妇的2.35倍。类似，吸烟、有早产史、高血压、子宫应激症的孕妇分娩低体重儿的优势都将增加。

#回归系数的置信区间
confint(glm2)
#把上面系数的置信区间指数化
exp(confint(glm2))


6.预测

对于大多数人来说概率比优势更好理解。给定一组自变量的值，我们可以使用建立的模型得到预测值。例如，建立如下虚拟的孕妇样本。

newdata <- data.frame(lwt=120,race='black',smoke='yes',ptl='no',ht='yes',ui='no')
#使用函数predict()得到模型的预测值
predict(glm2,newdata = newdata,type = 'response')

7.模型的检查

在评价模型时，Hosmer-Lemeshow检验可以用于判断观测值和预测值之间差异的显著性。ResourceSelection包中的函数hoslem.test()可以执行该检验。

#Hosmer-Lemeshow检验
library(ResourceSelection)
hoslem.test(low,fitted(glm2))


Hosmer-Lemeshow原假设是观测值和预测值之间的差异无统计学意义。上面的结果表明差异没有统计学意义，p=0.3593,因此不能拒绝关于模型很好地拟合了模型的假设。

8.模型结果的汇总输出

使用epiDisplay包中的函数logistic.display()函数可以用于汇总Logistic回归模型的主要结果。

#模型结果汇总
logistic.display(glm2)

上述输出结果包含了各个自变量调整前和调整后的优势比和95%置信区间，以及针对回归系数的Wald检验的P值和针对自变量的似然比检验的P值，这是大多数论文报告的Logistic回归的结果。我们可以将输出结果导出到一个.csv格式的文件中。

write.csv(logistic.display(glm2)$table,file='D:/model.csv')


参考：《R语言医学数据分析实战》