当响应类别有序的时候，logit可以利用类别顺序得到响应的模型。这样的模型比基线-类别模型有更简单的解释且有潜在的更高的功效。

YY 的累计概率是指 YY 落在或者低于一个特定点的概率。对于结果类别 jj，累计概率为：

P(Y≤j)=π1+...+πjP(Y\leq j)=\pi_1+...+\pi_j ，且 P(Y≤1)≤...≤P(Y≤J)=1P(Y\leq1)\leq...\leq P(Y\leq J)=1 ，累计概率的模型并不会利用最后一个概率，因为它必然等于1。

累积概率的logit为 logit[P(Y≤j)]=log(P(Y≤j)1−P(Y≤j))=log(π1+...+πjπj+1+...+πJ)logit[P(Y\leq j)]=log(\frac{P(Y\leq j)}{1-P(Y\leq j)})=log(\frac{\pi_1+...+\pi_j}{\pi_{j+1}+...+\pi_{J}})

具有比例优势特性的累计logit模型

第 jj 个累积logit模型中，类 11 到类 jj 合起来作为一个新的类别，类 j+1j+1 与类 JJ 合起来作为第二个类别。当有一个解释变量的时候，模型可以写成：

logit[P(Y≤j)]=αj+βx,j=1,...,J−1logit[P(Y\leq j)]=\alpha_j+\beta x,j=1,...,J-1 ， β\beta 描述的是 xx 对响应变量落在类 jj或者小于类jj 的类的对数优势的效应。

注意！这个式子中 β\beta 没有下标 jj ——模型假设 xx 的效应对所有 J−1J-1和累积变量都是相等的。预测变量的效应在不同累积logit的方程中都是相同的。当模型拟合的好时，它只需要一个参数而不是J−1J-1 个参数来描述 xx 的效应，这种情况下这个模型比基线-类别logit模型更容易概括和解释效应。

模型的解释可以利用累计概率和它们的补的优势比。对于 xx 的两个值 x1x_1 和 x2x_2 ，比较两个累计概率的优势比为 j|x=x_2)}{P(X\leq j | x=x_1)/P(X>j|x=x_1)}">P(X≤j|x=x2)/P(X>j|x=x2)P(X≤j|x=x1)/P(X>j|x=x1)\frac{P(X\leq j| x=x_2)/P(X>j|x=x_2)}{P(X\leq j | x=x_1)/P(X>j|x=x_1)} ，而这个优势比的对数就是在那两个 xx处的累计logit的差，它等于β(x2−x1)\beta(x_2-x_1) ，与两个 xx 之间的距离成比例。特别地，如果 x2−x1=1x_2-x_1=1 ，那么 xx每增加一个单位，响应变量在任意给定类别下的优势就应当乘以eβe^{\beta} 。在这个对数优势比当中，相同的比例常数适用于每一个累积概率，这样的特性被称为模型的比例优势假设。

一个栗子：政治意识形态与隶属党派的关系

> library(VGAM) > library(tidyr) > library(cdabookdb) > data("ideology") > ftable(ideology) Ideology VLib SLib Mod SCon VCon Gender Party Female Dem 44 47 118 23 32 Rep 18 28 86 39 48 Male Dem 36 34 53 18 23 Rep 12 18 62 45 51

因变量政治意识形态一共有5个类别，范围为很自由到很保守（VLib=Very Liberal; VCon = Very Conservative）。自变量书上选取的是政党——民主党x=1x=1 ，共和党 x=0x=0 。我们拟合有序响应变量的logit模型:

> ide_margin=margin.table(ideology,c(2,3))#忽略性别变量的边缘表 #spread(data,key,value) data就是整个数据，key表示要拆哪一列，value表示拆的那一列填什么值 > ide_margin_df=spread(as.data.frame(ide_margin),Ideology,Freq) > ide_m=vglm(cbind(VLib,SLib,Mod,SCon,VCon)~Party=="Dem",data=ide_margin_df,family=cumulative(parallel=TRUE)) > summary(ide_m) Call: vglm(formula = cbind(VLib, SLib, Mod, SCon, VCon) ~ Party == "Dem", family = cumulative(parallel = TRUE), data = ide_margin_df) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept):1 -2.46899 0.13182 -18.730 < 2e-16 *** (Intercept):2 -1.47454 0.10908 -13.517 < 2e-16 *** (Intercept):3 0.23712 0.09485 2.500 0.0124 * (Intercept):4 1.06954 0.10457 10.228 < 2e-16 *** Party == "Dem"TRUE 0.97452 0.12905 7.551 4.31e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]), logitlink(P[Y<=3]), logitlink(P[Y<=4]) Residual deviance: 3.6877 on 3 degrees of freedom Log-likelihood: -24.6201 on 3 degrees of freedom Number of Fisher scoring iterations: 3 No Hauck-Donner effect found in any of the estimates Exponentiated coefficients: Party == "Dem"TRUE 2.649882

由ML拟合的输出结果可以看出， J=5J=5个类别的响应变量一共有4个截距项，但是除了估计响应概率我们对这些截距项并不感兴趣。

独立性检验。为了检验独立性 ( H0:β=0H_0:\beta=0 ) ，最简单的是使用Wald检验。政治党派的估计效应是 β^=0.97(SE=0.13)\hat{\beta}=0.97(SE=0.13) ，Wald检验 p−value=4.31e−14p-value=4.31e^{-14}，结果显示政治党派的效应是显著的。一个更有效的方法是似然比检验，毕竟我们只有一个参数。> anova(ide_m,type="I",test="LRT") Analysis of Deviance Table (Type I tests: terms added sequentially from first to last) Model: cumulative, VGAMordinal, VGAMcategorical Links: logitlink Response: cbind(VLib, SLib, Mod, SCon, VCon) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 4 62.333 Party == "Dem" 1 58.645 3 3.688 1.888e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看出似然比统计量为 58.64558.645 ， df=1df=1 ， p−value=1.888e−14p-value=1.888e-14 ，给出了关联性极强的证据。

对 β\beta 的优势比解释。对任意固定的 jj ，民主党人的响应方向为自由方向而非保守方向的估计优势等于对于共产党人的估计优势的 exp(0.975)=2.65exp(0.975)=2.65 倍。我们可以说：“政治党派与政治意识形态之间存在着很强的关联，民主党人比共和党人的政治意识形态更趋向自由。”置信区间。 β\beta 的 95%95\% 置信区间为 0.975±1.96×0.129=(0.72,1.23)0.975\pm 1.96\times 0.129=(0.72,1.23) 。累计概率的优势比的置信区间为 (exp(0.72),exp(1.23))=(2.05,3.42)(exp(0.72),exp(1.23))=(2.05,3.42) 。民主党人在自由一端的优势至少是共和党人的两倍。模型拟合的检验。整体的拟合检验一般采用的是 X2X^2 统计量与 G2G^2统计量。当至多只有几个解释变量而且这些解释变量全部都是属性变量，而且几乎所有的单元频数不少于55时，这两个检验统计量就近似服从卡方分布。> X2=sum(resid(ide_m,type="pearson")^2) > x2_pvalue=1-pchisq(X2,df) > c(X2=X2,pvalue=x2_pvalue) X2 pvalue 3.6628363 0.3002487 > G2=deviance(ide_m) > g2_pvalue=1-pchisq(G2,df) > c(G2=G2,pvalue=g2_pvalue) G2 pvalue 3.6876784 0.2972241

可以看出这个模型的拟合还是充分的。

概率的表达。累计概率的模型表达为：P(Y≤j)=exp(αj+βx)(1+exp(αj+βx))P(Y\leq j)=\frac{exp(\alpha_j+\beta x)}{(1+exp(\alpha_j+\beta x))} ，如果想计算属于某一个政党的概率，将相邻的两个概率相减即可，例如 π^2=P^(Y=2)=P^(Y≤2)−P^(Y≤1)\hat{\pi}_2=\hat{P}(Y=2)=\hat{P}(Y\leq 2)-\hat{P}(Y\leq 1)> pred=predict(ide_m,data.frame(Party="Dem"),type="response") > pred VLib SLib Mod SCon VCon 1 0.1832507 0.1942834 0.3930548 0.1147564 0.1146547

另一个栗子：对心理健康的建模

> library(VGAM) > library(cdabookdb) > data("impairment") > head(impairment) Subject Impairment SES LifeEvents 1 1 Well 1 1 2 2 Well 1 9 3 3 Well 1 4 4 4 Well 1 3 5 5 Well 0 2 6 6 Well 1 0

数据是对弗罗里达州某县成年居民的一个随机样本的心理健康研究。心理伤害 (Impairment)是有序的，类别为（健康 well、轻度症状形成 mild symptom formation、中等症状形成 moderate symptom formation、受损 impaired）。这个例子的因变量Y=Y= 心理伤害，自变量是 X1X_1 社会经济地位（SES）与 X2X_2 重要生活事件的数量（LifeEvents）。其中社会经济地位为二元变量， 11 表示高， 00表示低。重要生活事件指生活中孩子的出生、新工作、离婚或者家庭成员的去世等。

具有比例优势形式的主效应模型（main effects model）为：

logit[P(Y≤j)]=αj+β1x1+β2x2logit[P(Y\leq j)]=\alpha_j+\beta_1x_1+\beta_2x_2

我们对模型进行拟合可以得到

> impairment_m=vglm(Impairment~SES+LifeEvents,family=cumulative(parallel=TRUE),data=impairment) > summary(impairment_m) Call: vglm(formula = Impairment ~ SES + LifeEvents, family = cumulative(parallel = TRUE), data = impairment) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept):1 -0.2819 0.6231 -0.452 0.65096 (Intercept):2 1.2128 0.6511 1.863 0.06251 . (Intercept):3 2.2094 0.7171 3.081 0.00206 ** SES 1.1112 0.6143 1.809 0.07045 . LifeEvents -0.3189 0.1194 -2.670 0.00759 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]), logitlink(P[Y<=3]) Residual deviance: 99.0979 on 115 degrees of freedom Log-likelihood: -49.5489 on 115 degrees of freedom Number of Fisher scoring iterations: 5 No Hauck-Donner effect found in any of the estimates Exponentiated coefficients: SES LifeEvents 3.0380707 0.7269742 对β\beta 的优势比解释。可以看出 β^1=1.111,β^2=−0.319\hat{\beta}_1=1.111,\hat{\beta}_2=-0.319 ，说明了从量表健康一端出发的累计概率随着生活事件的增加而减少，且在 SES=1SES=1的水平处更大。在给定生活事件得分下，高社会经济地位者低于任一固定水平的心理伤害程度的优势的估计值是低社会经济地位者的优势的估计值的e1.111=3.0e^{1.111}=3.0倍。高社会经济地位者相较于低社会经济地位者而言，所受的心理伤害程度会轻一些。对模型拟合的检验。皮尔逊X2X^2 与偏差 G2G^2 统计量只对非稀疏列联表才是有效的，这个数据中 00较多，这两个统计量就不适用了。我们可以通过比较这个模型与更复杂的模型来检验模型拟合的好坏。加入交互项之后我们可以得到一个模型，而且这个模型的MLML拟合为：> impairment_m_1=vglm(Impairment~SES*LifeEvents,family=cumulative(parallel=TRUE),data=impairment) > summary(impairment_m_1) Call: vglm(formula = Impairment ~ SES * LifeEvents, family = cumulative(parallel = TRUE), data = impairment) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept):1 0.09807 0.81102 0.121 0.90375 (Intercept):2 1.59248 0.83717 1.902 0.05714 . (Intercept):3 2.60660 0.90966 2.865 0.00416 ** SES 0.37090 1.13022 0.328 0.74279 LifeEvents -0.42045 0.19031 -2.209 0.02715 * SES:LifeEvents 0.18131 0.23611 0.768 0.44255 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]), logitlink(P[Y<=3]) Residual deviance: 98.5044 on 114 degrees of freedom Log-likelihood: -49.2522 on 114 degrees of freedom Number of Fisher scoring iterations: 5 No Hauck-Donner effect found in any of the estimates Exponentiated coefficients: SES LifeEvents SES:LifeEvents 1.4490350 0.6567529 1.1987822

logit[P^(Y≤j)]=α^j−0.420x1+0.371x2+0.181x1x2logit[\hat{P}(Y\leq j)]=\hat{\alpha}_j-0.420x_1+0.371x_2+0.181x_1x_2 。生活事件的优势比效应估计对于低社会经济地位组而言为 −0.420-0.420 ，对于高社会经济地位组而言是 −0.239-0.239。生活事件的影响看起来对于低社会经济地位组更剧烈一些，但是两组生活事件效应之间的区别是不显著的（p-value 为 0.44255）。

比较累计概率的解释。为描述定量变量的效应，我们比较这些变量的分位数处的累积概率；为描述定性变量的效应，我们比较这些变量的不同类别处的累积概率。对于定量变量通过令它们取均值来进行控制；对于属性变量通过固定类别来进行控制。如果有几个类别，就令它们取它们的指示变量的平均值。现在我们以上数据集心理健康为“健康”的结果作为栗子。首先考虑SESSES效应，我们需要控制生活事件数，于是我们取SES的均值4.275。> (x=mean(impairment$LifeEvents)) [1] 4.275 > round(predict(impairment_m,data.frame(LifeEvents=rep(x,2),SES=c(1,0)),type="response")[,1],4) 1 2 0.3696 0.1618

对于高社会经济地位者， P^(Y≤1)=P^(Y=1)=0.37\hat{P}(Y\leq1)=\hat{P}(Y=1)=0.37，对于低社会经济地位者P^(Y≤1)=P^(Y=1)=0.16\hat{P}(Y\leq1)=\hat{P}(Y=1)=0.16 。

再考虑生活事件效应。生活事件的上下四分位数分别为 2,6.252,6.25 ：

> quantile(impairment$LifeEvents) 0% 25% 50% 75% 100% 0.00 2.00 4.00 6.25 9.00 > round(predict(impairment_m,data.frame(LifeEvents=c(2,2,6.5,6.5),SES=c(1,0,1,0)),type="response")[,1],4) 1 2 3 4 0.5478 0.2850 0.2239 0.0867

对于高 SESSES ，在两个分位数处，健康的概率从0.55变为0.22,；对于低 SESSES，在两个分位数处，健康的概率从0.29变为0.09。这印证了我们之前的观点——给定社会经济地位，生活事件数的增加会对心理健康有负向效应。而且，给定生活事件数，社会经济地位低者比社会经济地位高者更容易受到心理伤害。