在 R 中将对数似然写成函数(什么是 theta)

基本问题是，必须将beta(问题一个分量的截距和斜率)和theta作为单个参数向量的一部分进行传递.我修复了括号位置的其他问题，并稍微重新组织了表达式.

代码中还有几个更重要的错误.

• 第一项为not a fraction；这是一个二项式系数.(即，您应该使用lchoose()，如下所示)
• 您在第一学期将a+1改为a-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)的起始参数有点危险-了解meaning个参数的一些信息并猜测生物学上合理的值会更有意义-但似乎没有问题.

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

由于我们没有将θ限制为正值，因此我们得到了一些关于记录负数的警告，但这些警告可能是无害的(例如，请参阅here)

Select

R有很多内置的机器来拟合负二项模型(我应该知道你在做什么！)

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

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

bbmle自动化程度稍低，但更灵活(这里我在对数刻度上拟合theta，以避免try 任何负值)

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

所有这三种方法(修正的负对数似然函数+optim、MASS::glm.nb、bbmle::mle2)都给出了相同的结果.