Direct Standardization

To calculate confounder-adjusted survival curves fo each group of this dataset using Direct Standardization (also known as G-Computation or Corrected Group Prognosis method, see Makuch (1982) or Chang et al. (1982)), we first have to fit a model. Here we will use a coxph model, but we could also use pretty much any other kind of model.

# it is important to use X=TRUE in the coxph function call
outcome_model <- survival::coxph(Surv(time, event) ~ x1 + x2 + x3 + x4 + x5 + 
                                   x6 + group, data=data_1, x=TRUE)

This model can then be used in a call to the adjustedsurv() function, as shown below:

adjsurv <- adjustedsurv(data=data_1,
                        variable="group",
                        ev_time="time",
                        event="event",
                        method="direct",
                        outcome_model=outcome_model,
                        conf_int=TRUE)

The argument data simply refers to our data.frame, the argument variable specifies our grouping variable of interest and the ev_time and event variable specify the time-to-event variables in our data.frame. Setting the method to "direct" will result in G-Computation estimates, based on the previously fit cox-regression model supplied using the outcome_model argument.

Doing this returns a list with some needed output objects. Most important however is the adj data.frame in there, containing the adjusted survival curves and corresponding confidence intervals (because we used conf_int=TRUE in the original function call). We can take a look at this object using the following code:

head(adjsurv$adj)

##         time      surv group          se  ci_lower  ci_upper
## 1 0.00000000 1.0000000     0 0.000000000 1.0000000 1.0000000
## 2 0.01457659 0.9959356     0 0.004039858 0.9880177 1.0000000
## 3 0.04062032 0.9917639     0 0.005767859 0.9804591 1.0000000
## 4 0.04232678 0.9876128     0 0.007034862 0.9738247 1.0000000
## 5 0.06129069 0.9834495     0 0.008123567 0.9675276 0.9993714
## 6 0.06271999 0.9792709     0 0.009066121 0.9615017 0.9970402

More importantly however, we can plot the survival curves directly using the plot method:

plot(adjsurv)

This plot function comes with many options which are listed in the documentation. To plot the point-wise confidence intervals, we can set the argument conf_int to TRUE:

plot(adjsurv, conf_int=TRUE)

Many more custom settings, such as median survival time indicator lines, censoring indicator lines and risk tables are available. For more details and examples see ?plot.adjustedsurv or the associated vignette, which goes in even more detail on how the plots can be customized.

Inverse Probability of Treatment Weighting

The adjustedsurv() function essentially works the same with every available method. Since the methods are however vastly different in nature, some additional arguments have to be supplied by the user. For example, using the Inverse Probability of Treatment Weighting method (IPTW), we need to model the treatment assignment mechanism instead of the outcome mechanism (see Xie and Liu (2005)). This can be done in many ways. Below we use a logistic regression model for exemplary purposes:

treatment_model <- glm(group ~ x1 + x2 + x3 + x4 + x5 + x6,
                       data=data_1, family="binomial"(link="logit"))

adjsurv <- adjustedsurv(data=data_1,
                        variable="group",
                        ev_time="time",
                        event="event",
                        method="iptw_km",
                        treatment_model=treatment_model,
                        conf_int=TRUE)

The resulting curves can be plotted as before:

plot(adjsurv, conf_int=TRUE)

Since both methods were used correctly here, there are only slight differences in the results. Big differences between the two methods usually indicate that either the outcome model or the treatment-assignment model are incorrectly specified.

Augmented Inverse Probability of Treatment Weighting

Doubly-Robust methods can be helpful in such cases (Ozenne et al. 2020). The standard Augmented Inverse Probability of Treatment Weighting estimator utilized both kinds of models at the same time. If either of the models is correctly specified, the resulting estimates will be unbiased. It can be used in the same way as the other methods, only that this time both models have to be supplied. Using the same models as above we can use the following syntax:

adjsurv <- adjustedsurv(data=data_1,
                        variable="group",
                        ev_time="time",
                        event="event",
                        method="aiptw",
                        treatment_model=treatment_model,
                        outcome_model=outcome_model,
                        conf_int=TRUE)

plot(adjsurv, conf_int=TRUE)

Direct Standardization

Let’s again start with the Direct Standardization method. Instead of using a simple coxph method, we need to use a model for the time-to-event process which takes the multiple event-types into account. One such method is the Cause-Specific Cox-Regression model. A simple implementation of this model is contained in the riskRegression R-Package.

outcome_model <- riskRegression::CSC(Hist(time, event) ~ x1 + x2 + x3 + x4 + 
                                       x5 + x6 + group, data=data_2)

This model can then be used in a call to the adjustedcif() function, as shown below:

adjcif <- adjustedcif(data=data_2,
                      variable="group",
                      ev_time="time",
                      event="event",
                      method="direct",
                      outcome_model=outcome_model,
                      cause=1,
                      conf_int=TRUE)
plot(adjcif, conf_int=TRUE)

This shows the confounder-adjusted CIFs for cause = 1. By setting cause to 2 we get the confounder-adjusted CIFs for the other cause:

adjcif <- adjustedcif(data=data_2,
                      variable="group",
                      ev_time="time",
                      event="event",
                      method="direct",
                      outcome_model=outcome_model,
                      cause=2,
                      conf_int=TRUE)
plot(adjcif, conf_int=TRUE)

Inverse Probability of Treatment Weighting

The IPTW estimator can be used exactly the same way as we did with the adjustedsurv() function:

treatment_model <- glm(group ~ x1 + x2 + x3 + x4 + x5 + x6,
                       data=data_2, family="binomial"(link="logit"))

adjcif <- adjustedcif(data=data_2,
                      variable="group",
                      ev_time="time",
                      event="event",
                      method="iptw",
                      treatment_model=treatment_model,
                      cause=1,
                      conf_int=TRUE)
plot(adjcif, conf_int=TRUE)

Implemented Methods

For users who want to learn more about the implemented methods, we recommend consulting the documentation of the main functions first (?adjustedsurv or ?adjustedcif). This page will also guide the user to more detailed method-specific information. An overview of all implemented methods is given in the associated vignette (vignette(topic="method_overview", package="adjustedCurves")).

Plot Customization

The plot.adjustedsurv() function contains a lot of arguments to customize the adjusted survival curve plots, which may be daunting at first. The associated vignette (vignette(topic="plot_customization", package="adjustedCurves")) describes most of these options in great detail with multiple examples.

Comparing Groups

Users may also be interested in formally comparing the (treatment) groups after adjustment has been made. There are multiple ways to perform such comparisons using this package for both standard time-to-event data and data with competing events. A vignette explaining what these different options are and how to use them is also included in this package (vignette(topic="comparing_groups", package="adjustedCurves")). We recommend consulting this vignette if group comparisons are of interest.

FAQ

Some of the most frequently asked questions are answered in the FAQ vignette (vignette(topic="faq", package="adjustedCurves")).