在 R 中将对数似然写为函数（什么是 theta？）

Question

我的模型中有以下对数似然，我试图将其作为 R 中的函数编写。

我的问题是因为我不知道如何根据函数编写 theta。我已经进行了几次尝试，如下所示，任何 tips/advice 如果这些接近正确，我们将不胜感激。

theta 写成 theta 的函数

#my likelihood function
mylikelihood = function(beta) {

#log-likelihood
result = sum(log(dengue$cases + theta + 1 / dengue$cases)) +
  sum(theta*log(theta / theta + exp(beta[1]+beta[2]*dengue$time))) + 
  sum(theta * log(exp(beta[1]+beta[2]*dengue$time / dengue$cases + exp(beta[1]+beta[2]*dengue$time))))

#return negative log-likelihood
return(-result)
}

我的下一次尝试将 thetas 替换为我数据集中的 Xi，这里是 (dengue$time)

#my likelihood function attempt 2
mylikelihood = function(beta) {

#log-likelihood
result = sum((log(dengue$Cases + dengue$Time + 1 / dengue$Cases))) +
sum(dengue$Time*log(dengue$time / dengue$Time + exp(beta[1]+beta[2]*dengue$Time))) + 
sum(dengue$Cases * log(exp(beta[1]+beta[2]*dengue$Time / dengue$Cases + 
exp(beta[1]+beta[2]*dengue$Time))))

 #return negative log-likelihood
 return(-result)
 }

数据

   head(dengue)
  Cases Week Time
1   148   36    1
2   275   37    2
3   205   38    3
4   133   39    4
5   123   40    5
6   138   41    6

其中任何一个接近正确吗？如果不正确，我哪里错了？

更新了对数似然的来源；

型号；

具有均值 µ 和色散参数 θ 的负二项分布具有 pmf；

Answer 1

根本问题是您必须将 beta（问题的一个分量的截距和斜率）和 theta 作为单个参数向量的一部分传递。我修复了括号放置的其他问题，我稍微重新组织了表达式。

您的代码中有几个更重要的错误。

第一项不是分数；它是一个二项式系数。（即，您应该使用 lchoose()，如下所示）
您在第一学期将 +1 更改为 -1。

nll <- function(pars) {                                                                                      
    beta <- pars[1:2]                                                                                        
    theta <- pars[3]                                                                                         
                                                                                                             
    ##log-likelihood                                                                                         
    yi <- dengue$Cases                                                                                       
    xi <- dengue$Time                                                                                        
    ri <- exp(beta[1]+beta[2]*xi)                                                                            
    result <- sum(lchoose(yi + theta - 1,yi)) +                                                              
        sum(theta*log(theta / (theta + ri))) +                                                               
        sum(yi * log(ri/(theta+ri)))                                                                         
    ##return negative log-likelihood                                                                         
    return(-result)                                                                                          
}

读取数据

dengue <- read.table(row.names = 1, header = TRUE, text = "                                                  
 Cases Week Time                                                                                             
1   148   36    1                                                                                            
2   275   37    2                                                                                            
3   205   38    3                                                                                            
4   133   39    4                                                                                            
5   123   40    5                                                                                            
6   138   41    6                                                                                            
")

合身

猜测 (1,1,1) 的起始参数有点危险 - 了解参数的含义并猜测生物学上合理的值会更有意义- 不过好像还可以。

nll(c(1,1,1))                                                                                                
optim(par = c(1,1,1), nll)

由于我们没有将 theta 限制为正数，因此我们会收到一些关于对负数进行对数的警告，但这些可能是无害的（例如，参见）

备选方案

R 有很多 built-in 机制来拟合负二项式模型（我应该知道你在做什么！）

MASS::glm.nb 自动为您设置所有内容，您只需指定预测变量（它使用对数 link 并添加截距，因此指定 ~Time 将使均值等于 exp(beta0 + beta1*Time)).

library(MASS)
glm.nb(Cases ~ Time, data = dengue)

bbmle 的自动化程度稍低，但更灵活（这里我将 theta 安装在对数刻度上以避免尝试任何负值）

library(bbmle)
mle2(Cases ~ dnbinom(mu = exp(logmu), size = exp(logtheta)),
                     parameters = list(logmu ~ Time),
                     data = dengue,
                     start = list(logmu = 0, logtheta = 0))

所有这三种方法（修正负 log-likelihood 函数 + optim、MASS::glm.nb、bbmle::mle2）给出相同的结果。

在 R 中将对数似然写为函数（什么是 theta？）

Writing a log-likelihood as a function in R (what is theta?)

r

function

log-likelihood

读取数据

合身

备选方案