R语言使用二进制回归将序数数据建模为多元GLM

 此结果的预测变量仅以一种方式对其产生影响，因此 为每个预测变量获得一个回归系数。但是该模型有几个截距，它们代表将变量切分以创建观察到的分类表现的点。

就像在普通回归模型中一样，每个预测变量都会以一种方式影响结果，这就是比例赔率假设或约束。或者，可以让每个预测变量在每个切入点对结果产生不同的影响。

如何使用单变量GLM软件对此建模？UCLA idre页面上有关于多元随机系数模型的文章。在这里很重要，因为他们使用nlme（单变量线性混合模型软件）对多元结果进行建模。基本思想是将数据堆叠起来，使其成为一种重复测量，但是找到一种向软件发出信号的信号，即结果是不同的，从而对预测变量要求不同的截距和斜率。

因此，我们要做的是将数据从宽转换为长，将其建模为常规二项式，但是我们需要告诉模型为每个级别估计不同的截距。为此，我使用具有unstructured工作相关性结构的通用估计方程（GEE）。3

示范

library(ordinal) # For ordinal regression to check our results
library(geepack) # For GEE with binary data

数据集。

soup <- ordinal::soup
soup$ID <- 1:nrow(soup) # Create a person ID variable
str(soup)

'data.frame':	1847 obs. of  13 variables:
 $ RESP    : Factor w/ 185 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ PROD    : Factor w/ 2 levels "Ref","Test": 1 2 1 2 1 2 2 2 2 1 ...
 $ PRODID  : Factor w/ 6 levels "1","2","3","4",..: 1 2 1 3 1 6 2 4 5 1 ...
 $ SURENESS: Ord.factor w/ 6 levels "1"<"2"<"3"<"4"<..: 6 5 5 6 5 5 2 5 5 2 ...
 $ DAY     : Factor w/ 2 levels "1","2": 1 1 1 1 2 2 2 2 2 2 ...
 $ SOUPTYPE: Factor w/ 3 levels "Self-made","Canned",..: 2 2 2 2 2 2 2 2 2 2 ...
 $ SOUPFREQ: Factor w/ 3 levels ">1/week","1-4/month",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ COLD    : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...
 $ EASY    : Factor w/ 10 levels "1","2","3","4",..: 7 7 7 7 7 7 7 7 7 7 ...
 $ GENDER  : Factor w/ 2 levels "Male","Female": 2 2 2 2 2 2 2 2 2 2 ...
 $ AGEGROUP: Factor w/ 4 levels "18-30","31-40",..: 4 4 4 4 4 4 4 4 4 4 ...
 $ LOCATION: Factor w/ 3 levels "Region 1","Region 2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ ID      : int  1 2 3 4 5 6 7 8 9 10 ...

我使用SURENESS变量。它有6个级别。使用DAY和GENDER变量对其进行建模。

# Select variables to work with
soup <- dplyr::select(soup, ID, SURENESS, DAY, GENDER)
# I like dummy variables with recognizable names
soup$girl <- ifelse(soup$GENDER == "Female", 1, 0) # Make male reference group
soup$day2 <- ifelse(soup$DAY == "2", 1, 0) # Make day 1 reference group

下一步是将顺序结果转换为代表每个阈值的5个结果

完成此操作后，我们准备对这5个新的结果变量进行转换。

head(soup.long) # Let's look at the data

     ID SURENESS DAY GENDER girl day2 SURE VAL SURE.f
1     1        6   1 Female    1    0    2   1      2
1848  1        6   1 Female    1    0    3   1      3
3695  1        6   1 Female    1    0    4   1      4
5542  1        6   1 Female    1    0    5   1      5
7389  1        6   1 Female    1    0    6   1      6
2     2        5   1 Female    1    0    2   1      2

让我们看看没有选择最高响应类别的人：

     ID SURENESS DAY GENDER girl day2 SURE VAL SURE.f
22   22        4   1 Female    1    0    2   1      2
1869 22        4   1 Female    1    0    3   1      3
3716 22        4   1 Female    1    0    4   1      4
5563 22        4   1 Female    1    0    5   0      5
7410 22        4   1 Female    1    0    6   0      6

这个人选择了SURENESSVAL中的a 。她的前三个分数是1，她的最后两个分数是0，因为4小于4-5阈值和5-6阈值。

下一步是为阈值创建虚拟变量。这些变量将用于表示模型中的截距。

请注意，我将虚拟变量乘以-1。在序数回归中，这样做使解释更容易。总之，它确保正系数增加了从较低类别（例如3）移至较高类别（4）或对较高响应类别做出响应的几率。

现在，我们准备运行模型。我们使用GEE。相关结构为unstructured。

接下来，我使用标准序数回归估算模型：

让我们比较系数和标准误差：

        Estimate Estimate.1 Std.err Std. Error     Wald  z value Pr(>|W|) Pr(>|z|)
SURE.f2 -2.13244   -2.13155 0.10454    0.10450 416.0946 -20.3971   0.0000   0.0000
SURE.f3 -1.19345   -1.19259 0.09142    0.09232 170.4284 -12.9179   0.0000   0.0000
SURE.f4 -0.89164   -0.89079 0.08979    0.09011  98.5995  -9.8857   0.0000   0.0000
SURE.f5 -0.65782   -0.65697 0.08945    0.08898  54.0791  -7.3833   0.0000   0.0000
SURE.f6 -0.04558   -0.04477 0.08801    0.08789   0.2682  -0.5093   0.6046   0.6105
girl    -0.04932   -0.04917 0.09036    0.09074   0.2980  -0.5419   0.5851   0.5879
day2    -0.26172   -0.26037 0.08584    0.08579   9.2954  -3.0351   0.0023   0.0024

可以看到结果非常接近。

但是，使用估计glm()不能建立一个人的结果之间的依存关系的估计会产生不同的结果。

        Estimate Std. Error z value Pr(>|z|)
SURE.f2 -2.15144    0.08255 -26.062   0.0000
SURE.f3 -1.21271    0.06736 -18.004   0.0000
SURE.f4 -0.91149    0.06472 -14.084   0.0000
SURE.f5 -0.67782    0.06327 -10.713   0.0000
SURE.f6 -0.06523    0.06178  -1.056   0.2911
girl    -0.07326    0.04961  -1.477   0.1398
day2    -0.26898    0.04653  -5.780   0.0000

估计值和标准误均不足。

我们可以轻松地放宽pom.bin模型中的比例赔率约束。让我们通过放宽对预测变量的约束来运行某些人所说的偏比例赔率模型day2。我们通过估计阈值虚拟变量和day2预测变量之间的相互作用来做到这一点。

我还使用名义参数运行了相同的模型进行比较day2。

             Estimate Estimate.1 Std.err Std. Error     Wald  z value Pr(>|W|) Pr(>|z|)
SURE.f2      -2.02982   -2.03106 0.11800    0.11834 295.8986 -17.1630  0.00000  0.00000
SURE.f3      -1.22087   -1.22213 0.09829    0.09857 154.2801 -12.3980  0.00000  0.00000
SURE.f4      -0.92773   -0.92899 0.09458    0.09443  96.2112  -9.8375  0.00000  0.00000
SURE.f5      -0.65744   -0.65870 0.09246    0.09188  50.5554  -7.1693  0.00000  0.00000
SURE.f6      -0.04733   -0.04859 0.08955    0.08965   0.2793  -0.5420  0.59714  0.58784
SURE.f2:day2  0.07359    0.07360 0.14148    0.14155   0.2705   0.5199  0.60298  0.60312
SURE.f3:day2  0.31691    0.31697 0.10607    0.10613   8.9270   2.9867  0.00281  0.00282
SURE.f4:day2  0.33301    0.33308 0.09970    0.09973  11.1551   3.3398  0.00084  0.00084
SURE.f5:day2  0.26330    0.26339 0.09618    0.09616   7.4938   2.7391  0.00619  0.00616
SURE.f6:day2  0.26741    0.26748 0.09347    0.09345   8.1842   2.8622  0.00423  0.00421
girl         -0.04809   -0.04994 0.09048    0.09077   0.2825  -0.5502  0.59507  0.58221

结果是可比较的。

现在，我们可以将比例比例赔率二进制模型与比例赔率二进制模型进行比较，以测试day2变量的约束条件。geepack允许anova()对两种模型进行Wald测试：

Analysis of 'Wald statistic' Table

Model 1 VAL ~ 0 + SURE.f2 + SURE.f3 + SURE.f4 + SURE.f5 + SURE.f6 + girl + SURE.f2:day2 + SURE.f3:day2 + SURE.f4:day2 + SURE.f5:day2 + SURE.f6:day2
Model 2 VAL ~ 0 + SURE.f2 + SURE.f3 + SURE.f4 + SURE.f5 + SURE.f6 + girl + day2
  Df   X2 P(>|Chi|)
1  4 6.94      0.14

两种模型之间的差异在统计上均不显着，表明day2变量的比例约束已足够。

我们可以使用或使用函数ordinal进行比较pom.ord和npom.ord建模anova()，从而进行相同的测试nomimal_test()。两者都是似然比检验，比上述GEE的Wald检验更充分。

Likelihood ratio tests of cumulative link models:

         formula:               nominal: link: threshold:
pom.ord  SURENESS ~ girl + day2 ~1       logit flexible  
npom.ord SURENESS ~ girl        ~day2    logit flexible  

         no.par  AIC logLik LR.stat df Pr(>Chisq)
pom.ord       7 5554  -2770                      
npom.ord     11 5555  -2766    6.91  4       0.14

nominal_test(pom.ord)

Tests of nominal effects

formula: SURENESS ~ girl + day2
       Df logLik  AIC  LRT Pr(>Chi)  
<none>     -2770 5554                
girl    4  -2766 5554 8.02    0.091 .
day2    4  -2766 5555 6.91    0.141  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这两个测试收敛到相同的结果，并且在比较GEE模型的Wald测试中也给出了相同的p值。然而，Wald- χ 2χ2 测试统计数据略高。

---

完成此操作后，使用序数数据包当然要容易得多。但是，将模型视为二进制可能会有一些好处，但是所有这些都是出于好奇而非必要。由于某种原因，我仍未弄清楚，当一个人尝试使用fitted()函数从模型中获得预测的概率时，它仅返回一组拟合的概率。理想情况下，它应该为每个阈值返回拟合概率。使用geepack，可以直接获得每个级别的预测概率。但是，这种优势是微不足道的。

---

而且，如果熟悉最大似然估计，则可以简单地对似然函数进行编程。

上面的例子在比例赔率情况下的语法为：

coef(summary(res))
      Estimate Std. Error
a1 -2.13155603 0.10450286
a2 -1.19259266 0.09232077
a3 -0.89079068 0.09010891
a4 -0.65697671 0.08898063
a5 -0.04477565 0.08788869
bg -0.04917604 0.09073602
bd -0.26037369 0.08578617

coef(summary(pom.ord))
        Estimate Std. Error     z value     Pr(>|z|)
1|2  -2.13155281 0.10450291 -20.3970663 1.775532e-92
2|3  -1.19259171 0.09232091 -12.9178937 3.567748e-38
3|4  -0.89078590 0.09010896  -9.8856524 4.804418e-23
4|5  -0.65697465 0.08898068  -7.3833401 1.543671e-13
5|6  -0.04476553 0.08788871  -0.5093434 6.105115e-01
girl -0.04917245 0.09073601  -0.5419287 5.878676e-01
day2 -0.26037360 0.08578617  -3.0351465 2.404188e-03

结果非常相似，对于比较模型的更确定的方式，我们总是可以比较对数似然：

logLik(res)
'log Lik.' -2769.784 (df=7)

logLik(pom.ord)
'log Lik.' -2769.784 (df=7)

---

1. Agresti，A。（2013）。分类数据分析。Wiley-Interscience。 ↩