广义线性模型
梁雪枫
2014年12月22日
广义线性模型（GLM）是正态线性模型的直接推广，使用于连续数据和离散数据。广义线性模型通过链接函数，将因变量的期望与线性自变量相联系，通过误差函数描述误差的分布，使得许多线性模型的方法能被用于一般的问题。下表为广义线性模型中常见的链接函数和误差函数。
| 连接函数 | 典型误差函数 | |
|---|---|---|
| 恒等 | \(x^{T}\beta =E(y)\) | 正态分布 |
| 对数 | \(x^{T}\beta =lnE(y)\) | Poisson分布 |
| Logit | \(x^{T}\beta =logitE(y)\) | 二项分布 |
| 逆 | \(x^{T}\beta =\frac{1}{E(y)}\) | Gamma分布 |
R中常用的分布族和连接函数见下表
| 分布族 | 连接函数 | 默认连接函数 |
|---|---|---|
| binomial | logit,probit,cloglog | link=“logit” |
| gaussian | identity | link=“identity” |
| Gamma | identity,inverse,log | link=“inverse” |
| inverse.gaussian | 1/mu^2 | link=“1/mu^2”“ |
| poisson | identity,log,sqrt | link=“log”“ |
| quasi | logit,probit,cloglog,identity, | link=“identity”,variance=“constant” |
| inverse,log,1/mu^2,sqrt |
R中通过glm(formula, family=family.generator,data=data.frame) 函数用来做广义线性回归。正态分布族的使用方法：glm(formula, family = gaussian(link = identity),data = data.frame) , link指定了连接函数，正态分布族的连接函数缺省值是恒等的，link = identity可以不写。分布族缺省值是正态分布，family = gaussian也可以不写。glm(formula,data=data.frame)与lm(formula,data=data.frame)等价。本章重点关注常用的两种模型：Logistic回归和Poisson回归。
Logistic回归
因变量为二分类或多分类时，Logistics回归是非常重要的模型。Logistics回归由于对资料的正态性和方差齐性不做要求、对自变量类型也不做要求等，使得Logistic回归模型在医学研究各个领域被广泛用，可探索某疾病的危险因素，根据危险因素预测某疾病发生的概率，等等。例如，想探讨胃癌发生的危险因素，可以选择两组人群，一组是胃癌组，一组是非胃癌组，两组人群肯定有不同的体征和生活方式等。这里的因变量就是是否胃癌，即“是”或“否”，为两分类变量，自变量就可以包括很多了，例如年龄、性别、饮食习惯、幽门螺杆菌感染等。自变量既可以是连续的，也可以是分类的。通过logistic回归分析，就可以大致了解到底哪些因素是胃癌的危险因素。Logistics回归模型的表达形式为 \[logit(P)=ln(\frac{P}{1-P})=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\cdots +\beta_{p}X_{p}\] P为暴露于某种状态下的结局概率。\(logit(P)\)是一种变量变换方式，表示对P进行logit变换。\(beta_{i}\)为偏回归系数，表示在其他自变量不变的条件下，\(X_{i}\)每变化一个单位\(logit(P)\)的估计值。对P进行了\(logit(P)\)变换后，\(ln(\frac{P}{1-P})\)的值可以取任意值。Logistics回归是通过最大似然估计（maximum likelihood estimation,MLE）求解常数项和偏回归系数,基本思想时当从总体中随机抽取n个样本后，最合理的参数估计量应该使得这n个样本观测值的概率最大。最大似然法的基本思想是先建立似然函数与对数似然函数，再通过使对数似然函数最大求解相应的参数值，所得到的估计值称为参数的最大似然估计值。
在R语言中，进行logistic回归的命令是通过广义线性模型进行的：fm <- glm(formula, family = binomial(link = logit),data=data.frame) 。Logistic回归的基本方法是极大似然方法，其前提是样本较大。在样本量较小、数据结构较偏时，可以用精确Logistic回归（Exact logistic regression）来解决这一问题，该方法通过建立条件似然函数，进一步求出参数的充分统计量的分布函数。随着计算方法的发展和优化，也出现了使用马尔可夫链蒙特卡罗算法来模拟精确Logistic回归。R语言中的elrm包就可以实现这种算法。glm()拟合二项模型时对于因变量，如果是向量， 则假定操作二元（binary）数据，因此要求是0/1向量。 如果因变量是双列矩阵，则假定第一列为试验成功的次数第二列为试验失败的次数。如果因变量是因子，则第一水平作为失败 (0) 考虑而其他的作为`成功’(1) 考虑。
单因素Logistics回归
某项研究观察一种基因对于胃癌的诊断价值，选择了115名胃癌患者和115名非胃癌患者，检测他们的基因表达状态，欲分析该基因对胃癌是否有一定的诊断价值。
| 胃癌 | 基因+ | 基因- |
|---|---|---|
| 是 | 50 | 65 |
| 否 | 4 | 111 |
本例研究的因变量为二分类变量，分析基因的影响可以用\(\chi ^{2}\)和Logistics回归。
x<-c(50, 4, 65, 111)
dim(x)<-c(2,2)
chisq.test(x,correct = FALSE)##
## Pearson's Chi-squared test
##
## data: x
## X-squared = 51.208, df = 1, p-value = 8.308e-13P值较小，可以认为该基因对胃癌的诊断具有统计学意义。
#将表转化为扁平格式
table2flat <- function(mytable){
df <- as.data.frame(mytable)
rows <- dim(df)[1]
cols <- dim(df)[2]
x <- NULL
for(i in 1:rows){
for(j in 1:as.integer(as.character(df$Freq[i]))){
row <- df[i,c(1:(cols-1))]
x <- rbind(x,row)
}
}
row.names(x) <- c(1:dim(x)[1])
return(x)
}
gene <- rep(c(1,0),times=2)
cancer <- rep(c("0","1"),each=2)
Freq <- c(50,4,65,111)
mytable <- as.data.frame(cbind(gene,cancer,Freq))
mydata <- table2flat(mytable)
#绘制条件密度图，查看线性关系
cdplot(cancer~gene,data=mydata)
fit.glm<- glm(cancer~gene,family = binomial, data = mydata)
summary(fit.glm)##
## Call:
## glm(formula = cancer ~ gene, family = binomial, data = mydata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5918 0.2661 0.2661 1.0682 1.0682
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.3232 0.5089 6.530 6.58e-11 ***
## gene1 -3.0609 0.5426 -5.641 1.69e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 250.70 on 229 degrees of freedom
## Residual deviance: 192.19 on 228 degrees of freedom
## AIC: 196.19
##
## Number of Fisher Scoring iterations: 6回归诊断
拟合优度（goodness of fit）
拟合优度度量的是预测值和观测值之间的总体一致性。但是在评价模型时，实际上测量的是预测值和观测值之间的差别，也就是实际上检测的是模型预测的“劣度”不是”优度“，即拟合不佳检验 （lack of fit test）常用的两个指标是 Hosmer-Lemeshow指标（HL）和信息测量指标（information measure）(IM). Hosmer Lemeshow拟合优度指标(通常简写为H-L),对应的统计假设\(H_{0}\)是预测值概率和观测值之间无显著差异，所以如果HL指标显示较大的P-value，说明统计结果不显著，因此，不能拒绝关于模型拟合数据很好的假设，换句话说，模型很好的拟合了数据。 IM指标中比较常用的是AIC，在其他条件不变的情况下，较小的AIC值表示拟合模型较好。
模型卡方统计（model chi-square statistic）
模型卡方统计检测的是模型中所包含的统计量对因变量有显著的解释能力，也就是说所设模型比零假设模型（即只包含常数项的模型）要好，在多元线性回归和ANOVA中，常用F检验达到目的。在logistic中用似然比检验（likelihood ratio test）,相当于F检验。需要注意的是，模型卡方值和拟合优度是两个完全不同的概念：模型卡方值度量的是自变量是否与因变量的odds自然对数线性相关，而拟合优度度量的是预测值与观测值之间的一致性。所以按照理想情况，最好是模型的卡方检验统计性显著而拟合优度的统计性不显著。如果发生不一致，实践中更优先关注前者。
预测准确性
模型卡方统计关注的只是对于零假设模型而言，所设模型显著不显著，它只是从总体上考虑了模型的显著性，但是所有X变量到底能解释多少Y变量的波动？这是预测准确性的问题，有两种方法：(1)类RSQUARE指标：在线性回归中，可以用RSQUARE来度量，显然RSQUARE越高说明预测越好，在logistic中，也有类似的指标。logistic中的RSQUARE也有许多重要的性质：与经典的RSQUARE定义一致，它可以被理解为Y变异中被解释的比例。(2)AUC值(C统计量)：拟合优度只是给出了观测值和预测概率直接的差别程度，然后给出了一个总体评价的指标，但是在实际应用中，往往更关心观测值和模型预测的条件事件概率的关联强度，这类指标被称为序列相关指标，指标值越高，表示预测概率与观测反应变量直接的关联越密切。通常用ROC图来和ROC图的曲线下面积（AUC）进行，AUC可以定量地评价模型的效果，AUC越大则模型效果越好。ROC曲线下的面积值在1.0和0.5之间。在AUC>0.5的情况下，AUC越接近于1，说明诊断效果越好。AUC在 0.5～0.7时有较低准确性，AUC在0.7～0.9时有一定准确性，AUC在0.9以上时有较高准确性。AUC=0.5时，说明诊断方法完全不起作用，无诊断价值。AUC<0.5不符合真实情况，在实际中极少出现。大于或等于0.75一般认为认为模型是可靠的。
ROC（receiver operating characteristic curve，受试者工作特征曲线）曲线，横轴是1-Specificity（特异度），纵轴是Sensitivity（灵敏度）。45度线是作为参照（baseline model）出现的，就是说，ROC的好坏，乃是跟45度线相比的。选择最佳的诊断界限值。ROC曲线越靠近左上角,试验的准确性就越高。最靠近左上角的ROC曲线的点是错误最少的最好阈值，其假阳性和假阴性的总数最少。两种或两种以上不同诊断试验对疾病识别能力的比较。在对同一种疾病的两种或两种以上诊断方法进行比较时，可将各试验的ROC曲线绘制到同一坐标中，以直观地鉴别优劣，靠近左上角的ROC曲线所代表的受试者工作最准确。亦可通过分别计算各个试验的ROC曲线下的面积(AUC)进行比较，哪一种试验的 AUC最大，则哪一种试验的诊断价值最佳。
对于0-1变量的二分类问题，分类的最终结果可以用表格表示为：
| 预测值0 | 预测值1 | |
|---|---|---|
| 实际值0 | a | b |
| 实际值1 | c | d |
其中，d是“实际为1而预测为1”的样本个数，c是“实际为1而预测为0”的样本个数，其余依此类推。显然地，主对角线所占的比重越大，则预测效果越佳，这也是一个基本的评价指标——总体准确率(a+d)/(a+b+c+d)。TPR（真阳性率、灵敏度）：True Positive Rate，将实际的1正确地预测为1的概率，d/(c+d)。FPR：False Positive Rate（假阳性率，1-特异度），将实际的0错误地预测为1的概率，b/(a+b)。TPR与FPR相互影响的重要因素就是“阈值”。当阈值为0时，所有的样本都被预测为正例，因此TPR=1，而FPR=1。此时的FPR过大，无法实现分类的效果。随着阈值逐渐增大，被预测为正例的样本数逐渐减少，TPR和FPR各自减小，当阈值增大至1时，没有样本被预测为正例，此时TPR=0，FPR=0。
统计量最为关注的是AUC值，其次是似然卡方统计量，然后才是HL统计量，对AIC 和RSQUARE 极少关注，这一点和多元线性回归有很大的不同，根本原因是多元线性回归是一个预测模型，目标变量的值具有实际的数值意义；而logistic是一个分类模型，目标变量的值是一个分类标识，因此更关注观测值和预测值之间的相对一致性，而不是绝对一致性。
rms包lrm()函数可以计算相关统计量。
model <- lrm(cancer~gene,data=mydata)
#Nagelkerke等其他拟合优度指标
goodfit <- function(glmFit)
{
N <- nobs(glmFit)
glm0 <- update(glmFit, . ~ 1)
LLf <- logLik(glmFit)
LL0 <- logLik(glm0)
McFadden <- as.vector(1 - (LLf / LL0))
CoxSnell <- as.vector(1 - exp((2/N) * (LL0 - LLf)))
Nagelkerke <- as.vector((1 - exp((2/N) * (LL0 - LLf))) / (1 - exp(LL0)^(2/N)))
result <- list(McFadden=McFadden, CoxSnell= CoxSnell,Nagelkerke=Nagelkerke)
return (result)
}
model$stats## Obs Max Deriv Model L.R. d.f. P
## 2.300000e+02 2.212457e-06 5.850576e+01 1.000000e+00 2.031708e-14
## C Dxy Gamma Tau-a R2
## 7.783039e-01 5.566077e-01 9.104991e-01 2.008734e-01 3.383625e-01
## Brier g gr gp
## 1.396597e-01 1.537119e+00 4.651169e+00 2.008734e-01goodfit(fit.glm)## $McFadden
## [1] 0.2333735
##
## $CoxSnell
## [1] 0.2245974
##
## $Nagelkerke
## [1] 0.3383625AUC值（C统计量）为0.778，可以认为模型比较可靠。似然比检验结果比较显著，观测值和预测值的一致性有差异。HL统计量,在多因素统计中予以计算。Deducer包rocplot可以绘制ROC曲线并计算AUC值。CoxSnell\(R_{2}\)系数与线性回归分析中的决定系数\(R_{2}\)有相似指出，也是回归方程对因变量变异解释程度的反应，由于CoxSnell\(R_{2}\)系数取值范围不易确定，不易直接判断拟合效果。Nagelkerke\(R_{2}\)系数是对CoxSnell\(R_{2}\)的修正，取值范围在0~1之间，越接近于1，说明模型的拟合优度越高。但对Logistic回归而言，伪决定系数不像线性回归中决定系数那么重要。
rocplot(fit.glm)
#### 影响分析 对于异常值识别仍然可用influence.measures()函数获得。
influencePlot(fit.glm)## Warning in plot.window(...): relative range of values = 51 * EPS, is small
## (axis 1)
## StudRes Hat CookD
## 51 -2.638313 0.008695652 0.3503971influence.measures(fit.glm)## Influence measures of
## glm(formula = cancer ~ gene, family = binomial, data = mydata) :
##
## dfb.1_ dfb.gen1 dffit cov.r cook.d hat inf
## 1 4.57e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 2 8.94e-15 -0.0459 -0.1325 1.000 0.005752 0.0087
## 3 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 4 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 5 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 6 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 7 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 8 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 9 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 10 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 11 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 12 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 13 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 14 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 15 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 16 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 17 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 18 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 19 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 20 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 21 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 22 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 23 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 24 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 25 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 26 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 27 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 28 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 29 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 30 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 31 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 32 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 33 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 34 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 35 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 36 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 37 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 38 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 39 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 40 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 41 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 42 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 43 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 44 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 45 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 46 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 47 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 48 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 49 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 50 1.47e-16 -0.0459 -0.1325 1.000 0.005752 0.0087
## 51 -2.70e-01 0.2530 -0.2698 0.947 0.122778 0.0087 *
## 52 -2.70e-01 0.2530 -0.2698 0.947 0.122778 0.0087 *
## 53 -2.70e-01 0.2530 -0.2698 0.947 0.122778 0.0087 *
## 54 -2.70e-01 0.2530 -0.2698 0.947 0.122778 0.0087 *
## 55 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 56 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 57 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 58 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 59 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 60 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 61 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 62 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 63 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 64 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 65 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 66 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 67 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 68 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 69 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 70 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 71 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 72 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 73 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 74 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 75 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 76 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 77 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 78 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 79 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 80 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 81 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 82 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 83 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 84 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 85 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 86 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 87 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 88 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 89 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 90 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 91 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 92 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 93 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 94 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 95 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 96 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 97 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 98 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 99 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 100 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 101 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 102 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 103 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 104 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 105 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 106 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 107 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 108 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 109 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 110 -1.05e-16 0.0380 0.1095 1.006 0.003403 0.0087
## 111 -1.05e-16 0.0380 0.1095 1.006 0.003403 0.0087
## 112 -1.05e-16 0.0380 0.1095 1.006 0.003403 0.0087
## 113 -1.05e-16 0.0380 0.1095 1.006 0.003403 0.0087
## 114 -1.05e-16 0.0380 0.1095 1.006 0.003403 0.0087
## 115 -1.05e-16 0.0380 0.1095 1.006 0.003403 0.0087
## 116 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 117 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 118 -9.44e-17 0.0380 0.1095 1.006 0.003403 0.0087
## 119 -1.05e-16 0.0380 0.1095 1.006 0.003403 0.0087
## 120 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 121 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 122 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 123 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 124 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 125 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 126 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 127 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 128 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 129 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 130 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 131 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 132 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 133 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 134 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 135 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 136 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 137 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 138 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 139 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 140 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 141 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 142 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 143 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 144 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 145 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 146 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 147 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 148 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 149 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 150 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 151 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 152 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 153 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 154 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 155 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 156 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 157 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 158 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 159 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 160 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 161 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 162 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 163 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 164 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 165 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 166 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 167 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 168 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 169 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 170 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 171 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 172 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 173 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 174 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 175 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 176 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 177 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 178 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 179 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 180 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 181 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 182 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 183 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 184 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 185 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 186 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 187 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 188 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 189 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 190 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 191 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 192 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 193 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 194 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 195 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 196 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 197 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 198 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 199 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 200 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 201 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 202 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 203 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 204 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 205 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 206 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 207 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 208 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 209 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 210 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 211 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 212 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 213 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 214 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 215 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 216 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 217 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 218 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 219 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 220 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 221 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 222 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 223 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 224 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 225 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 226 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 227 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 228 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 229 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087
## 230 2.72e-02 -0.0255 0.0272 1.017 0.000159 0.0087多重共线性
可用vif（方差膨胀因子）进行判断，vif开平方是否大于2，若大于2，则存在多重共线性问题。
过度离散
因变量的方差大于期望的二项分布的方差，过度离散会导致奇异标准误检验和不精确的的显著性检验。检验过度离散的一种方法是比较二项分布模型的残差偏差与残差自由度，如果比值比1大很多，可以认为存在过度离散。对过度离散的假设检验需要用family = “quasibinomial”再进行一次模型拟合。
overdispersion <- function(fit.glm){
Phi <- fit.glm$deviance/fit.glm$df.residual
fit.od <- glm(fit.glm$formula,family = quasibinomial, data = fit.glm$data)
p <- pchisq(summary(fit.od)$dispersion*fit.glm$df.residual,fit.glm$df.residual,lower=F)
return (list(Phi=Phi,p.value=p))
}
overdispersion(fit.glm)## $Phi
## [1] 0.8429389
##
## $p.value
## [1] 0.4504226比值在1附近，并且P值大于0.05，不能拒绝比值为1的假设，可以认为不存在过度离散。
模型参数解释
logistic.display(fit.glm)##
## Logistic regression predicting cancer : 1 vs 0
##
## OR(95%CI) P(Wald's test) P(LR-test)
## gene: 1 vs 0 0.05 (0.02,0.14) < 0.001 < 0.001
##
## Log-likelihood = -96.095
## No. of observations = 230
## AIC value = 196.1901exp(coef(fit.glm))即为OR值，表示自变量增加一个单位，因变量则乘以OR值。OR值具有风险的含义，在危险因素研究中具有重要意义。LR-test（likelihood ration test）为似然比检验,Wald’s test为Wadld检验，P值小于0.05均说明回归系数具有统计学意义，自变量与因变量有统计学联系。OR值大于1，为危险因素。OR值小于1，为保护因素。
多因素Logistics回归
AER包中包含一个Affairs数据，记录了一组婚外情数据，其中包括参与者性别、年龄、婚龄、是否有小孩、宗教信仰程度（5分制，1表示反对，5表示非常信仰）、学历、职业和婚姻的自我评分（5分制，1表示非常不幸福，5表示非常幸福）。
data(Affairs,package="AER")
Affairs$ynaffair[Affairs$affairs > 0] <- 1
Affairs$ynaffair[Affairs$affairs==0] <- 0
Affairs$ynaffair <- factor(Affairs$ynaffair,levels=c(0,1),labels=c("No","Yes"))与线性回归相似，bestglm包中bestglm函数可以完成logistic回归的全子集的自变量筛选。
Affairs <- Affairs[,c("gender","age","yearsmarried","children","religiousness","education","occupation","rating","ynaffair")]
best.logistic <-bestglm(Affairs,family = binomial,IC = "AIC",method = "exhaustive") ## Morgan-Tatar search since family is non-gaussian.best.logistic$BestModels## gender age yearsmarried children religiousness education occupation
## 1 TRUE TRUE TRUE FALSE TRUE FALSE FALSE
## 2 TRUE TRUE TRUE TRUE TRUE FALSE FALSE
## 3 FALSE TRUE TRUE TRUE TRUE FALSE TRUE
## 4 FALSE TRUE TRUE FALSE TRUE FALSE FALSE
## 5 FALSE TRUE TRUE TRUE TRUE FALSE FALSE
## rating Criterion
## 1 TRUE 621.8590
## 2 TRUE 622.1529
## 3 TRUE 623.2897
## 4 TRUE 623.3578
## 5 TRUE 623.4076summary(best.logistic$BestModel)##
## Call:
## glm(formula = y ~ ., family = family, data = Xi, weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5623 -0.7495 -0.5664 -0.2671 2.3975
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.94760 0.61234 3.181 0.001470 **
## gendermale 0.38612 0.20703 1.865 0.062171 .
## age -0.04393 0.01806 -2.432 0.015011 *
## yearsmarried 0.11133 0.02983 3.732 0.000190 ***
## religiousness -0.32714 0.08947 -3.656 0.000256 ***
## rating -0.46721 0.08928 -5.233 1.67e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 675.38 on 600 degrees of freedom
## Residual deviance: 611.86 on 595 degrees of freedom
## AIC: 623.86
##
## Number of Fisher Scoring iterations: 4也可用逐步法完成自变量的筛选。
fit.full <- glm(ynaffair~gender+age+yearsmarried+children+religiousness+education+occupation+rating,data=Affairs,family=binomial())
step(fit.full)## Start: AIC=627.51
## ynaffair ~ gender + age + yearsmarried + children + religiousness +
## education + occupation + rating
##
## Df Deviance AIC
## - education 1 609.68 625.68
## - occupation 1 609.70 625.70
## - gender 1 610.89 626.89
## - children 1 611.40 627.40
## <none> 609.51 627.51
## - age 1 615.67 631.67
## - yearsmarried 1 618.34 634.34
## - religiousness 1 622.92 638.92
## - rating 1 636.75 652.75
##
## Step: AIC=625.68
## ynaffair ~ gender + age + yearsmarried + children + religiousness +
## occupation + rating
##
## Df Deviance AIC
## - occupation 1 610.15 624.15
## - gender 1 611.29 625.29
## - children 1 611.62 625.62
## <none> 609.68 625.68
## - age 1 615.78 629.78
## - yearsmarried 1 618.46 632.46
## - religiousness 1 623.27 637.27
## - rating 1 636.93 650.93
##
## Step: AIC=624.15
## ynaffair ~ gender + age + yearsmarried + children + religiousness +
## rating
##
## Df Deviance AIC
## - children 1 611.86 623.86
## <none> 610.15 624.15
## - gender 1 613.41 625.41
## - age 1 616.00 628.00
## - yearsmarried 1 619.07 631.07
## - religiousness 1 623.98 635.98
## - rating 1 637.23 649.23
##
## Step: AIC=623.86
## ynaffair ~ gender + age + yearsmarried + religiousness + rating
##
## Df Deviance AIC
## <none> 611.86 623.86
## - gender 1 615.36 625.36
## - age 1 618.05 628.05
## - religiousness 1 625.57 635.57
## - yearsmarried 1 626.23 636.23
## - rating 1 639.93 649.93##
## Call: glm(formula = ynaffair ~ gender + age + yearsmarried + religiousness +
## rating, family = binomial(), data = Affairs)
##
## Coefficients:
## (Intercept) gendermale age yearsmarried religiousness
## 1.94760 0.38612 -0.04393 0.11133 -0.32714
## rating
## -0.46721
##
## Degrees of Freedom: 600 Total (i.e. Null); 595 Residual
## Null Deviance: 675.4
## Residual Deviance: 611.9 AIC: 623.9全自集和逐步法对自变量的筛选均应建立在对自变量专业考虑的基础上进行。本例中两种方法结果类似，对其结果进行诊断。
fit <- glm(ynaffair ~ gender + age + yearsmarried + religiousness + rating, family = binomial(), data = Affairs)
lrm(ynaffair ~ gender + age + yearsmarried + religiousness + rating,data=Affairs)##
## Logistic Regression Model
##
## lrm(formula = ynaffair ~ gender + age + yearsmarried + religiousness +
## rating, data = Affairs)
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 601 LR chi2 63.52 R2 0.149 C 0.707
## No 451 d.f. 5 g 0.896 Dxy 0.414
## Yes 150 Pr(> chi2) <0.0001 gr 2.451 gamma 0.416
## max |deriv| 7e-07 gp 0.156 tau-a 0.155
## Brier 0.166
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept 1.9476 0.6123 3.18 0.0015
## gender=male 0.3861 0.2070 1.87 0.0622
## age -0.0439 0.0181 -2.43 0.0150
## yearsmarried 0.1113 0.0298 3.73 0.0002
## religiousness -0.3271 0.0895 -3.66 0.0003
## rating -0.4672 0.0893 -5.23 <0.0001AUC=0.70,模型尚可。
fit0 <- glm(formula = ynaffair ~ 1, family = binomial(), data = Affairs)
anova(fit0,fit,test="Chisq")## Analysis of Deviance Table
##
## Model 1: ynaffair ~ 1
## Model 2: ynaffair ~ gender + age + yearsmarried + religiousness + rating
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 600 675.38
## 2 595 611.86 5 63.518 2.274e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1模型的likelihood-ratio检验，P值小于0.05，可以认为模型的自变量与因变量的odds自然对数线性相关。
Hosmer-Lemeshowz指标
hosmerlem <- function(y, yhat, g=10) {
cutyhat = cut(yhat, breaks = quantile(yhat, probs=seq(0,1,1/g)), include.lowest=TRUE)
obs = xtabs(cbind(1 - y, y) ~ cutyhat)
expect = xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
chisq = sum((obs - expect)^2/expect)
P = 1 - pchisq(chisq, g - 2)
return(list(chisq=chisq,p.value=P))
}
hosmerlem(y=Affairs$ynaffair, yhat=fitted(fit))## Warning in Ops.factor(1, y): '-' not meaningful for factors## $chisq
## [1] 601
##
## $p.value
## [1] 0Hosmer Lemeshow拟合优度指标检验P值小于0.05，可以认为预测值和观测值之间差异显著。 过度离散诊断
overdispersion(fit)## $Phi
## [1] 1.028334
##
## $p.value
## [1] 0.3291751过度离散检验P值大于0.05，可以认为不存在过度离散。
多重共线性
sqrt(vif(fit)) > 2## gendermale age yearsmarried religiousness rating
## FALSE FALSE FALSE FALSE FALSE自变量之间不存在多重共线性。可用influence.measures()函数进行影响分析，logistic.display()对自变量进行解释。
稳健Logistic回归
robust包中的glmRob（）函数可用来拟合稳健的广义线性模型，包括稳健Logistic回归；当拟合回归模型数据出现离群点和强影响点时，便可应用稳健Logistic回归。对influence.measures(fit)进行影响分析后，发现存在强影响点，应用稳健Logistic回归。
fit.rob <- glmRob(ynaffair ~ gender + age + yearsmarried + religiousness + rating, family = binomial(), data = Affairs)
summary(fit.rob)##
## Call: glmRob(formula = ynaffair ~ gender + age + yearsmarried + religiousness +
## rating, family = binomial(), data = Affairs)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5981 -0.7458 -0.5617 -0.2632 2.4162
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.10089 0.62033 3.387 7.073e-04
## gendermale 0.41140 0.20873 1.971 4.872e-02
## age -0.04438 0.01842 -2.409 1.599e-02
## yearsmarried 0.10972 0.03012 3.643 2.693e-04
## religiousness -0.32271 0.09022 -3.577 3.475e-04
## rating -0.50900 0.09028 -5.638 1.724e-08
##
## (Dispersion Parameter for binomial family taken to be 1 )
##
## Null Deviance: 833.2 on 600 degrees of freedom
##
## Residual Deviance: 612.0911 on 595 degrees of freedom
##
## Number of Iterations: 6
##
## Correlation of Coefficients:
## (Intercept) gendermale age yearsmarried religiousness
## gendermale 0.0289160
## age -0.6448048 -0.2633502
## yearsmarried 0.3169328 0.2046799 -0.7641181
## religiousness -0.3586014 -0.0006552 0.0065429 -0.1706658
## rating -0.6226794 -0.0592698 0.0763706 0.0626459 -0.0032677条件logistic回归
条件logistic回归假设自变量在各配对组中对结果变量的作用是相同的，即自变量的回归系数与配对组无关。配对设计的Logistic回归模型不含常数项，参数估计是根据条件概率得到的。对病例和对照进行配比能控制影响实验效应的主要非处理因素，可以提高统计分析的效能，通常可分为1:1，1:n，m:n配对。epicalc包中的VC1to1来自于验证吸烟、酗酒和橡胶行业工作是否是食管癌的危险因素的病例对照研究。
data(VC1to1,package = "epicalc")
pander(VC1to1)| matset | case | smoking | rubber | alcohol |
|---|---|---|---|---|
| 1 | 1 | 1 | 0 | 0 |
| 1 | 0 | 1 | 0 | 0 |
| 2 | 1 | 1 | 0 | 1 |
| 2 | 0 | 1 | 1 | 0 |
| 3 | 1 | 1 | 1 | 0 |
| 3 | 0 | 1 | 1 | 0 |
| 4 | 1 | 1 | 0 | 0 |
| 4 | 0 | 1 | 1 | 1 |
| 5 | 1 | 0 | 0 | 1 |
| 5 | 0 | 1 | 0 | 0 |
| 6 | 1 | 1 | 0 | 1 |
| 6 | 0 | 0 | 0 | 0 |
| 7 | 1 | 1 | 0 | 1 |
| 7 | 0 | 1 | 0 | 0 |
| 8 | 1 | 1 | 0 | 0 |
| 8 | 0 | 1 | 0 | 0 |
| 9 | 1 | 1 | 1 | 1 |
| 9 | 0 | 1 | 1 | 0 |
| 10 | 1 | 0 | 0 | 0 |
| 10 | 0 | 1 | 1 | 0 |
| 11 | 1 | 1 | 0 | 1 |
| 11 | 0 | 1 | 0 | 1 |
| 12 | 1 | 1 | 0 | 0 |
| 12 | 0 | 1 | 0 | 1 |
| 13 | 1 | 1 | 1 | 0 |
| 13 | 0 | 0 | 0 | 0 |
| 14 | 1 | 1 | 0 | 1 |
| 14 | 0 | 1 | 0 | 1 |
| 15 | 1 | 1 | 0 | 1 |
| 15 | 0 | 1 | 0 | 0 |
| 16 | 1 | 1 | 0 | 1 |
| 16 | 0 | 0 | 0 | 1 |
| 17 | 1 | 1 | 1 | 1 |
| 17 | 0 | 1 | 0 | 1 |
| 18 | 1 | 1 | 0 | 1 |
| 18 | 0 | 1 | 0 | 1 |
| 19 | 1 | 0 | 0 | 0 |
| 19 | 0 | 1 | 0 | 0 |
| 20 | 1 | 1 | 1 | 1 |
| 20 | 0 | 0 | 0 | 0 |
| 21 | 1 | 1 | 1 | 1 |
| 21 | 0 | 1 | 1 | 1 |
| 22 | 1 | 0 | 0 | 1 |
| 22 | 0 | 1 | 1 | 1 |
| 23 | 1 | 1 | 1 | 1 |
| 23 | 0 | 1 | 1 | 1 |
| 24 | 1 | 1 | 1 | 1 |
| 24 | 0 | 1 | 0 | 0 |
| 25 | 1 | 0 | 0 | 0 |
| 25 | 0 | 1 | 1 | 0 |
| 26 | 1 | 1 | 0 | 1 |
| 26 | 0 | 0 | 0 | 0 |
use(VC1to1)
matchTab(case,smoking,strat=matset)##
## Exposure status: smoking = 1
##
## Total number of match sets in the tabulation = 26
##
## Number of controls = 1
## No. of controls exposed
## No. of cases exposed 0 1
## 0 0 5
## 1 5 16
##
## Odds ratio by Mantel-Haenszel method = 1
##
## Odds ratio by maximum likelihood estimate (MLE) method = 1
## 95%CI= 0.29 , 3.454case变量1表示患病，0表示未患病。matset变量表示对子号。epicalc包中matchTab()函数用以计算条件优势比（McNemar’s优势比），表示病例间不一致部分的计数比值，其95%置信区间如果包含1，则表示变量没有统计学意义。
fit.c <- clogit(case~smoking+alcohol+rubber+strata(matset),data=VC1to1,method = "exact")
summary(fit.c)## Call:
## coxph(formula = Surv(rep(1, 52L), case) ~ smoking + alcohol +
## rubber + strata(matset), data = VC1to1, method = "exact")
##
## n= 52, number of events= 26
##
## coef exp(coef) se(coef) z Pr(>|z|)
## smoking 0.04321 1.04416 0.86916 0.050 0.9604
## alcohol 1.66701 5.29629 0.82878 2.011 0.0443 *
## rubber -0.68078 0.50622 0.94518 -0.720 0.4714
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## smoking 1.0442 0.9577 0.1901 5.736
## alcohol 5.2963 0.1888 1.0435 26.880
## rubber 0.5062 1.9754 0.0794 3.228
##
## Rsquare= 0.101 (max possible= 0.5 )
## Likelihood ratio test= 5.55 on 3 df, p=0.1355
## Wald test = 4.11 on 3 df, p=0.25
## Score (logrank) test = 5.05 on 3 df, p=0.168fit.c$loglik## [1] -18.02183 -15.24534survival包中clogit()函数可以完成条件Logistic回归，结果显示模型与空模型比较，差异无显著性。自变量smoking和rubber均无显著性差异，自变量alcohol差异显著。条件Logistic回归模型不能得到对数似然比和AIC值，但能得到条件对数似然比，以表示模型的拟合水平。
无序多分类Logistic回归
若因变量包含两个以上的无序类别（比如，已婚/寡居/离婚），便可使用mlogit包中的mlogit（）函数拟合多项Logistic回归。epicalc中Ectopic数据集，其中outc变量中Deci表示正常分娩，IA表示发生人工流产，EP表示发生宫外孕。hia变量表示以前是否有IA（人工流产史），gravi表示怀孕的次数。
data(Ectopic,package = "epicalc")
pander(head(Ectopic))| id | outc | hia | gravi |
|---|---|---|---|
| 1 | Deli | ever IA | 1-2 |
| 2 | Deli | ever IA | 3-4 |
| 3 | Deli | never IA | 1-2 |
| 4 | Deli | never IA | 1-2 |
| 5 | Deli | never IA | 1-2 |
| 6 | IA | ever IA | 1-2 |
ep <- Ectopic$outc=="EP"
ia <- Ectopic$outc=="IA"
deli <- Ectopic$outc=="Deli"
mnFit <- multinom(cbind(deli,ep,ia)~hia+gravi, data=Ectopic)## # weights: 15 (8 variable)
## initial value 794.296685
## iter 10 value 745.073806
## final value 744.587307
## convergedsummary(mnFit)## Call:
## multinom(formula = cbind(deli, ep, ia) ~ hia + gravi, data = Ectopic)
##
## Coefficients:
## (Intercept) hiaever IA gravi3-4 gravi>4
## ep -1.0194479 1.491331 0.4659169 0.695505
## ia -0.5088417 0.382695 0.8527850 1.164605
##
## Std. Errors:
## (Intercept) hiaever IA gravi3-4 gravi>4
## ep 0.154495 0.2216795 0.2398888 0.3658957
## ia 0.130889 0.2146635 0.2367876 0.3691153
##
## Residual Deviance: 1489.175
## AIC: 1505.175mlogit.display(mnFit)##
## Outcome =cbind(deli, ep, ia); Referent group = deli
## ep ia
## Coeff./SE RRR(95%CI) Coeff./SE RRR(95%CI)
## (Intercept) -1.02/0.154*** - -0.51/0.131*** -
## hiaever IA 1.49/0.222*** 4.44(2.88,6.86) 0.38/0.215 1.47(0.96,2.23)
## gravi3-4 0.47/0.24 1.59(1,2.55) 0.85/0.237*** 2.35(1.48,3.73)
## gravi>4 0.7/0.366 2(0.98,4.11) 1.16/0.369** 3.2(1.55,6.61)
##
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual Deviance: 1489.17
## AIC = 1505.17vglmFitMN <- vglm(outc~hia+gravi, family=multinomial(refLevel=3), data=Ectopic)
exp(VGAM::coef(vglmFitMN))## (Intercept):1 (Intercept):2 hiaever IA:1 hiaever IA:2 gravi3-4:1
## 0.3607945 0.6011905 4.4429622 1.4662256 1.5934754
## gravi3-4:2 gravi>4:1 gravi>4:2
## 2.3461652 2.0047063 3.2046317dfMNL <- mlogit.data(Ectopic, choice="outc", shape="wide", varying=NULL)
mlogitFit <- mlogit(outc ~ 0 | hia+gravi,, reflevel="Deli", data=dfMNL)
summary(mlogitFit)##
## Call:
## mlogit(formula = outc ~ 0 | hia + gravi, data = dfMNL, reflevel = "Deli",
## method = "nr", print.level = 0)
##
## Frequencies of alternatives:
## Deli EP IA
## 0.33333 0.33333 0.33333
##
## nr method
## 4 iterations, 0h:0m:0s
## g'(-H)^-1g = 4.96E-08
## gradient close to zero
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## EP:(intercept) -1.01945 0.15449 -6.5986 4.151e-11 ***
## IA:(intercept) -0.50884 0.13089 -3.8876 0.0001012 ***
## EP:hiaever IA 1.49132 0.22168 6.7274 1.727e-11 ***
## IA:hiaever IA 0.38269 0.21466 1.7828 0.0746266 .
## EP:gravi3-4 0.46592 0.23989 1.9422 0.0521099 .
## IA:gravi3-4 0.85278 0.23679 3.6015 0.0003164 ***
## EP:gravi>4 0.69550 0.36589 1.9008 0.0573263 .
## IA:gravi>4 1.16460 0.36911 3.1551 0.0016044 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -744.59
## McFadden R^2: 0.062583
## Likelihood ratio test : chisq = 99.419 (p.value = < 2.22e-16)exp(mlogitFit$coefficients)## EP:(intercept) IA:(intercept) EP:hiaever IA IA:hiaever IA EP:gravi3-4
## 0.3607945 0.6011905 4.4429622 1.4662256 1.5934754
## IA:gravi3-4 EP:gravi>4 IA:gravi>4
## 2.3461652 2.0047063 3.2046317
## attr(,"fixed")
## EP:(intercept) IA:(intercept) EP:hiaever IA IA:hiaever IA EP:gravi3-4
## FALSE FALSE FALSE FALSE FALSE
## IA:gravi3-4 EP:gravi>4 IA:gravi>4
## FALSE FALSE FALSEnnet包中multinom()函数、VGAM包中的vglm()函数、mlogit包中mlogit()函数均得到了相似的结果。mlogit()函数对数据格式于其他两个函数的要求有所不同，其中formula：mlogit提供了条件logit，多项logit，混合logit多种模型，对于多项logit的估计模型应写为：因变量~0|自变量,data：使用mlogit.data函数使得数据结构符合mlogit函数要求。Choice：确定分类变量是什么Shape：如果每一行是一个观测，选择wide，如果每一行是表示一个选择，那么选择long。alt.var：对于shape为long的数据，需要标明所有的选择名称。由于mlogit包可以做的logit模型更多。
本例中是以outc变量的Deli(分娩)做为参考水平的，有人工流产史的病例(hia ever IA)发生宫外孕（EP）的危险增加4.44，有人工流产史的病例(hia ever IA)发生人工流产（IA）的危险增加1.47(置信区间包括1，无显著性意义)。multinom()函数默认是第一水平，可通过levels(Ectopic$outc)方法查看。vglm()和mlogit()函数是可以指定参考水平。
模型拟合评价
PhatCateg <- VGAM::predict(vglmFitMN, type="response")
categHat <- levels(Ectopic$outc)[max.col(PhatCateg)]
facHat <- factor(categHat,levels=levels(Ectopic$outc))
cTab <- xtabs(~ outc+ facHat, data=Ectopic)
addmargins(cTab)## facHat
## outc EP IA Deli Sum
## EP 180 10 51 241
## IA 131 29 81 241
## Deli 83 14 144 241
## Sum 394 53 276 723CCR <- sum(diag(cTab)) / sum(cTab)
CCR## [1] 0.4882434上述方法可获得模型的正确分类率，本例的正确分类率为0.4882434，正确分类率偏低。
偏差、对数似然值和AIC值
deviance <- VGAM::deviance(vglmFitMN)
logLik<- VGAM::logLik(vglmFitMN)
AIC <- VGAM::AIC(vglmFitMN)
deviance## [1] 1489.175logLik## [1] -744.5873AIC## [1] 1505.175McFadden, Cox & Snell and Nagelkerke \(R^{2}\)伪决定系数
vglm()函数拟合结果并没有直接给出伪决定系数，可通过如下方法计算相关统计量。
vglm0 <- vglm(outc~ 1, family=multinomial(refLevel=3), data=Ectopic)
LLf <- VGAM::logLik(vglmFitMN)
LL0 <- VGAM::logLik(vglm0)
N <- nobs(vglmFitMN)
McFadden <- as.vector(1 - (LLf / LL0))
CoxSnell<- as.vector(1 - exp((2/N) * (LL0 - LLf)))
Nagelkerke<- as.vector((1 - exp((2/N) * (LL0 - LLf))) / (1 - exp(LL0)^(2/N)))
McFadden## [1] 0.06258288CoxSnell## [1] 0.1284732Nagelkerke## [1] 0.1284732Nagelkerke伪决定系数为0.1284732，表明自变量对因变量的解释程度不高。
系数及模型的检验
vglm函数结果中并没有系数及模型的检验情况。对模型的系数及其95%置信区间可从如下方法获得。
sumMN <- VGAM::summary(vglmFitMN)
coefMN <- VGAM::coef(sumMN)
zCrit <- qnorm(c(0.05/2, 1 - 0.05/2))
ciCoef <- t(apply(coefMN, 1, function(x) x["Estimate"] - zCrit*x["Std. Error"] ))
coefMN## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -1.0194468 0.1544949 -6.598579 4.151179e-11
## (Intercept):2 -0.5088434 0.1308891 -3.887593 1.012432e-04
## hiaever IA:1 1.4913213 0.2216792 6.727385 1.727388e-11
## hiaever IA:2 0.3826915 0.2146632 1.782753 7.462647e-02
## gravi3-4:1 0.4659174 0.2398885 1.942225 5.210985e-02
## gravi3-4:2 0.8527822 0.2367874 3.601468 3.164251e-04
## gravi>4:1 0.6954976 0.3658890 1.900843 5.732255e-02
## gravi>4:2 1.1645972 0.3691088 3.155160 1.604101e-03ciCoef## [,1] [,2]
## (Intercept):1 -0.7166424 -1.322251265
## (Intercept):2 -0.2523055 -0.765381185
## hiaever IA:1 1.9258045 1.056838098
## hiaever IA:2 0.8034236 -0.038040613
## gravi3-4:1 0.9360901 -0.004255348
## gravi3-4:2 1.3168769 0.388687455
## gravi>4:1 1.4126268 -0.021631628
## gravi>4:2 1.8880371 0.441157314似然比检验通过如下方法获得。
vglm0 <- vglm(outc~ 1, family=multinomial(refLevel=3), data=Ectopic)
VGAM::lrtest(vglmFitMN, vglm0)## Likelihood ratio test
##
## Model 1: outc ~ hia + gravi
## Model 2: outc ~ 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 1438 -744.59
## 2 1444 -794.30 6 99.419 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1似然比检验结果表明含有两个自变量的模型和仅有截距项的模型相比有显著性差异。对系数的检验结果表明有人工流产史的病例(hia ever IA)是发生宫外孕(EP)危险因素。
预测分类
vglm拟合结果可通过如下的方法得到每个分类的预测概率。
PhatCateg <- VGAM::predict(vglmFitMN, type="response")
head(PhatCateg)## EP IA Deli
## 1 0.4600392 0.2529737 0.2869871
## 2 0.4543112 0.3678299 0.1778589
## 3 0.1838926 0.3064195 0.5096879
## 4 0.1838926 0.3064195 0.5096879
## 5 0.1838926 0.3064195 0.5096879
## 6 0.4600392 0.2529737 0.2869871还可以通过如下两种方法分别得到针对multinom()、mlogit()每个分类的预测概率。
predict(mnFit, type="probs")
fitted(mlogitFit, outcome=FALSE)对分类结果的预测有如下两种方法。
PhatCateg <- VGAM::predict(vglmFitMN, type="response")
categHat <- levels(Ectopic$outc)[max.col(PhatCateg)]
head(categHat)## [1] "EP" "EP" "Deli" "Deli" "Deli" "EP"predCls <- predict(mnFit, type="class")
head(predCls)有序多分类Logistic回归
若因变量是一有序的类别（比如，无效/有效/显效/控制），使用无序多分类Logistic回归处理因变量，不但会丧失变量间联系的功效，而且会曲解因变量和自变量之间的相关方式。程序包MASS提供polr()函数、ordinal提供clm()函数、rms提供orm()函数、VGAM提供vglm()函数可以进行ordered logit或probit回归。累积Logistic回归模型(cumulative logit model)如下，\({logit}(p(Y \geq g)) = \ln \frac{P(Y \geq g)}{1 - P(Y \geq g)} = \beta_{0_{g}} + \beta_{1} X_{1} + \dots + \beta_{p} X_{p} \quad(g = 2, \ldots, k)\)。成比例比数比累计Logistic模型(proportional-adds cumulative logit mode)简化上述模型，使自变量\(X_{i}\)所对应的回归系数\(\beta_{i}\)都是相等。在此假设条件下，不同累计Logistic的回归线相互平行，只是截距\(\beta_{i}\)不同。 例 epicalc中HW93数据集是1993年泰国南部钩虫感染的调查资料，其中intense变量表示感染的严重程度为有序多分类变量，shoes表示是否穿鞋，agegr是年龄分组。
data(HW93,package = "epicalc")
intense.ord <- ordered(HW93$intense)在自变量较多的时候，可以采用R中自动逐步变量筛选step()函数，仅MASS包中polr()函数能够支持自变量的筛选。
polrFit <- polr(intense.ord~agegr+shoes,method="logistic",data=HW93)
exp(MASS:::confint.polr(polrFit))## Waiting for profiling to be done...##
## Re-fitting to get Hessian## 2.5 % 97.5 %
## agegr15-59 yrs 1.5173279 3.1156330
## agegr60+ yrs 1.9133925 6.7875609
## shoesyes 0.3413853 0.6862894ordinal.or.display(polrFit)##
## Re-fitting to get Hessian## Waiting for profiling to be done...##
## Re-fitting to get Hessian## Ordinal OR lower95ci upper95ci P value
## agegr15-59 yrs 2.169 1.517 3.116 1.39e-05
## agegr60+ yrs 3.596 1.913 6.788 4.07e-05
## shoesyes 0.485 0.341 0.686 2.71e-05VGAM包
vglmFit <- vglm(intense.ord~agegr+shoes, family=propodds, data=HW93)VGAM包能进行所有类型的logistic回归的计算，并且能进行累计Logistic回归模型的平行性假设检验，其他包则不能。模型中family=cumulative(parallel=TRUE, reverse=TRUE)指定拟合累计Logistic回归模型，而且parallel=T指定模型按平行性假定进行拟合，该选项可简写为amily=propodds。
vglm(intense.ord~agegr+shoes, family=cumulative(parallel=TRUE, reverse=TRUE),data=HW93)
vglm(intense.ord~agegr+shoes, family=acat(parallel=TRUE), data=HW93)
vglm(intense.ord~agegr+shoes, family=sratio(parallel=TRUE), data=HW93)rms包
ormFit <- orm(intense~agegr+shoes, data=HW93)ordinal包
clmFit <- clm(intense~agegr+shoes, link="logit", data=HW93)结果显示，上述有序多分类Logisitic回归模型有两个截距，每一个都是结果的一个切割点，这些截距项的值没有实际意义。年龄的系数通过两个切割点进行了分割，两个系数均为正数表示危险度随年龄的增加而增加，穿鞋的系数为负数表示穿鞋对两种感染水平均有保护作用。
模型评价
vglmFit <- vglm(intense.ord~agegr+shoes, family=propodds, data=HW93)
PhatCateg <- VGAM::predict(vglmFit, type="response")
categHat <- levels(HW93$intense)[max.col(PhatCateg)]
facHat <- factor(categHat, levels=levels(HW93$intense))
cTab <- xtabs(~ intense + facHat, data=HW93)
addmargins(cTab)## facHat
## intense 0 1-1,999 2,000+ Sum
## 0 13 184 0 197
## 1-1,999 20 329 0 349
## 2,000+ 4 87 0 91
## Sum 37 600 0 637(CCR <- sum(diag(cTab)) / sum(cTab))## [1] 0.5368917上述方法可获得模型的正确分类率，本例的正确分类率为0.5368917，正确分类率偏低。 偏差、对数似然值和AIC值
deviance <- VGAM::deviance(vglmFit)
logLik<- VGAM::logLik(vglmFit)
AIC <- VGAM::AIC(vglmFit)
deviance## [1] 1204.92logLik## [1] -602.4602AIC## [1] 1214.92McFadden, Cox & Snell and Nagelkerke \(R^{2}\)伪决定系数
vglm0 <- vglm(intense.ord~ 1, family=propodds, data=HW93)
LLf <- VGAM::logLik(vglmFit)
LL0 <- VGAM::logLik(vglm0)
McFadden <- as.vector(1 - (LLf / LL0))
CoxSnell<- as.vector(1 - exp((2/N) * (LL0 - LLf)))
Nagelkerke<- as.vector((1 - exp((2/N) * (LL0 - LLf))) / (1 - exp(LL0)^(2/N)))
McFadden## [1] 0.02555969CoxSnell## [1] 0.0427723Nagelkerke## [1] 0.05221336系数及模型的检验
sumOrd <- summary(vglmFit)
coefOrd <- coef(sumOrd)
exp(coefOrd[,1])## (Intercept):1 (Intercept):2 agegr15-59 yrs agegr60+ yrs shoesyes
## 1.8778574 0.1256241 2.1694020 3.5956164 0.4850639zCrit <- qnorm(c(0.05/2, 1 - 0.05/2))
ciCoef <- t(apply(coefOrd, 1, function(x) x["Estimate"] - zCrit*x["Std. Error"] ))MASS包建立的模型可直接使用confint()函数计算OR值及其可信区间。
summary(polrFit)##
## Re-fitting to get Hessian## Call:
## polr(formula = intense.ord ~ agegr + shoes, data = HW93, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## agegr15-59 yrs 0.7744 0.1834 4.222
## agegr60+ yrs 1.2797 0.3227 3.966
## shoesyes -0.7235 0.1780 -4.064
##
## Intercepts:
## Value Std. Error t value
## 0|1-1,999 -0.6301 0.1293 -4.8726
## 1-1,999|2,000+ 2.0745 0.1579 13.1363
##
## Residual Deviance: 1204.92
## AIC: 1214.92exp(cbind(OR=coef(polrFit),t(confint(polrFit))))## Warning in cbind(OR = coef(polrFit), t(confint(polrFit))): number of rows
## of result is not a multiple of vector length (arg 1)## OR agegr15-59 yrs agegr60+ yrs shoesyes
## 2.5 % 2.169385 1.517328 1.913392 0.3413853
## 97.5 % 3.595648 3.115633 6.787561 0.6862894ordinal包建立的模型用summary()函数即可输出系数。
summary(clmFit)## formula: intense ~ agegr + shoes
## data: HW93
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 637 -602.46 1214.92 5(0) 1.65e-08 2.6e+01
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## agegr15-59 yrs 0.7745 0.1834 4.222 2.42e-05 ***
## agegr60+ yrs 1.2797 0.3227 3.966 7.30e-05 ***
## shoesyes -0.7235 0.1780 -4.064 4.82e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1-1,999 -0.6301 0.1293 -4.873
## 1-1,999|2,000+ 2.0745 0.1579 13.136模型比较
vglmR <- vglm(intense.ord~ shoes, family=propodds, data=HW93)
VGAM::lrtest(vglmFit, vglmR)## Likelihood ratio test
##
## Model 1: intense.ord ~ agegr + shoes
## Model 2: intense.ord ~ shoes
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 1269 -602.46
## 2 1271 -615.25 2 25.587 2.779e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1VGAM::lrtest(vglmFit, vglm0)## Likelihood ratio test
##
## Model 1: intense.ord ~ agegr + shoes
## Model 2: intense.ord ~ 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 1269 -602.46
## 2 1272 -618.26 3 31.605 6.338e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1vglmFit与其他两个模型比较均有显著性差异，选择LogLik值较大的，还有两个自变量的模型。选择更优模型还可以比较两个模型的信息统计量AIC和BIC，信息统计量小的模型更优。
AIC(vglmR)## [1] 1236.507AIC(vglm0)## [1] 1240.526AIC(vglmFit)## [1] 1214.92平行性假设检验
为了检验平行性假设，需要建立非平行的模型，将平行性模型与非平行性模型进行似然比检验，检验平行性假设
vglmP <- vglm(intense.ord~agegr+shoes, family=cumulative(parallel=TRUE, reverse=TRUE),data=HW93)
vglmNP <- vglm(intense.ord~agegr+shoes, family=cumulative(parallel=FALSE, reverse=TRUE),data=HW93)
VGAM::lrtest(vglmP, vglmNP)## Likelihood ratio test
##
## Model 1: intense.ord ~ agegr + shoes
## Model 2: intense.ord ~ agegr + shoes
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 1269 -602.46
## 2 1266 -598.44 -3 8.0455 0.04508 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1clmP <- clm(intense~agegr+shoes, link="logit", data=HW93)
clmNP <- clm(intense~shoes, nominal=~agegr, data=HW93)
anova(clmP, clmNP)## Likelihood ratio tests of cumulative link models:
##
## formula: nominal: link: threshold:
## clmP intense ~ agegr + shoes ~1 logit flexible
## clmNP intense ~ shoes ~agegr logit flexible
##
## no.par AIC logLik LR.stat df Pr(>Chisq)
## clmP 5 1214.9 -602.46
## clmNP 7 1212.1 -599.07 6.7849 2 0.03363 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1平行性假设检验结果表明，P值小于0.05，可以认为平行性假设不成立。检验结果可用is.parallel()函数获得。 ### 精确Logistic回归 例 elrm包的drugDat数据集记录不同性别人群在某种药物治疗的结果，recovered表示恢复数量，n表示总人数。
data(drugDat,package = "elrm")
pander(drugDat)| sex | treatment | recovered | n |
|---|---|---|---|
| 1 | 1 | 16 | 27 |
| 0 | 1 | 10 | 19 |
| 1 | 0 | 13 | 32 |
| 0 | 0 | 7 | 21 |
data(drugDat)
drug.elrm=elrm(formula=recovered/n~sex+treatment,interest=~sex+treatment,iter=100000,burnIn=1000,dataset=drugDat)
summary(drug.elrm)Possion回归
Poisson回归的因变量是计数型的变量，自变量是连续性或类别型变量。Poisson回归因变量通常局限在一个固定长度时间段内进行测量（如过去一年交通事故数），整个观测集中时间长度都是不变的。Poisson回归主要有两个假设，首先，具有相同特征和同时的不同对象的人时风险是同质的，其次，当样本量越来越大时，频数的均数趋近于方差。
例 robust包中Breslow癫痫数据记录了治疗初期八周内，抗癫痫药物对癫痫发病数的影响，因变量sumY为随机后8周内癫痫发病数，自变量治疗Trt，年龄Age和治疗前8周的癫痫发病数Base。
data(breslow.dat,package="robust")
opar <- par(no.readonly=T)
par(mfrow = c(1,2))
attach(breslow.dat)
hist(sumY,breaks = 20,xlab = "Seizure Count",main="Distribution of Seizure")
boxplot(sumY~Trt,xlab="Treatment",main="Group Coomparisons")
par(opar)从图中可以清楚的看到因变量的偏移特性及可能的离群点。药物治疗下癫痫的发病数似乎变小，且方差也变小了。
fit <- glm(sumY~Base+Age+Trt,data=breslow.dat,family = poisson())
summary(fit)##
## Call:
## glm(formula = sumY ~ Base + Age + Trt, family = poisson(), data = breslow.dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.0569 -2.0433 -0.9397 0.7929 11.0061
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.9488259 0.1356191 14.370 < 2e-16 ***
## Base 0.0226517 0.0005093 44.476 < 2e-16 ***
## Age 0.0227401 0.0040240 5.651 1.59e-08 ***
## Trtprogabide -0.1527009 0.0478051 -3.194 0.0014 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2122.73 on 58 degrees of freedom
## Residual deviance: 559.44 on 55 degrees of freedom
## AIC: 850.71
##
## Number of Fisher Scoring iterations: 5VGAM包vglm()方法获得类似结果
vglmFit <- vglm(sumY~Base+Age+Trt, family=poissonff, data=breslow.dat)
summary(vglmFit)结果输出了偏差、回归参数、标准误和参数为0的检验。
拟合优度检验
检验建立Poisson模型的拟合优度
poisgof(fit)## $results
## [1] "Goodness-of-fit test for Poisson assumption"
##
## $chisq
## [1] 559.4437
##
## $df
## [1] 55
##
## $p.value
## [1] 1.203254e-84P值较小，表明模型的拟合优度较差。
模型的系数及解释
exp(coef(fit))## (Intercept) Base Age Trtprogabide
## 7.0204403 1.0229102 1.0230007 0.8583864Base、Age、Trt和截距项检验均显示有意义，在保持其他变量不变，年龄增加1岁，癫痫发病数将乘以1.023。一单位的Trt变化（从安慰剂到治疗组），癫痫发病数将乘以0.86,也就是说治疗组想对于安慰剂组发病数下降了。危险比的95%置信区间可通过
idr.display(fit)##
## Poisson regression predicting sumY
##
## crude IDR(95%CI) adj. IDR(95%CI)
## Base (cont. var.) 1.02 (1.02,1.02) 1.02 (1.02,1.02)
##
## Age (cont. var.) 0.99 (0.98,1) 1.02 (1.01,1.03)
##
## Trt: progabide vs placebo 0.93 (0.85,1.01) 0.86 (0.78,0.94)
##
## P(Wald's test) P(LR-test)
## Base (cont. var.) < 0.001 < 0.001
##
## Age (cont. var.) < 0.001 < 0.001
##
## Trt: progabide vs placebo 0.001 0.001
##
## Log-likelihood = -421.3535
## No. of observations = 59
## AIC value = 850.7071过度离散
与Logistic回归类似，如果残差的偏差和和残差的自由度之比大于1,那么表明存在过度离散。Poisson分布的方差和均数相等，当因变量的方差比预测方差大时，Poisson分布可能会发生过度离散。过度离散可能会对结果的解释造成影响，可能会得到很小的标准误和置信区间，并且显著性检验也比较宽松。发生过度离散可能是遗漏了某个重要变量或者是计数事件并不独立。过度离散检验可用qcc包的qcc.overdispersion.test()方法。
qcc.overdispersion.test(breslow.dat$sumY,type = "poisson")##
## Overdispersion test Obs.Var/Theor.Var Statistic p-value
## poisson data 62.87013 3646.468 0P值小于0.05,表明确实存在过度离散。通过用family=“quasipoisson” 替换family=“poisson”,以完成对过度离散数据的拟合。
fit.od <- glm(sumY~Base+Age+Trt,data=breslow.dat,family = quasipoisson())
summary(fit.od)##
## Call:
## glm(formula = sumY ~ Base + Age + Trt, family = quasipoisson(),
## data = breslow.dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.0569 -2.0433 -0.9397 0.7929 11.0061
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.948826 0.465091 4.190 0.000102 ***
## Base 0.022652 0.001747 12.969 < 2e-16 ***
## Age 0.022740 0.013800 1.648 0.105085
## Trtprogabide -0.152701 0.163943 -0.931 0.355702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 11.76075)
##
## Null deviance: 2122.73 on 58 degrees of freedom
## Residual deviance: 559.44 on 55 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5VGAM包vglm()方法获得类似结果
vglm <- vglm(sumY~Base+Age+Trt, family=quasipoissonff, data=breslow.dat)
summary(vglm)使用类Poisson方法估计的参数与Poisson相同，但标准误变大。当考虑过度离散，Base、Trt和Age均没有显著意义。
异方差一致的标准误差
可通过如下方法获得
hcSE <- vcovHC(fit, type="HC0")
coeftest(fit, vcov=hcSE)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.9488259 0.3651482 5.3371 9.445e-08 ***
## Base 0.0226517 0.0012357 18.3316 < 2.2e-16 ***
## Age 0.0227401 0.0115804 1.9637 0.04957 *
## Trtprogabide -0.1527009 0.1711089 -0.8924 0.37217
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1时间段变化的Poisson回归
当观测时间长度不同时，可以拟合时间段变化的Poisson回归模型，次住假设结果变量是比率。为分析比率，数据中需包含一个记录每个观测时间长度的变量(如time)。然后模型将从\(ln(\lambda )=\beta _{0}+\sum_{j=1}^{p}\beta _{j}X_{j}\)修改为\(ln\begin{pmatrix}\frac{\lambda}{time} \end{pmatrix}=\beta _{0}+\sum_{j=1}^{p}\beta _{j}X_{j}\)。为拟合新模型，需要使用glm()函数中的offset选项。假设Breslow中有一个time变量，记录了病人随机分组后监测时间长度的变化，拟合模型如下
fit <- glm(sumY~Base+Age+Trt,data=breslow.dat,offset=log(time),family = poisson())
vglmFit <- vglm(sumY~Base+Age+Trt,offset=log(time), family=poissonff, data=breslow.dat)零膨胀的Poisson回归
当因变量中，0计数的数目比Poisson回归预测的数据多时，即总体的一个子群体无任何被计数的行为时，就可能发生这种问题。
set.seed(123)
N <- 200
sigma <- matrix(c(4,2,-3, 2,16,-1, -3,-1,8), byrow=TRUE, ncol=3)
mu <- c(-3, 2, 4)
XY <- rmvnorm(N, mean=mu, sigma=sigma)
Y <- round(XY[ , 3] - 1.5)
Y[Y < 0] <- 0
dfCount <- data.frame(X1=XY[ , 1], X2=XY[ , 2], Y)
ziFitP <- zeroinfl(Y ~ X1 + X2 | 1, dist="poisson", data=dfCount)
vglm(Y ~ X1 + X2, family=zipoissonff, data=dfCount)## Call:
## vglm(formula = Y ~ X1 + X2, family = zipoissonff, data = dfCount)
##
## Coefficients:
## (Intercept):1 (Intercept):2 X1 X2
## 0.456862233 1.886747933 -0.216977846 -0.002069534
##
## Degrees of Freedom: 400 Total; 396 Residual
## Log-likelihood: -402.4037稳健Poisson回归
influence.measures()对拟合的模型完成影响分析后，如存在离群点和强影响点，可用robust包中glmRob()方法拟合稳健广义线性模型。
fit.rob <- glmRob(sumY~Base+Age+Trt, family = poisson(), data=breslow.dat)
summary(fit.rob)##
## Call: glmRob(formula = sumY ~ Base + Age + Trt, family = poisson(),
## data = breslow.dat)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -42.27051 -1.62371 0.04412 0.82431 9.02635
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.869172 0.241628 7.736 1.028e-14
## Base 0.037434 0.004145 9.031 1.700e-19
## Age 0.009499 0.007450 1.275 2.023e-01
## Trtprogabide -0.280084 0.096155 -2.913 3.582e-03
##
## (Dispersion Parameter for poisson family taken to be 1 )
##
## Null Deviance: 11983 on 58 degrees of freedom
##
## Residual Deviance: 2939.072 on 55 degrees of freedom
##
## Number of Iterations: 4
##
## Correlation of Coefficients:
## (Intercept) Base Age
## Base -0.40654
## Age -0.91177 0.11421
## Trtprogabide 0.01327 -0.43165 -0.08329负二项回归(Negative binomial regression)
Poisson回归假定因变量是均数和方差相等，如果出现方差比均数大，就会形成过度离散，Poisson回归会低估预测变量的标准误。当过度离散比较明显时，指定误差项服从负二项分布，得到的负二项回归系数与Poisson回归相同，但标准误更大，结果的解释与Poisson回归相同。
例 epicalc包DHF99数据集是一实地调查的滋生蚊子幼虫的水容器的数据，因变量containers是有蚊子幼虫滋生的容器的频数，education 和viltype是可能对因变量有影响的自变量。
data(DHF99,package="epicalc")
opar <- par(no.readonly=T)
par(mfrow = c(1,2))
attach(DHF99)## The following object is masked from package:robustbase:
##
## educationhist(containers,breaks = 20)
boxplot(containers~viltype)
par(opar)
qcc.overdispersion.test(DHF99$containers,type = "poisson")##
## Overdispersion test Obs.Var/Theor.Var Statistic p-value
## poisson data NA NA NA因变量的偏移特性比较明显,因变量有缺失值，Poisson回归的过度离散情况不能够检验。
负二项回归拟合用MASS包中glm.nb()方法
glmFitNB <- glm.nb(containers ~ education + viltype, data=DHF99)
summary(glmFitNB)##
## Call:
## glm.nb(formula = containers ~ education + viltype, data = DHF99,
## init.theta = 0.3014602889, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.7995 -0.7559 -0.5958 -0.3473 2.9291
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.741697 0.211679 -3.504 0.000459 ***
## educationSecondary 0.007613 0.487305 0.016 0.987536
## educationHigh school 0.177650 0.485567 0.366 0.714468
## educationBachelor 0.137298 0.568759 0.241 0.809245
## educationOther -0.002669 0.494737 -0.005 0.995696
## viltypeurban -1.965173 0.527656 -3.724 0.000196 ***
## viltypeslum -0.678422 0.428132 -1.585 0.113055
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.3015) family taken to be 1)
##
## Null deviance: 175.25 on 298 degrees of freedom
## Residual deviance: 156.46 on 292 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 434.05
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.3015
## Std. Err.: 0.0768
##
## 2 x log-likelihood: -418.0450VGAM包中vglm()方法如下
vglmFitNB <- vglm(containers ~ education + viltype, family=negbinomial, data=DHF99)
summary(vglmFitNB)拟合优度检验
检验建立负二项回归模型的拟合优度
poisgof(glmFitNB)## $results
## [1] "Goodness-of-fit test for Poisson assumption"
##
## $chisq
## [1] 156.4599
##
## $df
## [1] 292
##
## $p.value
## [1] 1P值较大，表明模型的拟合优度较好。
模型的系数及解释
exp(coef(glmFitNB))## (Intercept) educationSecondary educationHigh school
## 0.4763050 1.0076418 1.1944077
## educationBachelor educationOther viltypeurban
## 1.1471704 0.9973346 0.1401316
## viltypeslum
## 0.5074172viltype的P值小于0.05,意义比较显著，结果解释与Poisson回归类似。对于自变量的选择可以采用step()和AIC值。
零膨胀的负二项回归回归
与Poisson回归类似，因变量中0计数的频数较多时，应采用零膨胀的负二项回归回归
ziFitNB <- zeroinfl(containers ~ education + viltype | 1, dist="negbin", data=DHF99)
summary(ziFitNB)##
## Call:
## zeroinfl(formula = containers ~ education + viltype | 1, data = DHF99,
## dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.4439 -0.4297 -0.3663 -0.2338 9.4931
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.740777 0.223940 -3.308 0.000940 ***
## educationSecondary 0.007712 0.486969 0.016 0.987365
## educationHigh school 0.177771 0.487075 0.365 0.715129
## educationBachelor 0.137002 0.560956 0.244 0.807053
## educationOther -0.002949 0.513227 -0.006 0.995415
## viltypeurban -1.965128 0.530590 -3.704 0.000213 ***
## viltypeslum -0.678531 0.428890 -1.582 0.113635
## Log(theta) -1.197902 0.273266 -4.384 1.17e-05 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.049 80.423 -0.088 0.93
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 0.3018
## Number of iterations in BFGS optimization: 63
## Log-likelihood: -209 on 9 DfVGAM包中vglm()方法如下，此方法的自变量不能为分类变量。
vglm(containers ~ village, family=zinegbinomial, data=DHF99)出用AIC比较模型外，Quang Vuong提出如果一个模型比另一个模型更接近真实的函数，那么从这个模型得到的每个个体的对数似然值也应该显著的大于从另一个模型得到的每个个体的对数似然值。pscl包vuong() 方法实现了Vuong检验。
vuong(ziFitNB, glmFitNB)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -0.157152 model2 > model1 0.43756
## AIC-corrected -1015.939022 model2 > model1 < 2e-16
## BIC-corrected -2895.360768 model2 > model1 < 2e-16Vuong 检验的统计量则成标准的N(0, 1)正态分布。Vuong 值大于1.96，则模型1好于模型2，小于-1.96，则结论相反。Vuong Test = -0.15表明两个模型同样地接近真实函数。