如何用R编写有序logit模型对数似然函数的优化代码？

问题背景 #

给定有序logit模型及对应对数似然函数，已知k的取值集合为{0,1,2,3,4}，其中k=0时α为-∞，k=4时α为+∞。需编写R函数优化该对数似然函数，求解β=β₁和k={1,2,3}对应的α的极大似然估计（MLE）。

以下是尝试的代码，但无法运行，需修复：

J <- length(Poss)

log_likelihood <- function(params, Poss = Poss){
  beta <- params[1]
  alpha <- c(-Inf, params[2:4], +Inf)
  
  for (k in 2:4) {
    for (j in J) {
      delta[k] <- alpha[k] - Poss[j] * beta
      ll <- delta[k] - log(1 + exp(delta[k]) - delta[k] + log(1 + exp(delta[k])))
      return(ll)
    }
  }
  return(sum(ll))
}

guess_optim <- c(0,0,0,0) #for beta1 and alpha1,2,3

optim(guess_optim, log_likelihood, Poss = Poss)

代码问题分析与修复 #

核心问题点 #

• 循环逻辑错误：内层循环for (j in J)仅会执行一次（J是长度值而非序列），应改为for (j in 1:J)遍历所有样本
• 对数似然公式错误：原公式完全不符合有序logit的似然结构，混淆了累积概率的计算逻辑
• 变量未初始化：delta和ll未提前声明，会导致赋值报错
• 提前返回：循环内部直接return(ll)会让函数提前退出，无法计算所有样本的似然和
• 优化方向错误：optim默认极小化，但我们需要极大化对数似然，需调整方向

修复后的代码 #

注意：需补充因变量y（对应每个样本的k类别取值0-4），这是原代码缺失的关键输入：

# 示例数据（替换为你的实际数据）
set.seed(123)
n <- 100
Poss <- rnorm(n)
y <- sample(0:4, n, replace = TRUE)

log_likelihood <- function(params, x, y) {
  beta <- params[1]
  alpha <- c(-Inf, params[2:4], Inf)
  
  ll_total <- 0
  n <- length(y)
  
  for (i in 1:n) {
    k_i <- y[i]
    x_i <- x[i]
    
    # 计算logistic累积概率
    F_k_prev <- plogis(alpha[k_i] - x_i * beta)
    F_k_curr <- plogis(alpha[k_i + 1] - x_i * beta)
    
    # 计算当前样本的对数似然
    ll_i <- log(F_k_curr - F_k_prev)
    ll_total <- ll_total + ll_i
  }
  
  # 返回负对数似然适配optim的极小化逻辑
  return(-ll_total)
}

# 初始值：beta + alpha1、alpha2、alpha3
initial_guess <- c(0, -1, 0, 1)

# 执行优化，添加alpha单调性约束（有序logit要求α1<α2<α3）
optim_result <- optim(initial_guess, log_likelihood, x = Poss, y = y, 
                      method = "L-BFGS-B", 
                      lower = c(-Inf, -Inf, -Inf, -Inf),
                      upper = c(Inf, Inf, Inf, Inf),
                      control = list(fnscale = 1))

# 查看估计结果
optim_result

关键说明 #

• 单调性约束：有序logit要求阈值α满足α₁<α₂<α₃，若需严格约束，可在函数内部判断，若α不满足递增则返回极大值
• plogis函数：R中plogis()等价于logistic累积分布函数1/(1+exp(-z))，可直接简化累积概率计算
• 因变量y：必须输入每个样本对应的类别k值，否则无法计算对应区间的似然

内容的提问来源于stack exchange，提问作者Kemit4