前言

    接着上期Logistic回归的学习，这期我将主要介绍无序多分类Logistic回归及有序Logistic回归。内容主要包括：

1. 无序及有序多分类Logistic回归使用情形

2. 基本统计学原理及统计假设

3. R实现及可视化

4. 模型变量选择，模型检验、预测、评价等

概述

    不同于二分类logistic回归，有序多分类logistic回归(或无序多分类)适用于应变量为等级、程度差别的资料(或无序分类资料)，它基于累积概率构建回归模型。有序logistic回归模型可表示为：

   Y为一个等级变量包括g个类别，因此有序多分类logistic回归有g-1个方程，从公式中我们还可以注意到各方程回归系数相同，不同方程差别主要体现在常数项上，这主要是因为我们在做有序logistics回归前有个平行性假设，假定自变量g-1个模型中对累积概率的优势比（OR值）影响相同，即不同累积logit回归线相互平行，也叫成比例发生比假设。

    与有序logistics回归不同，无序多分类logistics回归方程中自变量回归系数可以不同，模型可表示为：

    因此无序多分类logistics回归可以说是二分类logistics回归的扩展。

有序Logistic回归

    为方便理解，我们还是引入一个例子来说明：

    研究者希望通过收集毕业学生性别、父母教育状况(low\high)、毕业时的学校性质(public\private)以及gpa成绩，了解(预测)学生对待申请毕业去向学校的态度(unlikely\somewhat likely\very likely)  。 显然态度倾向属于等级类型 ，为有序资料。

数据读取

require(foreign)
require(nnet)
require(ggplot2)
require(reshape2)
dat <- read.dta("http://www.ats.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

平行性假设检验

    有序Logistic回归分析主要用到MASS包和VGAM包，但平行性假设只有VGAM包可以完成，我们首先运用VGAM包实现平行性假设，然后通过平行性假设的原理，构建一个可视化的平行性验证方法。

用VGAM包进行平行性检验

library(VGAM)
#建立模型#parallel = TRUE,指模型满足平行性假设
model1 <- vglm(ordered(apply) ~ pared + public +gpa,data = dat,family = cumulative(parallel = TRUE))
summary(model1)

## Call:

## vglm(formula = ordered(apply) ~ pared + public + gpa, family = cumulative(parallel = TRUE), 

##     data = dat)

## Pearson residuals:

##                   Min      1Q Median     3Q    Max

## logit(P[Y<=1]) -1.671 -1.1309 0.6756 0.8247 1.8059

## logit(P[Y<=2]) -4.057  0.1806 0.2040 0.4554 0.7532

## Coefficients:

##               Estimate Std. Error z value Pr(>|z|)    

## (Intercept):1  2.20335    0.78440   2.809  0.00497 ** 

## (Intercept):2  4.29879    0.80915   5.313 1.08e-07 ***

## pared         -1.04766    0.26845  -3.903 9.52e-05 ***

## public         0.05867    0.28861   0.203  0.83891    

## gpa           -0.61575    0.26258  -2.345  0.01903 *  

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Number of linear predictors:  2 

## Names of linear predictors: logit(P[Y<=1]), logit(P[Y<=2])

## Dispersion Parameter for cumulative family:   1

## Residual deviance: 717.0249 on 795 degrees of freedom

## Log-likelihood: -358.5124 on 795 degrees of freedom

## Number of iterations: 4 

## Exponentiated coefficients:

##     pared    public       gpa 

## 0.3507593 1.0604268 0.5402335

#VGAM包可以直接给出各变量P值#parallel = FALSE,指模型不满足平行性假设，即拟合多项ligistic回归模型
model2 <- vglm(ordered(apply) ~ pared + public +gpa,data = dat,family = cumulative(parallel = FALSE))
#建立假设#H0：平行性模型与非平行性模型拟合效果一致#H1：平行性模型与非平行性模型拟合效果不一致#似然比检验平行性假设
lrtest(model1,model2)

## Likelihood ratio test

## Model 1: ordered(apply) ~ pared + public + gpa

## Model 2: ordered(apply) ~ pared + public + gpa

##   #Df  LogLik Df  Chisq Pr(>Chisq)

## 1 795 -358.51                     

## 2 792 -356.51 -3 4.0138       0.26

#结果表明，p>0.05,满足平行性假设，model1成立
is.parallel(model1)

## (Intercept)       pared      public         gpa 

##       FALSE        TRUE        TRUE        TRUE

    通过parallel = TRUE和FALSE,构建2个满足平行性假设和不满足平行性假设的模型，建立假设，通过lrtest()似然比检验平行性假设（两模型的一致性）。

#结果表明，p>0.05,满足平行性假设，model1成立.运用is.parallel(model1)可查看结果的平行性。

运用绘图方法查看模型平行性

    运用Hmisc包创建函数sf(),计算大于或等于目标变量(如y>=2)，累积概率log值；再通过summary()构建一个结果的表格。我们也可以构建一个如以y>=2，二分类Logistic回归模型，来验证表格数据的准确性。当pared1=1时，代入glm方程，-0.378+1.1438=0.7658，即表格中pared为1，y>=2对应的表格值。分别计算每行数据，累积概率P(y>=3)-p(y>=2)值，即得到变量分别为0，1时的log差值，通过判断差值相差距离，可视化平行性假设（差值距离相近）。由图可见，pared及gpa在图中“+”距离相近，而public的“+”相距较远，说明pared及gpa满足平行性，而public不满足平行性假设。

require(MASS)
require(Hmisc)
sf <- function(y) {

  c('Y>=1' = qlogis(mean(y >= 1)),    'Y>=2' = qlogis(mean(y >= 2)),    'Y>=3' = qlogis(mean(y >= 3)))}
(s <- with(dat, summary(as.numeric(apply) ~ pared + public + gpa, fun=sf)))

glm(I(as.numeric(apply) >= 2) ~ pared, family="binomial", data = dat)

## Call:  glm(formula = I(as.numeric(apply) >= 2) ~ pared, family = "binomial", 

##     data = dat)

## Coefficients:

## (Intercept)        pared  

##     -0.3783       1.1438  

## Degrees of Freedom: 399 Total (i.e. Null);  398 Residual

## Null Deviance:       550.5 

## Residual Deviance: 534.1     AIC: 538.1

glm(I(as.numeric(apply) >= 3) ~ pared, family="binomial", data = dat)

## Call:  glm(formula = I(as.numeric(apply) >= 3) ~ pared, family = "binomial", 

##     data = dat)

## Coefficients:

## (Intercept)        pared  

##      -2.441        1.094  

## Degrees of Freedom: 399 Total (i.e. Null);  398 Residual

## Null Deviance:       260.1 

## Residual Deviance: 252.2     AIC: 256.2

s[, 4] <- s[, 4] - s[, 3]
s[, 3] <- s[, 3] - s[, 3]
s 

plot
plot(s, which=1:3, pch=1:3, xlab='logit', main=' ', xlim=range(s[,3:4]))

参考资料

• idre Institute for Digital Research and Education:Data Analysis Exampleshttp://www.ats.ucla.edu/stat/dae/

• R软件-医学统计分析-有序logistic回归 http://wenku.baidu.com/link?url=5qliS_YiD2X2UH82wmbryB5GIvBhBsR5yh6z-MDtuC1NTJTV5cd6QmOhPKWsD_zTm9IL7Jpg71zL07w5_-D6S1Am7SM_aFF5f6a028aLUju

• 王汉生. (2008). 应用商务统计分析

• 其他相关资料