Есть ли способ, чтобы уже определенный аргумент мог отображаться как отсутствующий при использовании функции optim () в R?

Question

Я пытаюсь получить Оценки максимального правдоподобия логарифмического правдоподобия распределения Гамбеля для анализа выживания (я говорю, чтобы вы не отстали от логарифмической функции правдоподобия, я думаю, что это правильно).Чтобы сделать это, я должен максимизировать минус логарифмическую вероятность с помощью функции optim, я пытался это сделать, но консоль выдает мне сообщение об ошибке в fn (par, ...): отсутствует аргумент "b",без значения по умолчанию.

Я попытался сделать это так же, как в ответе на эту ссылку: Решить для максимальной вероятности с двумя параметрами при ограничениях , но консольная игра меняследующее: ошибка в optim (c (1, 1), function (x) log_lhood (x [1], x [2], d = lung $ status): целевая функция в optim оценивает длину 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

Answer 1

Здесь есть несколько проблем. Предполагается, что первым аргументом функции будет вектор параметров. Также вам нужно 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))
}