Interventional prediction score for longitudinal treatment strategies

Estimates the predictive performance of predictions under longitudinal interventions, by forming a weighted pseudopopulation in which every subject was assigned the longitudinal treatment strategy of interest.

Usage

ip_score_long(
  probabilities,
  data_outcome,
  data_long,
  treatment_formula,
  treatment_of_interest,
  metrics = c("auc", "brier", "oeratio", "calplot"),
  visit_times,
  time_horizon,
  cens_model = "KM",
  cens_formula = ~1,
  null_model = TRUE,
  bootstrap = 0,
  bootstrap_progress = TRUE,
  iptw,
  ipcw,
  quiet = FALSE,
  strip_ipt_models = TRUE
)

Arguments

probabilities: A numeric vector corresponding to risk estimates under the treatment level of interest, or a (named) list of multiple numeric vectors for validating and comparing multiple risk estimates.
data_outcome: A dataframe containing the survival outcomes. It must consist of 1 row per subject, with columns 'id', 'time', and 'status'.
data_long: A dataframe in 'long' format, containing the time dependent treatment variable, the time dependent confounders and other time dependent covariates, possibly required for modeling the censoring mechanism. Baseline covariates can also be included, in which case they would be repeated for every visit. It is in 'long' format, i.e. each subject has one row for each of their visits. Must have the columns 'id' and 'visit_time'. There should not be any visits after a subject's survival time. Subject should not skip any visit.
treatment_formula: A formula which identifies the treatment (left hand side) and the confounders (right hand side) in data_long. E.g. A ~ L + A_lag_1. The treatment can be either binary (0/1) or a categorical factor. The confounders are used to estimate the inverse probability of treatment weights (IPTW) model. The IPTW can also be specified themselves using the iptw argument, in which case the right hand side of this formula is ignored (the left hand side must still identify the treatment, i.e. A ~ 1).
treatment_of_interest: A treatment strategy for which the interventional performance measures should be evaluated. Should be in the form of a vector, with one element per visit. See details.
metrics: A character vector specifying which performance metrics to be computed. Options are c("auc", "brier", "oeratio", "calplot").
visit_times: A numeric vector, indicating the times of the visits. The first visit must always be at t = 0.
time_horizon: the prediction horizon of interest.
cens_model: Model for estimating inverse probability of censored weights (IPCW). Methods currently implemented are Kaplan-Meier ("KM") or Cox ("cox"), both applied to the censored times. KM is only supported when the right hand side of cens_formula is 1.
cens_formula: Formula for which the r.h.s. determines the censoring probabilities. Could consist of only baseline variables but also of time-dependent variables. Either way, all variables specified must be present in data_long
null_model: If TRUE fit a risk prediction model which ignores the covariates and predicts the same value for all subjects. The model is fitted using the data in which all subjects are 'counterfactually' assigned the treatment of interest (using the IPTW, as estimated using the treatment_formula or as given by the iptw argument). For time-to-event outcomes, the subjects are also 'counterfactually' uncensored (using the IPCW, as estimated using the cens_formula, or as given by the ipcw argument).
bootstrap: If this is an integer greater than 0, this indicates the number of bootstrap iterations, to compute 95% confidence intervals around the performance metrics.
bootstrap_progress: if set to TRUE, print a progress bar indicating the progress of bootstrap procedure.
iptw: A numeric vector, containing the inverse probability of treatment weights. These are normally computed using the treatment_formula, but they can be specified directly via this argument. A function can also be specified, which takes as input arguments 'data_outcome' and 'data_long' and returns a numeric vector of the IPTW weight. See details.
ipcw: A numeric vector, containing the inverse probability of censoring weights. These are normally computed using the cens_formula, but they can be specified directly via this argument. A function can also be specified, which takes as input arguments 'data_outcome' and 'data_long' and returns a numeric vector of the IPCW weight. See details.
quiet: If set to TRUE, don't print assumptions.
strip_ipt_models: If set to TRUE (default), the models for the IPT and IPC-weights are stripped of unnecessary data. Set to FALSE if you plan to do extensive diagnostics on the fitted IPT/IPC models. The resulting `ip_score` object will use quite a lot more memory.

Value

An object of class `ip_score`, for which the `print()` and `plot()` methods are implemented. The object is a nested list containing:

`$score`, which contains the predictive performance in the pseudopopulation.
`$bootstrap`, if applicable, the 95% confidence intervals of the performance metrics, and the performance metrics for each individual bootstrap iteration.
`$outcome`, the observed outcome of the original dataset.
`$treatment`, the observed outcome of the original dataset.
`$predictions`, the predictions to be evaluated, i.e. the probability of event for each patient, had their treatment been set to treatment_of_interest.
`$ipt`, method, model and inverse probability of treatment weights (IPTW). These are NA for patients that are not in the pseudopopulation.
`$ipc`, method, model and inverse probability of censoring weights (IPCW). These are NA for patients that were censored.
`$pseudopop`, binary vector indicating which subjects of the original population were in the pseudopopulation, by following the treatment of interest and remaining uncensored.

The print method summarizes the results and if (quiet = FALSE), prints the assumptions required for valid inference.

Details

To form a pseudopopulation that represents a setting in which everybody was treated for all 5 visits, set `treatment_of_interest` to `c(1,1,1,1,1)`. Any pattern can be chosen, i.e. `c(1,0,1,0,1)` is also valid, as long as the number of visits is correct. Alternatively, one could also set this argument to "always" or "never" as a shortcut for `rep(1, n_visits)` or `rep(0, n_visits)`. Treatment strategies where only certain parts are set fixed are also possible via `NA`, i.e. use `c(1,1,NA, NA, NA)` when you want to form a pseudopopulation where everybody's first 2 treatment levels are set to 1, and their remaining 3 can be whatever they would have been. If treatment is categorical, this should be a character vector denoting the treatment levels of interest.

Bootstrap is not possible when manually specifiying the IPTW/IPCW as numeric vectors. If specifying a function that computes the ITPW/IPCW given data, it is possible. The given function will be called on each bootstrapped dataset to compute the 95 such as thresholding extreme IP weights, can be implemented this way. The censoring weight returned should be the 1 / probability of remaining uncensored at the time horizon, or at their event time, whichever happens first.

References

Keogh RH, Van Geloven N. Prediction Under Interventions: Evaluation of Counterfactual Performance Using Longitudinal Observational Data. Epidemiology. 2024;35(3):329-339.

Examples

set.seed(5)
n <- 1000
data <- data.frame(id = 1:n)

data <- within(data, {
  # 2 visits at t = 0 and t = 2. Time dependent confounding between A and L
  L0 <- rnorm(n)
  A0 <- rbinom(n, 1, plogis(0.7*L0))
  L1 <- rnorm(n, 0.8*L0 - 0.2*A0)
  A1 <- rbinom(n, 1, plogis(0.7*L1 + 0.6*A0))
  # Quickly generate random survival times as an example
  time <- runif(n, 0, 6)
  status <- rbinom(n, 1, 0.5)
})

# If you have data in this 'wide' form, then usually you can use the convenience
# functions to be able to run ip_score_long without doing a lot of extra
# processing
head(data)
#>   id status     time A1         L1 A0          L0
#> 1  1      0 4.094243  1  0.8174053  0 -0.84085548
#> 2  2      0 3.095679  0  1.0777666  1  1.38435934
#> 3  3      0 3.936240  0  1.2297761  1 -1.25549186
#> 4  4      1 5.390756  1  1.5875598  0  0.07014277
#> 5  5      0 1.245473  1  2.1049182  1  1.71144087
#> 6  6      0 4.491761  0 -1.0007181  1 -0.60290798
# note that due to our simulation there may be measurements after a subject was
# censored or had event. This will be removed later.

data_outcome <- data[, c("id", "time", "status")]
data_long <- wide_to_long(
  df = data,
  baseline_variables = "id",
  wide_variables = list("A" = c("A0", "A1"),
                        "L" = c("L0", "L1")),
  visit_times = c(0,2),
  outcome_times = data$time
)
data_long <- add_lag_terms(data_long, "A")

head(data_long)
#>   id visit_time A          L A_lag_1
#> 1  1          0 0 -0.8408555       0
#> 2  1          2 1  0.8174053       0
#> 3  2          0 1  1.3843593       0
#> 4  2          2 0  1.0777666       1
#> 5  3          0 1 -1.2554919       0
#> 6  3          2 0  1.2297761       1
# note that the measurements that were generated after survival time have been
# removed by wide_to_long().

# To validate, we need predictions. To get started, here is a model that
# randomly predicts a number between 0.2 and 0.8
random_model <- runif(n, min = 0.2, max = 0.8)

# Our treatment at visit i depends on confounder value L at same visit and,
# if visit = 2, on treatment level at previous visit. Thus, we will estimate our
# IPTW model with A ~ L + A_lag_1 in our long data.

# Estimate performance of our random predictions in a pseudopopulation in
# which everybody was assigned treatment strategy {0, 0}
ip_score_long(
  probabilities = random_model,
  data_outcome = data_outcome,
  data_long = data_long,
  treatment_formula = A ~ L + A_lag_1,
  treatment_of_interest = c(0,0),
  visit_times = c(0,2),
  time_horizon = 4
)
#> Estimation of the performance of the prediction model in a
#>  pseudopopulation where everyone's treatment A was set to {0, 0}.
#> The pseudopopulation ($pseudopop) is constructed from 217 (21.7%)
#>  subjects who indeed received treatment level {0, 0} and remained
#>  uncensored till time=4. These subjects are reweighted to represent the
#>  full target population under a hypothetical intervention in which
#>  everyone received this treatment level and remained uncensored till
#>  time=4.
#> The following assumptions must be satisfied for correct inference:
#> 
#> Causal assumptions:
#> 
#> - Conditional exchangeability: after adjustment for the covariates used
#>  to construct the inverse probability of treatment weights (IPTW), i.e.,
#>  {L, A_lag_1}, there are no unmeasured confounders of treatment
#>  assignment and outcome.
#> - Conditional positivity: the probability of receiving treatment level
#>  {0, 0} should be greated than zero for each combination of the
#>  variables {L, A_lag_1} that is observed in the full population. The
#>  distribution of IPT-weights can be assessed with
#>  $ipt$weights[$pseudopop$ids].
#> - Consistency: the observed outcome under the received treatment equals
#>  the potential outcome under that treatment. This includes the
#>  assumption of no interference between subjects.
#> - Independent censoring. The censoring mechanism is completely
#>  independent of the outcome process.
#> - Positivity for censoring: requires that the probability of remaining
#>  uncensored till time=4 is greater than zero. The distribution of
#>  IPC-weights can be assessed with $ipc$weights[$pseudopop$ids].
#> 
#> Modeling assumptions:
#> 
#> - Correctly specified propensity model. Estimated treatment model is
#> logit(A) = 0.05 + 0.72*L + 0.18*A_lag_1. See also $ipt$model.
#> - The censoring distribution was estimated nonparametrically using the
#>  Kaplan-Meier estimator. The probability of remaining uncensored is
#> P(C > 4) = 0.59. See also $ipc$model.
#> 
#> Performance estimates:
#> 
#>         model auc brier oeratio
#>    null model 0.5 0.246   1.000
#>  random_model 0.5 0.281   0.872


# Performance of our random predictions in a pseudopopulation in
# which everybody was assigned treatment level 1 at visit 1, and whatever
# they would normally have at visit 2.
ip_score_long(
  probabilities = random_model,
  data_outcome = data_outcome,
  data_long = data_long,
  treatment_formula = A ~ L + A_lag_1,
  treatment_of_interest = c(1,NA),
  visit_times = c(0,2),
  time_horizon = 4
)
#> Estimation of the performance of the prediction model in a
#>  pseudopopulation where everyone's treatment A was set to {1, *}, where
#>  * can be any value as would normally be observed.
#> The pseudopopulation ($pseudopop) is constructed from 330 (33%)
#>  subjects who indeed received treatment level {1, *} and remained
#>  uncensored till time=4. These subjects are reweighted to represent the
#>  full target population under a hypothetical intervention in which
#>  everyone received this treatment level and remained uncensored till
#>  time=4.
#> The following assumptions must be satisfied for correct inference:
#> 
#> Causal assumptions:
#> 
#> - Conditional exchangeability: after adjustment for the covariates used
#>  to construct the inverse probability of treatment weights (IPTW), i.e.,
#>  {L, A_lag_1}, there are no unmeasured confounders of treatment
#>  assignment and outcome.
#> - Conditional positivity: the probability of receiving treatment level
#>  {1, *} should be greated than zero for each combination of the
#>  variables {L, A_lag_1} that is observed in the full population. The
#>  distribution of IPT-weights can be assessed with
#>  $ipt$weights[$pseudopop$ids].
#> - Consistency: the observed outcome under the received treatment equals
#>  the potential outcome under that treatment. This includes the
#>  assumption of no interference between subjects.
#> - Independent censoring. The censoring mechanism is completely
#>  independent of the outcome process.
#> - Positivity for censoring: requires that the probability of remaining
#>  uncensored till time=4 is greater than zero. The distribution of
#>  IPC-weights can be assessed with $ipc$weights[$pseudopop$ids].
#> 
#> Modeling assumptions:
#> 
#> - Correctly specified propensity model. Estimated treatment model is
#> logit(A) = 0.05 + 0.72*L + 0.18*A_lag_1. See also $ipt$model.
#> - The censoring distribution was estimated nonparametrically using the
#>  Kaplan-Meier estimator. The probability of remaining uncensored is
#> P(C > 4) = 0.59. See also $ipc$model.
#> 
#> Performance estimates:
#> 
#>         model   auc brier oeratio
#>    null model 0.500 0.249   1.000
#>  random_model 0.512 0.276   0.933