1

我试图获得用于生存分析的 Gumbel 分布的对数似然的最大似然估计量（我这样说是为了让你不会被对数似然函数迷惑，我认为它是正确的）。为了做到这一点，我必须通过使用 optim 函数来最大化负对数似然，我试图这样做，但控制台在 fn(par, ...) 中给了我一个错误：缺少参数“b”，没有默认值。

我也尝试以与此链接的答案类似的方式执行此操作：在约束条件下使用两个参数求解最大似然性，但控制台游戏我如下： optim(c(1, 1), function(x) 中的错误) log_lhood(x[1], x[2], d = lung$status, : 优化中的目标函数计算为长度 0 而不是 1。

log_lhood <- function(m,b,d,t){                           
  sum<-0
  for (i in 1:length(lung)){ 
    if (d[i] == 1){
      sum<- sum - log(1-exp(-exp(-(t[i]-m)/b)))
    } else {
      sum<- sum - log((1/b)*exp(-(t[i]-m/b+exp(-(t[i]-m/b)))))
    }
  }
}
#a,b parameter optimization

optim(c(0,0), fn = log_lhood, d = lung$status, t = KM_fit$time) #1st way

optim(c(1, 1), function(x) log_lhood(x[1], x[2],d=lung$status,t=KM_fit_test$time)) #2nd way as in the link
4

1 回答 1

2

这里有几个问题。函数的第一个参数应该是参数向量。你nrow(lung)也不需要length(lung)，最好使用length(d). 此外，您不应该在这里使用循环，它的使用效率非常低ifelse()（在 R 中，我们总是尝试对所有内容进行矢量化）。您还需要检查是否可以为所有参数值计算对数似然度（例如 b=0）。你也忘了return(sum)。这sum也是一个你不应该掩盖的有用功能。

这运行。

library(reprex)

lung <- data.frame(status=c(0,0,1,1))
KM_fit <- data.frame(time=c(0,1,2,3))

log_lhood <- function(x,d,t){
    m <- x[1]
    b <- x[2]
    sum <-0
    for (i in 1:nrow(lung)){
        if (d[i] == 1){
            sum <- sum - log(1-exp(-exp(-(t[i]-m)/b)))
        } else {
            sum <- sum - log((1/b)*exp(-(t[i]-m/b+exp(-(t[i]-m/b)))))
        }
    }
    return(sum)
}
#a,b parameter optimization

optim(par=c(0,1), fn = log_lhood, d = lung$status, t = KM_fit$time)
$par
[1] 1.661373 1.811780

$value
[1] 5.318068

$counts
function gradient 
      63       NA 

$convergence
[1] 0

$message
NULL

你可以像这样重写你的函数。

log_lhood <- function(x,d,t){
    m <- x[1]
    b <- x[2]
    s <- ifelse(d==1,
                  -log(1-exp(-exp(-(t-m)/b))),
                  -log((1/b)*exp(-(t-m/b+exp(-(t-m/b)))))
                  )
    return(sum(s, na.rm=TRUE))
}
于 2019-05-06T01:47:17.570 回答