用R将对数似然写成函数(θ是什么？)

Question

我从我的模型中得到了下面的日志可能性，我正试图将它写成R中的函数。

​

​

我的问题是我不知道如何用函数来写θ。我在这方面做了几次尝试，如下所示，任何关于这些建议是否接近正确的建议都是值得赞赏的。

函数与θ一起编写为θ

#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)
}

我下一次尝试用我的数据集中的Xi替换thetas，这里是(登革热$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()，如下所示)
• 您在第一个术语中将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)的起始参数有点危险--了解参数的含义并猜测生物上合理的值可能更有意义--但这似乎是可以的。

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

由于我们没有限制θ是正的，所以我们会收到一些关于取负数日志的警告，但是这些可能是无害的(例如，参见here)。

替代方案

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

MASS::glm.nb自动为您设置一切，您只需指定预测变量(它使用对数链接并添加一个截距，因此指定~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))

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