原文鏈接：http://tecdat.cn/?p=10204
用於分析序數數據的最常見模型是 邏輯模型 。本質上，您將結果視為連續潛在變量的分類表現。此結果的預測變量僅以一種方式對其產生影響，因此 為每個預測變量獲得一個回歸系數。但是該模型有幾個截距，它們代表將變量切分以創建觀察到的分類表現的點。
就像在普通回歸模型中一樣，每個預測變量都會以一種方式影響結果，這就是比例賠率假設或約束。或者，可以讓每個預測變量在每個切入點對結果產生不同的影響。
如何使用單變量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)