-
Notifications
You must be signed in to change notification settings - Fork 13
Open
Description
library('tidyverse')
library('survival')
# example from https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf
dim(cgd0)
#> [1] 128 20
newcgd <- tmerge(data1=cgd0[, 1:13], data2=cgd0, id=id, tstop=futime)
newcgd <- tmerge(newcgd, cgd0, id=id, infect = event(etime1))
newcgd <- tmerge(newcgd, cgd0, id=id, infect = event(etime2))
newcgd <- tmerge(newcgd, cgd0, id=id, infect = event(etime3))
newcgd <- tmerge(newcgd, cgd0, id=id, infect = event(etime4))
newcgd <- tmerge(newcgd, cgd0, id=id, infect = event(etime5))
newcgd <- tmerge(newcgd, cgd0, id=id, infect = event(etime6))
newcgd <- tmerge(newcgd, cgd0, id=id, infect = event(etime7))
newcgd <- tmerge(newcgd, newcgd, id, enum=cumtdc(tstart))
dim(newcgd)
#> [1] 203 17
newcgd[1:5,c(1, 4:6, 13:17)]
#> id treat sex age futime tstart tstop infect enum
#> 1 1 1 2 12 414 0 219 1 1
#> 2 1 1 2 12 414 219 373 1 2
#> 3 1 1 2 12 414 373 414 0 3
#> 4 2 0 1 15 439 0 8 1 1
#> 5 2 0 1 15 439 8 26 1 2
summary(newcgd)
#> Call:
#> tmerge(data1 = newcgd, data2 = newcgd, id = id, enum = cumtdc(tstart))
#>
#> early late gap within boundary leading trailing tied missid
#> infect 0 0 0 44 0 0 0 0 0
#> infect 0 0 0 16 0 0 1 0 0
#> infect 0 0 0 8 0 0 0 0 0
#> infect 0 0 0 3 0 0 0 0 0
#> infect 0 0 0 2 0 0 0 0 0
#> infect 0 0 0 1 0 0 0 0 0
#> infect 0 0 0 1 0 0 0 0 0
#> enum 0 0 0 0 75 128 0 0 0
newcgd <- newcgd %>%
mutate(outcome = Surv(tstart, tstop, infect))
set.seed(1234)
library('tidymodels')
tidymodels_prefer()
library('censored')
# split the data
ini_split <- group_initial_split(newcgd, group = id)
train <- training(ini_split)
test <- testing(ini_split)
# we can train cox model
coxph_model <- coxph(Surv(tstart, tstop, infect) ~ treat + inherit + steroids,
data =train,
cluster = id)
coxph_model %>%
broom::tidy()
#> # A tibble: 3 × 6
#> term estimate std.error robust.se statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 treat -1.17 0.309 0.361 -3.23 0.00123
#> 2 inherit -0.0762 0.265 0.355 -0.215 0.830
#> 3 steroids -0.902 0.528 0.386 -2.34 0.0195
# we can use glmnet
library('glmnet')
#> Loading required package: Matrix
#>
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#>
#> expand, pack, unpack
#> Loaded glmnet 4.1-9
X <- train %>%
select(treat, inherit, steroids) %>%
as.matrix()
Y_M <- train %>%
select(outcome) %>%
as.matrix()
Y <- Surv(Y_M[,'outcome.start'], Y_M[,'outcome.stop'], Y_M[,'outcome.status'])
glmnet_model <- cv.glmnet(x=X, y=Y, family = 'cox', nfolds = 5)
#> Warning: cox.fit: algorithm did not converge
plot(glmnet_model)
glmnet_model <- glmnet(x=X, y=Y, family = 'cox', lambda = 10^(-3))
#> Warning: cox.fit: algorithm did not converge
glmnet_model %>%
broom::tidy()
#> # A tibble: 3 × 5
#> term step estimate lambda dev.ratio
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 treat 1 -1.18 0.001 0.0377
#> 2 inherit 1 0.0395 0.001 0.0377
#> 3 steroids 1 -0.00850 0.001 0.0377
# now let's try tidymodel's approach:
library('workflowsets')
folds <- group_vfold_cv(train, v = 5, group = id)
rec <- recipe(outcome ~ treat + inherit + steroids, data = train)
cox_spec <- proportional_hazards(
mode = "censored regression",
engine = "glmnet",
penalty = tune(),
mixture = tune()
)
wflowset <- workflow_set(
preproc = list(recipe = rec),
models = list(cox = cox_spec)
)
grid_ctrl <-
control_grid(
save_pred = TRUE,
parallel_over = "everything",
save_workflow = TRUE
)
# does not work, expects Surv(time, status) type of outcome:
grid_results <-
wflowset %>%
workflow_map(
seed = 3110,
resamples = folds,
grid = 35,
control = grid_ctrl,
eval_time = c(15,30,60,90,120),
metrics = metric_set(brier_survival, roc_auc_survival, brier_survival_integrated)
)
#> → A | warning: cox.fit: algorithm did not converge
#> There were issues with some computations A: x1 → B | error: For this usage, the allowed censoring type is "right".
#> Warning: All models failed. Run `show_notes(.Last.tune.result)` for more
#> information.
Created on 2025-06-05 with reprex v2.1.1
Metadata
Metadata
Assignees
Labels
No labels