用purrr拟合多个tidy模型的生存模型

ecbunoof  于 2023-01-10  发布在  其他
关注(0)|答案(1)|浏览(130)

我正在尝试使用tidymodels,workflow和purr来拟合许多生存模型。我可以让这种方法适用于其他模型,例如线性回归,但不适用于生存。我已经将生存扩展加载到parsnip中。
以下代码用于
1.生成小数据集。
1.证明通常的cox-ph工作正常。
1.使用tidymodels和工作流运行线性回归,效果良好(models1)。
1.但是,models2和models3会导致错误消息:"fit_xy()"中的错误:!删失回归模型必须使用公式接口。
这个消息对我没有多大帮助,我想它与内嵌Surv()对象有关,但还没有看到如何用另一种方式构造它。
我知道防风草生存可能仍然是工作中的进展,但它是不清楚我是否已经实施,应该工作或没有。感谢任何帮助,谢谢。

library(tidyverse)
library(tidymodels)
library(survival)
library(censored)

set.seed(1973)
df <- tibble(id = seq(1:1000))
df <- df %>% mutate( survtime = floor(100*rgamma(n=1000, shape =1 , rate=1)) ,
                     fail = runif(n=1000) >0.33 ,
                      a1 = runif(n=1000) >0.1,
                     a2 = runif(n=1000) >0.5)
                     

head(df)

cox1 <- coxph(data = df, Surv(survtime, event=fail) ~a1)
summary(cox1)

cox2 <- coxph(data = df, Surv(survtime, event=fail) ~a2)
summary(cox2)

A <- c("a1", "a2")
models1 <- map(A,
               ~workflow() %>%
                 add_model(linear_reg()) %>%
                 add_formula(reformulate(.x, response = 'survtime')) %>%
                 fit(df)
)
models1

#not working
models2 <- map(A,
               ~workflow() %>%
                 add_model(proportional_hazards()) %>%
                 add_formula(reformulate(.x, response = 'Surv(survtime, event=fail)')) %>%
                 fit(df)
)

#not working
models3 <- map(A,
               ~workflow() %>%
                 add_model(proportional_hazards()) %>%
                 add_formula(as.formula(paste0('Surv(survtime, event=fail) ~ ', .x))) %>%
                 fit(df)
)
models3

编辑以将变量X、x1和x2更改为A、a1和a2,以便更清楚。错误消息相同:
"fit_xy()"中的错误:!删失回归模型必须使用公式接口。
但是,问题可能是R版本或包。错误发生在版本4.1.2与生存3.4 - 0 parsnip 1.0.1
但代码可在另一个安装了R 4.2.2 parsnip 1.0.3 survival 3.3 - 1的系统上运行
不幸的是,我在系统管理员的怜悯更新,所以我没有检查版本,不能很容易地排除故障。

agyaoht7

agyaoht71#

大写“X”而非小写“x”是问题所在。
你不是第一个这样做的人(我盯着它看了几分钟),所以不要难过。
然而......这个例子正是我们使用reprex软件包的原因。它会立即告诉你问题所在,因为它是从一个新的会话开始的(并且会向我们显示错误消息)
了解是成功的一半:-)

library(tidyverse)
library(tidymodels)
library(survival)
library(censored)

set.seed(1973)
df <- tibble(id = seq(1:1000))
df <- df %>% mutate( survtime = floor(100*rgamma(n=1000, shape =1 , rate=1)) ,
                     fail = runif(n=1000) >0.33 ,
                      x1 = runif(n=1000) >0.1,
                     x2 = runif(n=1000) >0.5)
                     

head(df)
#> # A tibble: 6 × 5
#>      id survtime fail  x1    x2   
#>   <int>    <dbl> <lgl> <lgl> <lgl>
#> 1     1       27 FALSE TRUE  FALSE
#> 2     2       77 FALSE TRUE  TRUE 
#> 3     3      180 FALSE TRUE  FALSE
#> 4     4      127 FALSE TRUE  TRUE 
#> 5     5       69 TRUE  TRUE  TRUE 
#> 6     6       20 FALSE TRUE  TRUE

cox1 <- coxph(data = df, Surv(survtime, event=fail) ~x1)
summary(cox1)
#> Call:
#> coxph(formula = Surv(survtime, event = fail) ~ x1, data = df)
#> 
#>   n= 1000, number of events= 692 
#> 
#>          coef exp(coef) se(coef)     z Pr(>|z|)  
#> x1TRUE 0.2854    1.3303   0.1250 2.284   0.0224 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>        exp(coef) exp(-coef) lower .95 upper .95
#> x1TRUE      1.33     0.7517     1.041     1.699
#> 
#> Concordance= 0.509  (se = 0.007 )
#> Likelihood ratio test= 5.62  on 1 df,   p=0.02
#> Wald test            = 5.22  on 1 df,   p=0.02
#> Score (logrank) test = 5.25  on 1 df,   p=0.02

cox2 <- coxph(data = df, Surv(survtime, event=fail) ~x2)
summary(cox2)
#> Call:
#> coxph(formula = Surv(survtime, event = fail) ~ x2, data = df)
#> 
#>   n= 1000, number of events= 692 
#> 
#>           coef exp(coef) se(coef)     z Pr(>|z|)  
#> x2TRUE 0.13210   1.14122  0.07655 1.726   0.0844 .
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>        exp(coef) exp(-coef) lower .95 upper .95
#> x2TRUE     1.141     0.8763    0.9822     1.326
#> 
#> Concordance= 0.513  (se = 0.011 )
#> Likelihood ratio test= 2.98  on 1 df,   p=0.08
#> Wald test            = 2.98  on 1 df,   p=0.08
#> Score (logrank) test = 2.98  on 1 df,   p=0.08

X <- c("x1", "x2")
models1 <- map(X,
               ~workflow() %>%
                 add_model(linear_reg()) %>%
                 add_formula(reformulate(.x, response = 'survtime')) %>%
                 fit(df)
)
models1
#> [[1]]
#> ══ Workflow [trained] ══════════════════════════════════════════════════════════
#> Preprocessor: Formula
#> Model: linear_reg()
#> 
#> ── Preprocessor ────────────────────────────────────────────────────────────────
#> survtime ~ x1
#> 
#> ── Model ───────────────────────────────────────────────────────────────────────
#> 
#> Call:
#> stats::lm(formula = ..y ~ ., data = data)
#> 
#> Coefficients:
#> (Intercept)       x1TRUE  
#>      121.24       -23.58  
#> 
#> 
#> [[2]]
#> ══ Workflow [trained] ══════════════════════════════════════════════════════════
#> Preprocessor: Formula
#> Model: linear_reg()
#> 
#> ── Preprocessor ────────────────────────────────────────────────────────────────
#> survtime ~ x2
#> 
#> ── Model ───────────────────────────────────────────────────────────────────────
#> 
#> Call:
#> stats::lm(formula = ..y ~ ., data = data)
#> 
#> Coefficients:
#> (Intercept)       x2TRUE  
#>     104.055       -7.706

#not working
models2 <- map(x,
               ~workflow() %>%
                 add_model(proportional_hazards()) %>%
                 add_formula(reformulate(.x, response = 'Surv(survtime, event=fail)')) %>%
                 fit(df)
)
#> Error in vctrs_vec_compat(.x, .purrr_user_env): object 'x' not found

#not working
models3 <- map(x,
               ~workflow() %>%
                 add_model(proportional_hazards()) %>%
                 add_formula(as.formula(paste0('Surv(survtime, event=fail) ~ ', .x))) %>%
                 fit(df)
)
#> Error in vctrs_vec_compat(.x, .purrr_user_env): object 'x' not found
models3
#> Error in eval(expr, envir, enclos): object 'models3' not found

创建于2023年1月6日,使用reprex v2.0.2

相关问题