模擬有序 logit 模型的數據

我目前正在從事一個項目，其中我的因變量以 1-5 的李克特量表的形式可用。但在我在真實數據上運行我的模型之前，我想在人工數據上對其進行測試，但我正在努力如何模擬這種數據的數據。

我想本著 Polson、Scott 和 Windle (2013) 以及他們在 R 中的包“BayesLogit”的精神建立一個貝葉斯有序 logit 模型（稍後我想引入隨機效應，但這並不重要）。

我希望任何人都可以提供幫助！

假設您有一個有序的分類因變量，它是一個從 1 到 5 的李克特量表（我將其編碼為y），以及兩個正態分佈和均勻分佈的自變量（分別為x1和）。x2

基本的有序邏輯回歸所做的是預測對數賠率，其中是因變量有多少個級別（在本例中為 4）：

首先，屬於 2、3、4 和 5 類的概率之和與屬於 1 類的概率的對數。然後，類別 3、4 和 5 超過 1 和 2 的總和，等等。這些看起來像：

,,,

為方便起見，我們將這些對數賠率分別稱為 1、2、3 和 4。

那麼我們如何得到這四個對數賠率預測呢？他們每個人都有自己的方程式：

您會注意到，這些方程的唯一不同之處在於它們具有不同的截距。這是比例優勢或平行回歸假設。現在，我們可以基於這個模型來模擬數據。既然是訂購的，（例如，因為 a、b、c 和 d 發生的機率不能小於 a、b 和 c 發生的機率）。

首先，讓我們製作x1and x2：

set.seed(1839)
n <- 150
x1 <- rnorm(n)
x2 <- runif(n)

現在，讓我們設置預測對數機率的總體參數：

b01 <- 1
b02 <- 0.05
b03 <- -0.05
b04 <- -1
b1 <- 0.8
b2 <- 0.5

現在，我們可以根據上面的等式預測每個對數賠率：

logodds1 <- b01 + b1 * x1 + b2 * x2
logodds2 <- b02 + b1 * x1 + b2 * x2
logodds3 <- b03 + b1 * x1 + b2 * x2
logodds4 <- b04 + b1 * x1 + b2 * x2

現在我們可以使用得到每個人的分子的概率：

inv_logit <- function(logit) exp(logit) / (1 + exp(logit))
prob_2to5 <- inv_logit(logodds1)
prob_3to5 <- inv_logit(logodds2)
prob_4to5 <- inv_logit(logodds3)
prob_5 <- inv_logit(logodds4)

現在，我們可以做減法來得到每個類別本身的概率：

prob_1 <- 1 - prob_2to5
prob_2 <- prob_2to5 - prob_3to5
prob_3 <- prob_3to5 - prob_4to5
prob_4 <- prob_4to5 - prob_5

只是為了確定，所有概率都高於零嗎？

> table(c(prob_1, prob_2, prob_3, prob_4, prob_5) > 0)

TRUE 
750 

涼爽的。現在，我們可以使用這些概率從數字 1 到 5（感興趣的類別）進行採樣，以獲得我們的因變量y：

y <- c()
for (i in 1:n) {
 y[i] <- sample(
   x = c(1:5), 
   size = 1, 
   prob = c(prob_1[i], prob_2[i], prob_3[i], prob_4[i], prob_5[i])
 )
}

我本可以在apply那裡發表聲明，但我很懶惰。

現在，我們可以將它們放入 adata.frame並運行序數回歸：

> dat <- data.frame(x1, x2, y = factor(y))
> 
> library(ordinal)
> summary(clm(y ~ x1 + x2, data = dat))
formula: y ~ x1 + x2
data:    dat

link  threshold nobs logLik  AIC    niter max.grad cond.H 
logit flexible  150  -204.23 420.47 8(1)  2.74e-07 2.2e+02

Coefficients:
  Estimate Std. Error z value Pr(>|z|)    
x1   0.8933     0.1786   5.003 5.65e-07 ***
x2   0.5810     0.5595   1.038    0.299    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
   Estimate Std. Error z value
1|2 -1.25461    0.35426  -3.541
2|3  0.02907    0.33321   0.087
3|4  0.25541    0.33261   0.768
4|5  0.91630    0.34070   2.689

請注意，我們恢復了參數 OK，但不完全正確。為什麼？該sample函數有隨機誤差。您可以將其寫入一個函數，執行多次，然後查看您恢復參數的次數等。

該代碼並不完美，但它應該可以幫助您入門。出於教學目的，它也更加冗長。

由於您是在貝葉斯框架中進行的，因此我非常感謝您回答我的 CrossValidated關於為序數回歸選擇半信息先驗的問題。

引用自：https://stats.stackexchange.com/questions/321770