robindenz1 / adjustedcurves Goto Github PK

Is it possible to create the complement of the survival probabilities i.e. to plot adjusted failure curves instead of the Kaplan Meier survival plot?

Thank you for your terrific work with this package; it is incredibly helpful and thorough! Was delighted when I found this package.

My question pertains to the adjustedcif function when one is using method=iptw. Is there any way to pass a vector of weights instead of a model to the treatment_model argument? This would be similar to the functionality offered by method=iptw for the adjustedsurv function. Asking because it would be great to estimate weights outside of the package and then just supply the user-created weights to the treatment_model argument of adjustedcif.

Reasons for wanting to externally generate weights: ability to try different packages for weight creation (WeightIt versus PSWeight), trim weights, assess covariate balance, etc.

Thanks for this contribution. Under my brief testing I needed to start with a clean R install with the latest R (4.1.3 Windows). It seemed I needed to rename time and event columns in my dataset to be "time" and "event". Is there a best way to compute the overall adjusted curve (collapsed across all variables)? Would the mean of all variables/categories be valid?

Hi,
many thanks for this package and for your work!

I'm using adjustedcif for plotting curves for a competing risk outcome with "iptw_pseudo" and providing custom weights as explained in cif_iptw_pseudo examples.

I have few (probably silly) questions for result visualization:

• could it be possibile to plot not the cif curve but 1-cif (kind of survival probability)? I do not see any fun argument in plot.adjustedcif that could do that
• i'd like to add risk table below the plot, but i do not think the function can address this issue. Could it be possibile to extract the code under plot.adjustedcif so that i can customize it accordingly (e.g. plotting curves with ggsurvfit package, that would easily allow to add risk tables, etc etc).
  Thanks for your support.

Hi all -- I am getting an error when trying to use a survreg object which contains an accelerated failure time (AFT) survival model. Any idea why this error may be thrown? Is adjustedCurves::adjustedsurv unable to handle survreg objects?

test_grouped_alt <- survreg(Surv(time_overall, sx) ~  age_firstvisit_c + 
                             gender + 
                             race_char +
                             group +
                             payer,
                           cluster=studyid,
                           data= as.data.frame(surv_df),
                           dist = 'lognormal')

survtest <- adjustedCurves::adjustedsurv(data = as.data.frame(surv_df), 
                            variable = "group",
                            ev_time = "time_overall",
                            event = "sx",
                            method = "direct",
                            outcome_model = test_grouped_alt,
                            cause = 1,
                            conf.int = FALSE)

Error noted:
Error in value[[3L]](cond) : The following error occured using the default S3 predict method: 'Error in eval(predvars, data, env): object 'age_firstvisit_c' not found ' Specify a valid 'predict_fun' or use a different model. See details.

Hi,

Is it possible to change the interval of x-axis?

Thank you.

Hi, is it possible to make this work with other ggsurvplot functionality such as theme and facet_grid(variable)?
Thanks for the confidence intervals (super handy!)

Hi, I am so happy to meet your good package.
But my data has weights in TWANG package like survey weights
Is this also applicable to adjustedsurv functions?

I want to use this weight directly to cumulative incidence plot.
Unfortunately, my data is too heavy.....

hope you have a nice day.
Sincerely,
Coconning Kang

I installed from github and am getting an error with the direct method if I create the Surv object beforehand:
Error: node stack overflow
Error during wrapup: node stack overflow
Error: no more error handlers available (recursive errors?); invoking 'abort' restart

to reproduce:

set.seed(31)
sim_dat <- adjustedCurves::sim_confounded_surv(n=2500, max_t=1.2, group_beta=0)
sim_dat$group <- as.factor(sim_dat$group)

#pre-creating the Surv object
sim_dat$survTbl <- survival::Surv(sim_dat$time, sim_dat$event)

# estimate a cox-regression for the outcome
cox_mod <- survival::coxph(survTbl ~ x1 + x2 + x4 + x5 + group,
                 data=sim_dat, x=TRUE)

# use it to calculate adjusted survival curves
adjsurv <- adjustedCurves::adjustedsurv(data=sim_dat,
                        variable="group",
                        ev_time="time",
                        event="event",
                        method="direct",
                        outcome_model=cox_mod,
                        conf_int=TRUE)

# plot with confidence intervals
plot(adjsurv, conf_int=TRUE)

SessionInfo:

R version 4.1.3 (2022-03-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=Estonian_Estonia.1257  LC_CTYPE=Estonian_Estonia.1257    LC_MONETARY=Estonian_Estonia.1257
[4] LC_NUMERIC=C                      LC_TIME=Estonian_Estonia.1257    
system code page: 1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
  [1] nlme_3.1-160               cmprsk_2.2-11              rms_6.3-0                  RColorBrewer_1.1-3        
  [5] numDeriv_2016.8-1.1        tools_4.1.3                backports_1.4.1            doRNG_1.8.2               
  [9] utf8_1.2.2                 R6_2.5.1                   rpart_4.1.19               Hmisc_4.7-0               
 [13] DBI_1.1.3                  colorspace_2.0-3           nnet_7.3-18                tidyselect_1.1.2          
 [17] gridExtra_2.3              compiler_4.1.3             cli_3.3.0                  quantreg_5.94             
 [21] htmlTable_2.4.1            SparseM_1.81               sandwich_3.0-2             scales_1.2.0              
 [25] checkmate_2.1.0            polspline_1.1.20           survMisc_0.5.6             mvtnorm_1.1-3             
 [29] stringr_1.4.0              digest_0.6.29              foreign_0.8-83             R.utils_2.12.0            
 [33] base64enc_0.1-3            jpeg_0.1-9                 pkgconfig_2.0.3            htmltools_0.5.3           
 [37] parallelly_1.32.1          adjustedCurves_0.10.1.9000 fastmap_1.1.0              htmlwidgets_1.5.4         
 [41] rlang_1.0.4                rstudioapi_0.13            generics_0.1.3             zoo_1.8-10                
 [45] R.oo_1.25.0                dplyr_1.0.9                car_3.1-0                  magrittr_2.0.3            
 [49] Formula_1.2-4              interp_1.1-3               Matrix_1.5-1               Rcpp_1.0.9                
 [53] munsell_0.5.0              fansi_1.0.3                abind_1.4-5                R.methodsS3_1.8.2         
 [57] lifecycle_1.0.1            multcomp_1.4-20            stringi_1.7.6              carData_3.0-5             
 [61] MASS_7.3-58.1              grid_4.1.3                 parallel_4.1.3             listenv_0.8.0             
 [65] crayon_1.5.1               mets_1.2.9                 survminer_0.4.9            deldir_1.0-6              
 [69] lattice_0.20-45            splines_4.1.3              knitr_1.39                 pillar_1.8.0              
 [73] timereg_2.0.2              ggpubr_0.4.0               ggsignif_0.6.3             rngtools_1.5.2            
 [77] future.apply_1.9.0         codetools_0.2-18           glue_1.6.2                 latticeExtra_0.6-30       
 [81] data.table_1.14.2          png_0.1-7                  vctrs_0.4.1                foreach_1.5.2             
 [85] MatrixModels_0.5-0         gtable_0.3.0               purrr_0.3.4                tidyr_1.2.0               
 [89] future_1.27.0              km.ci_0.5-6                assertthat_0.2.1           riskRegression_2022.03.22 
 [93] ggplot2_3.3.6              xfun_0.32                  prodlim_2019.11.13         xtable_1.8-4              
 [97] broom_1.0.0                rstatix_0.7.0              survival_3.4-0             tibble_3.1.8              
[101] iterators_1.0.14           KMsurv_0.1-5               writexl_1.4.2              cluster_2.1.4             
[105] lava_1.6.10                globals_0.16.0             ellipsis_0.3.2             TH.data_1.1-1  

Hi,

Are the estimands ATE or ATT?
Thank you for the great package!

How do I get the p-value for the comparison between groups after using adjustedsurv?

Hi,

First, many thanks for writing this package - it is of great value to me.

If I understand correctly the package can only handle survival data of the form Surv(time, status).

Would it be a lot of work to generalize that for Surv(start,stop,status) data?

Hi,
Is there a way to put the numbers at risk table below the adjusted curves and to check for the SMD differences for the "ps" based analysis?

I have a coxph model e.g. :

fit <- coxph(Surv(survival_months, died) ~ prot_expr + gender + disease_onset + age_at_onset, data = df)

Generate a grouping variable based on the protein expression (prot_expr) - high vs low expression

min_exp <- min(df$prot_expr)
max_exp <- max(df$prot_expr)
df$expr_lvl <- cut(df$prot_expr, c(min_exp,prot_cutpoint,max_exp))

I want the variable argument to be the expr_lvl :

survival_est <- adjustedsurv(data=df ,outcome_model= fit,ev_time="survival_months", event="died",variable ="expr_lvl",method="direct",conf_int=TRUE)

Unfortunatelly when I run the plot command I get this error:

Error in check_inputs_adjustedsurv(data = data, variable = variable, ev_time = ev_time,  : 
  'variable' has to be included in the cox-regression model.

Is there a way to bypass this issue ?

The plotting sometimes generates a warning that comes likely from the fact that an underling ggplot2 function has removed the linewith parameter in the latest version.

Thank you for creating a great package!

I wanted to ask if it is possible to add p-value in the plot?

I noticed that neither ggadjustedcurves of survminer nor adjustedsurv function has a parameter for adding p-value.

Is it feasible to extract the p-value of the Cox model and add it to the plot object and then render it?

Thank you for your help!

Hello, thank you very much for your preliminary work. We have a problem when using this package, it seems that this package cannot add risk table, can we request to add this function?

Hi,
Thanks for sharing this wonderful package.
I wonder if it is possible to get an adjusted survival curve with a number at risk table?
Or are there other ways to get an adjusted survival curve with a number at risk table?
Thank you.

Hi,
when using method = direct with a cox model obtained using cozph and strata() on a variable with non-proportional hazard, I obtain an outcome with NA value in the adjusted survival. Of interest, these NA values are reported at the end of the curves for each group (3 NA values, one for each group, last value of the curve).
Thanks,
Pier Paolo

Hello,
Many thanks for your contribution in developing such a great package!

I am facing an error when using the adjustedsurv() function, which says
"Error in app$vspace(new_style$margin-top %||% 0) :
attempt to apply non-function".

I am wondering what could be the reason for this error, and what can be done to fix it? The detailed code and outputs are listed below.

Thank you for your generous help!

Code:
mod4_Cox1 <- coxph(Surv(D28_OUTCOME_DAYS,D28_DEATH)~STATUS2+AGE+SEX+SMOKE2+CHD+HT+DM+ND+LD+CPD+SCORE+BARICITINIB+WH_AV,data = temp2,x=T)
summary(mod4_Cox1)

adj4 <- adjustedCurves::adjustedsurv(data=temp2, variable="STATUS2", ev_time="D28_OUTCOME_DAYS", event="D28_DEATH", method="direct", outcome_model=mod4_Cox1, conf_int = T)

SessionInfo:

Call:
coxph(formula = Surv(D28_OUTCOME_DAYS, D28_DEATH) ~ STATUS2 +
AGE + SEX + SMOKE2 + CHD + HT + DM + ND + LD + CPD + SCORE +
BARICITINIB + WH_AV, data = temp2, x = T)

adj4 <- adjustedCurves::adjustedsurv(data=temp2,
variable="STATUS2",
ev_time="D28_OUTCOME_DAYS",
event="D28_DEATH",
method="direct",
outcome_model=mod4_Cox1,
conf_int = T)

Loading required namespace: riskRegression
Error in app$vspace(new_style$margin-top %||% 0) :
attempt to apply non-function

Hello,

Thank you for creating this package, it has been immensely helpful!

I have been running the adjustedcurves package for a specific dataset. When running the outcome model and the treatment model separately, there does not seem to be any issues. However, when I try to run the AIPTW model combining both the outcome and treatment model, I encounter the following error:

Error in table(object$data[[all.vars(formula(object))[1]]], object$y) :
all arguments must have the same length

I was not able to find the specific source of the error and I was unsure how to go about debugging this. Any help would be appreciated.

Additionally, I wanted to know if it would be possible to extract the weighted dataset when applying the treatment model to adjust the survival curve.

Thank you

Hello (again),

Really enjoy using this package and thankful for your efforts with adjustedCurves! I've noticed that adjusted_curve_diff and adjusted_curve_ratio will give NA values for the quantities they estimate if the times argument has a value is beyond when the last event occurred. As long as at least one individual is still under follow-up, one should (ideally) be able to get estimates for these values. For instance running:

adjusted_curve_ratio(adjsurv, times=1825, group_1="1", group_2="0", conf_int=TRUE)

will produce NAs for the ratio because the last event was at day 1450. However, the longest time to censoring is 2000, so the risk ratio should be constant between 1450 and 2000 (but still calculatable at day 2000). This will happen anytime the max censoring time is past the max event time. The current workaround is to: 1) Find the last event time 2) calculate the RR at that event time, 3) find the last follow-up date, 4) carry the RR forward to the last follow-up date.

Not sure if this an issue or enhancement, but would be awesome when possible to look into or add this feature. :)

Thanks again!

Charley

Thank you for creating a great package!

I wanted to ask if it is possible to plot hazards instead of survival curves?

Thank you for your help!

This is a very exellent package. Could this package achieve covariate adjusted KM survival curves for competing risk events analysis?

Hello!

Thank you for the great package!
Unfortunately, I unable to plot confidence intervals for my flexsurv model despite specifying bootsrap=TRUE. My flexsurv model and script for adjustedsurv are as following:

m1=flexsurvspline(Surv(time.years,
ACS.cases) ~ SSc,
data = df2,
k = 6, scale="hazard")
p1=adjustedsurv(data=df2,variable="ssc",ev_time="time.years",
event="ACS.cases",method="direct",outcome_model=m1, conf_int=TRUE, bootstrap=TRUE,n_boot=1000, n_cores=18)
plot(p1, conf_int=TRUE, steps=FALSE)
Warning message:
Cannot draw confidence intervals. Either set 'conf_int=TRUE' in adjustedsurv() call or use bootstrap estimates.
Could you help me with this?

Thank you!

Hi,
Would it be possible to output a table of at risk / events below the adjusted cumulative incidence / survival plot. Similar to those that can be produced in tidycmprsk:
https://mskcc-epi-bio.github.io/tidycmprsk/
I guess these at risk / events numbers will be predicted values as opposed to the raw data if they were to "align" with the adjusted curve.
Thanks.

Rstudio is crashing alternatively got a warning on memory. its a dataset with approximately 250 000 subjects. The variable Exercise is 5 category factor variable. Total amount of events i about 10 000 and median followup time is 37 yrs.

It works when using a sliced dataset of 20 000 instead of 250 000.

Also tried setting R_MAX_VSIZE=300Gb

Code used

cox_model <- coxph(Surv(risk_time_20, ICD_I_incidenceORdeath) ~ Exercise + Sex, data = data, x=TRUE)

adjustedsurv(data=data,
                        variable="Exercise",
                        ev_time="risk_time_20",
                        event="ICD_I_incidenceORdeath",
                        method="direct",
                        outcome_model=cox_model,
                        conf_int=T,
                        na.action = "na.omit",
                        n_cores = 2)

I have a new macbook pro with 36gb ram. Is the only solution to find a more powerful computer?

Sorry for no reproducible example.

This is a plot with 20 000 individuals from the original 250 000 individuals: