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

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

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

r max analysis survival log-likelihood
person nnisgia    schedule 05.05.2019    source источник
comment
Это почти можно проверить, пожалуйста, укажите минимальную сумму, чтобы кто-то мог запустить и испытать вашу ошибку. lung из пакета?   -  person Evan Friedland    schedule 06.05.2019
comment
это из пакета выживания. ты прав прости.   -  person nnisgia    schedule 06.05.2019

Ответы (1)


arrow_upward
2
arrow_downward

Здесь есть несколько проблем. Предполагается, что первым аргументом функции является вектор параметров. Также вам нужно 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))
}
person Simon Woodward    schedule 06.05.2019
comment
Неспособность вернуть sum в этом случае — это нормально (хотя в целом это не очень хорошая практика), потому что это всегда последняя переменная, назначенная внутри функции. - person Frank; 06.05.2019
comment
Большое спасибо, теперь все работает нормально, я учту ваши комментарии. еще раз большое спасибо! - person nnisgia; 06.05.2019