R语言与统计-5：Logistic回归
R语言与统计-1：t检验与秩和检验
R语言与统计-2：方差分析
R语言与统计-3：卡方检验
R语言与统计-4：线性回归分析与模型诊断
Logistic回归和线性回归的本质区别是Logistic回归的因变量y值是分类变量（二分类或者多分类）。最常见的例子是患者的结局是生存还是死亡（二分类）。
包括两分类logistic回归、无序多分类logistic回归、有序多分类logistic回归、条件logistic回归。
1. 两分类logistic回归
在R语言中，使用glm()函数来进行两分类Logistic回归。
例1
library(HSAUR2)
data('plasma')
head(plasma)
# fibrinogen globulin ESR
# 1 2.52 38 ESR < 20
# 2 2.56 31 ESR < 20
# 3 2.19 33 ESR < 20
# 4 2.18 31 ESR < 20
# 5 3.41 37 ESR < 20
# 6 2.46 36 ESR < 20
以分类变量ESR作为反应变量，以前面两个连续变量进行自变量来建模
fit1 <- glm(ESR~fibrinogen+globulin,data = plasma,family = binomial())
# family = binomial()表示是二分类
summary(fit1)
# Call:
# glm(formula = ESR ~ fibrinogen + globulin, family = binomial(),
# data = plasma)
# Deviance Residuals:
# Min 1Q Median 3Q Max
# -0.9683 -0.6122 -0.3458 -0.2116 2.2636
# Coefficients:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -12.7921 5.7963 -2.207 0.0273 *
# fibrinogen 1.9104 0.9710 1.967 0.0491 *
# globulin 0.1558 0.1195 1.303 0.1925
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# (Dispersion parameter for binomial family taken to be 1)
# Null deviance: 30.885 on 31 degrees of freedom
# Residual deviance: 22.971 on 29 degrees of freedom
# AIC: 28.971
# Number of Fisher Scoring iterations: 5
可以看到fibrinogen是有统计学意义的，globulin是没有统计学意义的。
计算fibrinogen的OR值：也就是这个变量每增加一个单位，ESR由一个水平变为另一个水平的概率增加多少倍
exp(coef(fit1)['fibrinogen'])
# fibrinogen
# 6.755579
exp(1.9104) #exp()进行系数转换
# [1] 6.755791
计算置信区间
置信区间（Confidence interval）是指由样本统计量所构造的总体参数的估计区间。在统计学中，一个概率样本的置信区间是对这个样本的某个总体参数的区间估计。置信区间展现的是这个参数的真实值有一定概率落在测量结果的周围的程度，其给出的是被测量参数的测量值的可信程度，即前面所要求的“一个概率”。
exp(confint(fit1,parm='fibrinogen'))
# Waiting for profiling to be done...
# 2.5 % 97.5 %
# 1.404131 73.000836
例2
library(MPV)
#家庭收入与是否拥有独立住房(是1/否0)之间的关系
knitr::kable(head(p13.2))
# | x| y|
# |-----:|--:|
# | 38000| 0|
# | 51200| 1|
# | 39600| 0|
# | 43400| 1|
# | 47700| 0|
# | 53000| 0|
log_fit <- glm(y~x,family = binomial(),data = p13.2) #binomial表示二分类
# Call:
# glm(formula = y ~ x, family = binomial, data = p13.2)
# Deviance Residuals:
# Min 1Q Median 3Q Max
# -2.0232 -0.8766 0.5072 0.7980 1.6046
# Coefficients:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -8.7395139 4.4394326 -1.969 0.0490 *
# x 0.0002009 0.0001006 1.998 0.0458 *
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# (Dispersion parameter for binomial family taken to be 1)
# Null deviance: 27.526 on 19 degrees of freedom
# Residual deviance: 22.435 on 18 degrees of freedom
# AIC: 26.435
# Number of Fisher Scoring iterations: 4
0.0002009是beta1的估计值
在线性回归里面，0.0002009代表x平均变化一个单位，y相应的变化多少个单位。但是在logistic回归里面，我们对原始p进行了logistic转换，所以这里的0.0002009是一个相对的比值比的结果，要对其进行对数转换才能得到OR值（比值比）。
exp(coef(log_fit)[2])
# x
# 1.000201
OR值是1.000201，也就是在其他变量不变的情况下，x值每变化一个单位，y会相应的增加1.000201倍
exp(confint(log_fit)[2,])
# Waiting for profiling to be done...
# 2.5 % 97.5 %
# 1.000024 1.000437
OR值的95%置信区间（不包括1，有统计学意义，等价于p<0.05）
过度离势
抽样于二项分布的数据的期望方差是δ2=nπ(1-π)，n为观测数，π为属于Y=1组的概率。所谓过度离势，即观测到的响应变量的方差大于期望的二项分布的方差。过度离势会导致奇异的标准误检验和不精确的显著性检验。
检测过度离势的一种方法是比较二项分布模型的残差偏差与残差自由度，如下：
当其值远大于1时，认为存在过度离势。同样也可以对其进行检验，pchisq(summary(fit.od)$dispersion * fit$df.residual, fit$df.residual, lower = F)
当出现过度离势时，是需要进行校正的。可使用 family=quasibinomial()对family = binomial()的部分进行替换。
log_fit <- glm(y~x,family = quasibinomial(),data = p13.2)
summary(log_fit)
# Call:
# glm(formula = y ~ x, family = quasibinomial(), data = p13.2)
# Deviance Residuals:
# Min 1Q Median 3Q Max
# -2.0232 -0.8766 0.5072 0.7980 1.6046
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -8.7395139 4.7522603 -1.839 0.0825 .
# x 0.0002009 0.0001077 1.866 0.0784 .
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# (Dispersion parameter for quasibinomial family taken to be 1.145897)
# Null deviance: 27.526 on 19 degrees of freedom
# Residual deviance: 22.435 on 18 degrees of freedom
# AIC: NA
# Number of Fisher Scoring iterations: 4
2. 无序多分类logistic回归
对于无序多分类logistic回归，可将其看作非条件二分类logistic回归的扩展。分析过程是将多分类因变量拆分成多个二分类变量，拟合n-1个(n为因变量的分类数)二分类logistic回归，变量的回归系数β、效应量OR值的计算等在n-1个模型中分别进行。假设某一疾病有A、B、C三个表型，将其分别赋值为1、2、3 (哑变量)，若以A表型(1)为参照，那么该模型的表示如下：
上面一个公式是2和1去比，下面一个公式是3和1去比。
无序多分类logistic回归需要满足3个条件：
条件1：因变量唯一，且为无序多分类变量。
条件2：存在一个或多个自变量，可为定性与定量变量。
条件3：一般要求例数较少因变量类的观察例数为自变量个数的10~15倍(EPV原则)且经验上每组的人数最好多于30例，自变量的参照水平组不应少于30或50例。
做无序多分类logistic回归可以使用nnet包中的multinom函数。
例：演示数据集是高中生进入高中之后会选择的编程难度
library(foreign) #读取网页数据
library(nnet)
ml <- read.dta('https://stats.idre.ucla.edu/stat/data/hsbdemo.dta')
head(ml)
# id female ses schtyp prog read write math science socst honors awards cid
# 1 45 female low public vocation 34 35 41 29 26 not enrolled 0 1
# 2 108 male middle public general 34 33 41 36 36 not enrolled 0 1
# 3 15 male high public vocation 39 39 44 26 42 not enrolled 0 1
# 4 67 male low public vocation 37 37 42 33 32 not enrolled 0 1
# 5 153 male middle public vocation 39 31 40 39 51 not enrolled 0 1
# 6 51 female high public general 42 36 42 31 39 not enrolled 0 1
因变量prog包括：vocation业余的，general普通的，academic学术的(3分类)
# 将academic指定为reference
ml$prog2 <- relevel(ml$prog,ref = "academic")
# 自变量是ses（社会经济水平）和write
test <- multinom(prog2~ses+write,data = ml)
## weights: 15 (8 variable)
# initial value 219.722458
# iter 10 value 179.982880
# final value 179.981726
# converged
summary(test)
先看Coefficients：因为academic被作为参照，所以只有general和vocation两个指标。ses是社会经济水平，也是个三分类变量(low, middle, high)，这里只有sesmiddle和seshigh是因为把seslow作为了reference。write是一个连续变量。下面的数字都是系数，需要用exp()函数进行指数变换，才能得到我们想要的OR值。下面的Std.Errors就是标准误，反映模型拟合程度。
如何去解释上面这个模型的结果？首先要先把所有的系数进行exp也就是指数变换。
exp(coef(test))
# (Intercept) sesmiddle seshigh write
# general 17.32582 0.5866769 0.3126026 0.9437172
# vocation 184.61262 1.3382809 0.3743123 0.8926116
0.5866769就是说当社会地位ses从low变成middle时，学生选择general的可能性要比academic高0.5866769倍。1.3382809则是说当社会地位ses从low变成middle时，学生选择general的可能性要比vocation高0.5866769倍。（以此类推）
上面的模型只展示了系数值，没有P值，需要手动计算。
z <- summary(test)$coefficients/summary(test)$standard.errors
z
# (Intercept) sesmiddle seshigh write
# general 2.445214 -1.2018081 -2.261334 -2.705562
# vocation 4.484769 0.6116747 -1.649967 -5.112689
p <- (1-pnorm(abs(z),0,1))*2
p
# (Intercept) sesmiddle seshigh write
# general 0.0144766100 0.2294379 0.02373856 6.818902e-03
# vocation 0.0000072993 0.5407530 0.09894976 3.176045e-07
3. 有序多分类logistic回归
有序多分类logistic回归用于因变量为有序多分类的情况，如患者对药物的反应y共有三种等级：疗效差、一般和好。定义y=1(疗效差)、y=2(疗效一般)、y=3(疗效好)。
对于有序多分类logistic回归，模型将因变量的多个分类依次分割为多个二元logistic回归。如3种不同疗效的患者，分析时拆分为2个二元logistic回归，分别为 (1vs 2 3) 、(1 2 vs3)
logit(p1) =α1+β1x1+β2x2+…
logit(p1+p2) =α2+β1x1+β2x2+…
可以看到，有序多分类logistic回归的假设是，拆分后的几个二元logistic回归的自变量系数相等，仅常数项不等。因此，有序多分类的logistic回归模型中，必须进行平行线检验，即检验自变量系数是否相等；如果不满足平行线检验，可采用无序多分类logistic回归。
做有序多分类logistic回归可以使用
MASS包中的polr函数。
例：下表中因变量y是大学生申请大学的可能性大小（有序三分类）。
require（foreign)
require(MASS)
dat <- read.dta("https://stats.idre.ucla.edu/stat/data/ologit.dta")
head(dat)
# apply pared public gpa
# 1 very likely 0 0 3.26
# 2 somewhat likely 1 0 3.21
# 3 unlikely 1 1 3.94
# 4 somewhat likely 0 0 2.81
# 5 somewhat likely 0 0 2.53
# 6 unlikely 0 1 2.59
#粗略看一下数据
ftable(xtabs(~public+apply+pared,data = dat))
# pared 0 1
# public apply
# 0 unlikely 175 14
# somewhat likely 98 26
# very likely 20 10
# 1 unlikely 25 6
# somewhat likely 12 4
# very likely 7 3
m <- polr(apply ~ pared + public + gpa, data = dat, Hess=TRUE,method='logistic')
summary(m)
提取所有的系数
ctable <- coef(summary(m))
ctable
# Value Std. Error t value
# pared 1.04769010 0.2657894 3.9418050
# public -0.05878572 0.2978614 -0.1973593
# gpa 0.61594057 0.2606340 2.3632399
# unlikely|somewhat likely 2.20391473 0.7795455 2.8271792
# somewhat likely|very likely 4.29936315 0.8043267 5.3452947
计算p值和系数表格整合
p <- pnorm(abs(ctable[,'t value']),lower.tail=FALSE)*2
ctable <- cbind(ctable,'p value'=p)
ctable
# Value Std. Error t value p value
# pared 1.04769010 0.2657894 3.9418050 8.087072e-05
# public -0.05878572 0.2978614 -0.1973593 8.435464e-01
# gpa 0.61594057 0.2606340 2.3632399 1.811594e-02
# unlikely|somewhat likely 2.20391473 0.7795455 2.8271792 4.696004e-03
# somewhat likely|very likely 4.29936315 0.8043267 5.3452947 9.027008e-08
计算置信区间
ci <- confint(m)
# Waiting for profiling to be done...
ci
# 2.5 % 97.5 %
# pared 0.5281768 1.5721750
# public -0.6522060 0.5191384
# gpa 0.1076202 1.1309148
指数变换
exp(cbind(OR=coef(m),ci))
# OR 2.5 % 97.5 %
# pared 2.8510579 1.6958376 4.817114
# public 0.9429088 0.5208954 1.680579
# gpa 1.8513972 1.1136247 3.098490
如何解释：三个自变量pared是父母教育水平，OR=2.85是说pared每改变一个单位，结局变量从第一个水平变成第二个水平/从第二个水平变成第三个水平的概率是0.85。95%置信区间不包括1，说明它是显著的。
4. 条件logistic回归
条件logistic回归针对1:n匹配的病例对照研究
在配对的数据中，病例和对照会成对存在，如1个病例匹配1个或多个对照。设有n个匹配组(每一组可视为一个层)，每一组的第一个观察对象为病例，另有m个对照，用Xitj表示第i组第t个观察对象的第j个影响因素对应的观察值。
若用Pi表示第i层在m个影响因素作用下的发病概率，那么条件logistic模型可表示为
β0i表示各匹配组的效应，不同匹配组的数值可以不同，但假定其相同；β1、β2、β3等表示被估计的总体参数。
条件logistic回归需要满足以下6个条件：
条件1：因变量为二分类变量。
条件2：至少有1个自变量，可以是分类变量，也可以是连续变量。
条件3：因变量的观察结果为配对设计或具有相关性，即不满足独立性。
条件4：因变量对子数为自变量个数的10~15倍(EPV原则)，最好＞30对，自变量的参照水平组不应少于30或50例。
条件5：自变量之间无多重共线性。
条件6：自变量不存在明显的异常值。
data(amlxray,package='faraway')
head(amlxray)
# ID disease Sex downs age Mray MupRay MlowRay Fray Cray CnRay
# 1 7004 1 F no 0 no no no no no 1
# 2 7004 0 F no 0 no no no no no 1
# 3 7006 1 M no 6 no no no no yes 3
# 4 7006 0 M no 6 no no no no yes 2
# 5 7009 1 F no 8 no no no no no 1
# 6 7009 0 F no 8 no no no no no 1
str(amlxray)
# 'data.frame': 238 obs. of 11 variables:
# $ ID : Factor w/ 112 levels "7004","7006",..: 1 1 2 2 3 3 4 4 5 5 ...
# $ disease: num 1 0 1 0 1 0 1 0 1 0 ...
# $ Sex : Factor w/ 2 levels "F","M": 1 1 2 2 1 1 2 2 2 2 ...
# $ downs : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 2 1 1 1 ...
# $ age : int 0 0 6 6 8 8 1 1 4 4 ...
# $ Mray : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 2 1 ...
# $ MupRay : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
# $ MlowRay: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 2 1 ...
# $ Fray : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
# $ Cray : Factor w/ 2 levels "no","yes": 1 1 2 2 1 1 1 1 1 2 ...
# $ CnRay : Ord.factor w/ 4 levels "1"<"2"<"3"<"4": 1 1 3 2 1 1 1 1 1 2 ...
需要使用strata(ID)来指明哪个参数是匹配的
library(survival)
cmod <- clogit(disease~Sex+Mray+Fray+CnRay+strata(ID),data=amlxray)
summary(cmod)
