如何在 R 中编码，这个对数似然函数

Question

给出了有序logit模型；

考虑 Logit 模型的以下对数似然函数：

还给出 k = {0,1,2,3,4}，其中对于 k = 0 和 k = 4，我们的 alpha 分别为 -Inf 和 +Inf。现在我的问题是如何编写一个优化对数似然函数的函数，即找到 B = B1 的 MLE 和 k = {1,2,3} 的 alpha。

这是我的尝试；

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)

只有这样不行，我该如何解决这个问题？

Answer 1

我认为第一个问题是您在内循环中

return(ll)

，这将导致第一个记录的第一个值——您的过程不会迭代，直到它到达

return(sum(ll))

。第二个问题是，观察结果应该只对它们实际所属类别的可能性做出贡献，你不做这种区分。最后，

optim

最小化，所以你不应该返回你想要最大化的值。

为了测试实现，我们首先获取一些示例数据（来自here）：

## Example data
dat <- foreign::read.dta("https://stats.idre.ucla.edu/stat/data/ologit.dta")

Y <- as.integer(cut(dat$gpa, 4))  ## Create 4 quartiles
X <- dat$public                   ## Dummy predictor, could be slope too

## This is the result we'll be reproducing
MASS::polr(as.factor(Y) ~ X)

#> Coefficients:
#>        X 
#> 1.178314 
#> 
#> Intercepts:
#>         1|2         2|3         3|4 
#> -2.43862931 -0.06816316  2.05272249 

现在，这是一个可以为您提供大致相同的参数估计值的实现：

## Logit link inversion
expit <- function(x) 1/(1+exp(-x))

## -log-likelihood function
ll <- function(params, X, Y) {
  alpha <- params[-1]
  beta <- params[1]
  
  ll <- unlist(mapply(\(x, y) {
    eta <- alpha - x*beta
    mu <- expit(eta)
    ## This line calculates expit(delta) - expit(delta_-1) in probability scale
    probs <- pmin(pmax(diff(c(0, mu, 1)), .Machine$double.eps), 1)
    ## Only use the probability of the observed category!
    log(probs[y])
  }, x=X, y=Y))
  ## Return inverse so value is maximized
  -sum(ll)
}

guess <- c(1, -2, 0, 2)
optim(guess, ll, X=X, Y=Y)
#> $par
#> [1]  1.17991554 -2.43853589 -0.06775819  2.05298068

一些注意事项：

• 我计算了概率尺度的可能性，然后记录下来，而不是直接使用您给出的公式（但简单的代数会表明它们是相同的）。我没有将对数赔率设置为（负）无穷大，而是将概率设置为 0/1。
• 您需要提供合适的起始值才能使其发挥作用。特别是，截距应该已经是不同的并且顺序正确（总是根据您想要对类别进行排序的方式增加或减少），从所有它们相等开始会导致问题。
• 这依赖于

  Y

  是从 1 开始的数字 + 顺序，并且仅处理单个数字

  X

  预测器，正确的拟合例程应该允许更大的灵活性（例如，在内部将

  Y

  的级别转换为此有序序列） .