Формула R в survfit

У меня проблемы с формулами, средами и survfit().

Все работает нормально для lm(), но не работает для survfit().

Общая постановка проблемы:

Я подгоняю ряд формул к некоторым данным. Итак, я вызываю функцию моделирования с формулой, переданной в качестве переменной. Позже я хочу работать с формулой из подогнанного объекта.

(С моей наивной точки зрения проблема заключается в том, что survfit не записывает окружающую среду.)

Подробный пример

Ожидаемое поведение, как показано в lm():

library("plyr")

preds <- c("wt", "qsec")

f <- function() {
  lm(mpg ~ wt, data = mtcars)
}

fits <- alply(preds, 1, function(pred)
{
  modform <- reformulate(pred, response = "mpg")

  lm(modform, data = mtcars)
})

fits[[1]]$call$formula
##modform
formula(fits[[1]])
## mpg ~ wt
## <environment: 0x1419d1a0>

Несмотря на то, что fits[[1]]$call$formula разрешается в modform, я все еще могу получить исходную формулу с помощью formula(fits[[1]]).

Но у survfit() ничего не получается:

library("plyr")
library("survival")

preds <- c("resid.ds", "rx", "ecog.ps")

fits <- 
  alply(preds, 1, function(pred)
  {
    modform <- paste("Surv(futime, fustat)", pred, sep = " ~ ")
    modform <- as.formula(modform)
    print(modform)

    fit <- survfit(modform, data = ovarian)
  })

fits[[1]]$call$formula
## modform
formula(fits[[1]])
## Error in eval(expr, envir, enclos) : object 'modform' not found

Здесь (и в отличие от lm-фитов) formula(fits[[1]]) не работает!

Итак, мой конкретный вопрос: как я могу вернуть формулу, используемую для соответствия fits[[1]]?


r formula survival-analysis
person Andreas    schedule 16.10.2016    source источник
comment
fit <- eval(substitute(survfit(modform, data = ovarian), list(modform = modform))) или просто добавьте формулу обратно в объект survfit fit$call$formula <- modform, затем верните fit   -  person rawr    schedule 17.10.2016
comment
Спасибо! Это работает. Я так понимаю, нет возможности восстановить формулу из fits[[1]], поскольку она создается в примере кода вопроса? В таком случае: если вы напишете свой комментарий в качестве ответа, я приму это.   -  person Andreas    schedule 17.10.2016

Ответы (1)


arrow_upward
3
arrow_downward

Проблема в том, что когда x$formula равно NULL, для объекта lm существует резервный план получения формулы; этого не существует для survfit объектов

library("plyr")
library("survival")

preds <- c("wt", "qsec")
f <- function() lm(mpg ~ wt, data = mtcars)

fits <- alply(preds, 1, function(pred) {
  modform <- reformulate(pred, response = "mpg")
  lm(modform, data = mtcars)
})

fits[[1]]$formula
# NULL

Формула может быть извлечена с помощью formula(fits[[1]]), который использует общий formula. Метод lm S3 для formula таков:

stats:::formula.lm

# function (x, ...) 
# {
#   form <- x$formula
#   if (!is.null(form)) {
#     form <- formula(x$terms)
#     environment(form) <- environment(x$formula)
#     form
#   }
#   else formula(x$terms)
# }

Поэтому, когда fits[[1]]$formula возвращает NULL, forumla.lm ищет атрибут terms в объекте и таким образом находит формулу.

fits[[1]]$terms

У объектов survfit нет x$formula или x$terms, поэтому formula(x) выдает ошибку

preds <- c("resid.ds", "rx", "ecog.ps")
fits <-  alply(preds, 1, function(pred) {
    modform <- paste("Surv(futime, fustat)", pred, sep = " ~ ")
    modform <- as.formula(modform)
    fit <- survfit(modform, data = ovarian)
  })

fits[[1]]$formula
# NULL

formula(fits[[1]]) ## error

formula(fits[[1]]$terms)
# list()

Вы можете исправить это, вставив формулу в вызов и оценив ее

modform <- as.formula(paste("Surv(futime, fustat)", 'rx', sep = " ~ "))
substitute(survfit(modform, data = ovarian), list(modform = modform))
# survfit(Surv(futime, fustat) ~ rx, data = ovarian)

eval(substitute(survfit(modform, data = ovarian), list(modform = modform)))

# Surv(futime, fustat) ~ rx

# Call: survfit(formula = Surv(futime, fustat) ~ rx, data = ovarian)
# 
#      n events median 0.95LCL 0.95UCL
# rx=1 13      7    638     268      NA
# rx=2 13      5     NA     475      NA

Или вручную поместив формулу в x$call$formula

fit <- survfit(modform, data = ovarian)
fit$call$formula
# modform
fit$call$formula <- modform
fit$call$formula
# Surv(futime, fustat) ~ rx

fit
# Call: survfit(formula = Surv(futime, fustat) ~ rx, data = ovarian)
# 
#      n events median 0.95LCL 0.95UCL
# rx=1 13      7    638     268      NA
# rx=2 13      5     NA     475      NA
person rawr    schedule 17.10.2016
comment
Действительно хороший ответ. Очень признателен. - person Andreas; 18.10.2016