Movatterモバイル変換

Causal Inference with Continuous (Multiple Time Point) Interventions

Description

This package facilitates the estimation of counterfactual outcomes for multiple values of continuous interventions at different time points, and allows plotting of causal dose-response curves.

It implements the standard g-methods approach using the (semi-)parametricg-formula, as described in the Schomaker et al. (2024) reference listed below. Weighted dose-response curves that address positivity violations, and are fitted via sequentialg-computation, are currently only available on GitHub and are not (yet) integrated in this package.

The main function of the package is currentlygformula.

Details

Author(s)

Michael Schomaker

Maintainer: Michael Schomaker <michael.schomaker@stat.uni-muenchen.de>

References

Schomaker M, McIlleron H, Denti P, Diaz I. (2024)Causal Inference for Continuous Multiple Time Point Interventions, Statistics in Medicine, 43:5380-5400, see alsohttps://arxiv.org/abs/2305.06645.

Pharmacoepidemiological HIV treatment data

Description

A hypothetical, simulated dataset which is line with the data-generating process of Schomaker et al. (2024) and inspired by the data of Bienczak et al. (2017); see references below.

Usage

data(EFV)

Format

A data frame with 5000 observations on the following variables:

sex

The patient's sex

metabolic

Metabolism status (slow, intermediate, extensive) related to the single nucleotide polymorphisms in the CYP2B6 gene, which is relevantfor metabolizing evafirenz and directly affects its concentration in the body.

log_age

log(age) at baseline

NRTI

Nucleoside reverse transcriptase inhibitor (NRTI) component of HIV treatment, i.e.abacavir, stavudine or zidovudine.

weight.0

log(weight) at time 0 (baseline)

efv.0

Efavirenz concentration at time 0 (baseline)

VL.0

Elevated viral load (viral failure) at time 0 (baseline)

adherence.1

Adherence at time 1 (if 0, then signs of non-adherence)

weight.1

log(weight) at time 1

efv.1

Efavirenz concentration at time 1

VL.1

Elevated viral load (viral failure) at time 1

adherence.2

Adherence at time 2 (if 0, then signs of non-adherence)

weight.2

log(weight) at time 2

efv.2

Efavirenz concentration at time 2

VL.2

Elevated viral load (viral failure) at time 2

adherence.3

Adherence at time 3 (if 0, then signs of non-adherence)

weight.3

log(weight) at time 3

efv.3

Efavirenz concentration at time 3

VL.3

Elevated viral load (viral failure) at time 3

adherence.4

Adherence at time 4 (if 0, then signs of non-adherence)

weight.4

log(weight) at time 4

efv.4

Efavirenz concentration at time 4

VL.4

Elevated viral load (viral failure) at time 4

References

Schomaker M, McIlleron H, Denti P, Diaz I. (2024)Causal Inference for Continuous Multiple Time Point Interventions, Statistics in Medicine, 43:5380-5400, see alsohttps://arxiv.org/abs/2305.06645.

Bienczak et al. (2017)Determinants of virological outcome and adverse events in African children treated with paediatric nevirapine fixed-dose-combination tablets, AIDS,31:905-915

Examples

data(EFV)str(EFV)

Pharmacoepidemiological HIV treatment data

Description

A hypothetical, simulated dataset which is line with the data-generating process of Schomaker et al. (2024) and inspired by the data of Bienczak et al. (2017); see references below. Compared to the datasetEFV, it contains all variables of the DAG in Figure 3 of Schomaker et al. (2023), also those which are not needed for identification of the counterfactual quantity of interest; that is, the expected viral suppression (VL) under a specific intervention on efavirenz concentrations (efv.0, efv.1, ...).

Usage

data(EFVfull)

Format

A data frame with 5000 observations on the following variables:

sex

The patient's sex

metabolic

Metabolism status (slow, intermediate, extensive) related to the single nucleotide polymorphisms in the CYP2B6 gene, which is relevantfor metabolizing evafirenz and directly affects its concentration in the body.

log_age

log(age) at baseline

NRTI

Nucleoside reverse transcriptase inhibitor (NRTI) component of HIV treatment, i.e.abacavir, stavudine or zidovudine.

weight.0

log(weight) at time 0 (baseline)

comorbidity.0

Presence of co-morbidities at time 0 (baseline)

dose.0

Dose of efavirenz administered at time 0 (basline)

efv.0

Efavirenz concentration at time 0 (baseline)

VL.0

Elevated viral load (viral failure) at time 0 (baseline)

adherence.1

Adherence at time 1 (if 0, then signs of non-adherence)

weight.1

log(weight) at time 1

comorbidity.1

Presence of co-morbidities at time 1

dose.1

Dose of efavirenz administered at time 1

efv.1

Efavirenz concentration at time 1

VL.1

Elevated viral load (viral failure) at time 1

adherence.2

Adherence at time 2 (if 0, then signs of non-adherence)

weight.2

log(weight) at time 2

comorbidity.2

Presence of co-morbidities at time 2

dose.2

Dose of efavirenz administered at time 2

efv.2

Efavirenz concentration at time 2

VL.2

Elevated viral load (viral failure) at time 2

adherence.3

Adherence at time 3 (if 0, then signs of non-adherence)

weight.3

log(weight) at time 3

comorbidity.3

Presence of co-morbidities at time 3

dose.3

Dose of efavirenz administered at time 3

efv.3

Efavirenz concentration at time 3

VL.3

Elevated viral load (viral failure) at time 3

adherence.4

Adherence at time 4 (if 0, then signs of non-adherence)

weight.4

log(weight) at time 4

comorbidity.4

Presence of co-morbidities at time 4

dose.4

Dose of efavirenz administered at time 4

efv.4

Efavirenz concentration at time 4

VL.4

Elevated viral load (viral failure) at time

References

Schomaker M, McIlleron H, Denti P, Diaz I. (2024)Causal Inference for Continuous Multiple Time Point Interventions, Statistics in Medicine, 43:5380-5400, see alsohttps://arxiv.org/abs/2305.06645.

Bienczak et al. (2017)Determinants of virological outcome and adverse events in African children treated with paediatric nevirapine fixed-dose-combination tablets, AIDS,31:905-915

Examples

data(EFVfull)str(EFVfull)

Counterfactual contrast from the parametric or sequential g-formula for continuous multiple time point interventions

Description

Estimation of a contrast between counterfactual outcomes under different values of (continuous) interventions, or across different time points, using the parametric or sequential g-formula.

Usage

contrast(X, abar, nodes, contrastType = "difference", measure = mean,         cond = NULL, cilevel = 0.95, ...)

Arguments

Details

Causal effects are defined as contrasts between the distributions of counterfactual variables under different interventions, across different time points or across different covariate strata. The counterfactual distributions to compare must be uniquely determined, by either specifying two values ofabar at a singlenodes or twonodes at a single intervention levelabar or the natural course scenario withabar = 'natural' or two covariate stratacond. If the natural course scenario is selected and twonodes are specified, the natural intervention is compared across the two nodes. If onenodes is specified, the natural and observed scenarios are compared at a single node.

By default, the difference between the expectations of the two counterfactual outcome distributions is calculated. The difference can be exchanged for a ratio, odds ratio or custom contrast in thecontrastType argument, and expectations can be exchanged for custom measures in themeasure argument. Conditional measures can be specified through thecond argument. Custom contrasts, including those comparing more than two counterfactuals, can be defined by passing a function tocontrastType.

Confidence intervals are based on the nonparametric bootstrap withB samples.

Value

Returns a list of classcontrastResult:

See Also

gformula andsgf for estimating expected counterfactual outcomes under multiple intervention values andcustom.measure for measures other than expectations.

Examples

## Not run: data(EFV)gf1 <- gformula(  X = EFV, Anodes = c("efv.0", "efv.1", "efv.2", "efv.3", "efv.4"),  Ynodes = c("VL.0", "VL.1", "VL.2", "VL.3", "VL.4"),  Lnodes = c("adherence.1", "weight.1", "adherence.2", "weight.2",             "adherence.3", "weight.3", "adherence.4", "weight.4"),  abar = seq(1, 5), B = 10, ret = TRUE)# compare outcomes at last time point under (1,...,1) and (5,...,5) contrast(gf1, abar = c(1, 5), nodes = "VL.4")# compare outcomes at different time points, for same intervention (2,...)contrast(gf1, abar = 2, nodes = c("VL.3", "VL.2"))# compare own measure (rel. risk reduction) instead of mean# ... and conditional on subsetrelativeRiskReduction <- function(k, l){(k - l) / k}contrast(  gf1, abar = c(1, 2), nodes = "VL.4",  contrastType = relativeRiskReduction,  cond = quote(sex == 1))# Instead of the mean, any other measure can be taken,# and - of course - applied also to counterfactual Lnodescontrast(  gf1, abar = 2, nodes = c("weight.3", "weight.2"),  measure = median)## End(Not run)

Custom estimands after applyinggformula

Description

The default estimate returned bygformula is theexpected outcome under the respective intervention strategiesabar.custom.measure takes an object of classgformula and enables estimation of other estimands based on the counterfactual datasets produced bygformula (if the optionret=TRUE had been chosen), for example estimands conditional on baseline variables, quantiles instead of expectations, and others.

Usage

custom.measure(X, fun = NULL, cond = NULL, verbose = TRUE, with.se = FALSE, ...)

Arguments

Details

In settings with censoring, it will often be needed to pass on the optionna.rm=T, e.g. for the mean, median, quantiles, and others.

Calculation of the bootstrap standard error (i.e.,with.se=T) is typically not needed; but, for example, necessary for the calculations after multiple imputation and hence used bymi.boot.

Value

An object of classgformula. Seegformula for details.

See Also

see alsogformula

Examples

data(EFV)est <- gformula(X=EFV,                Lnodes  = c("adherence.1","weight.1",                            "adherence.2","weight.2",                            "adherence.3","weight.3",                            "adherence.4","weight.4"                ),                Ynodes  = c("VL.0","VL.1","VL.2","VL.3","VL.4"),                Anodes  = c("efv.0","efv.1","efv.2","efv.3","efv.4"),                abar=seq(0,2,1), ret=TRUE)estcustom.measure(est, fun=prop,categ=1) # identicalcustom.measure(est, fun=prop,categ=0)custom.measure(est, fun=prop, categ=0, cond="sex==1")# note: metabolic has been recoded internally (see output above)custom.measure(est, fun=prop, categ=0, cond="metabolic==0") # does not make sense here, just for illustration (useful for metric outcomes)custom.measure(est, fun=quantile, probs=0.1)

Fit models after screening

Description

Fits the models that have been generated with screening usingmodel.formulas.update.

Usage

fit.updated.formulas(formulas, X)

Arguments

Details

Fits generalized (additive) linear models based on the screened model formula list generated bymodel.formulas.update.

Value

Returns a list of length 2:

See Also

model.formulas.update

Examples

data(EFV)# first: generate generic model formulasm <- make.model.formulas(X=EFV,                         Lnodes  = c("adherence.1","weight.1",                                     "adherence.2","weight.2",                                     "adherence.3","weight.3",                                     "adherence.4","weight.4"                                    ),                         Ynodes  = c("VL.0","VL.1","VL.2","VL.3","VL.4"),                         Anodes  = c("efv.0","efv.1","efv.2","efv.3","efv.4"),                         evaluate=FALSE)                          # second: update these model formulas based on variable screening with LASSOglmnet.formulas <-  model.formulas.update(m$model.names, EFV)glmnet.formulas # then: fit and inspect the updated modelsfitted.models <- fit.updated.formulas(glmnet.formulas, EFV)fitted.models$all.summariesfitted.models$all.summaries$Ynames[1] # first outcome model

Parametric g-formula for continuous multiple time point interventions

Description

Estimation of counterfactual outcomes for multiple values of continuous interventions at different time points using the g-formula.

Usage

gformula(X, Anodes, Ynodes, Lnodes = NULL, Cnodes = NULL,         abar = NULL, cbar = "uncensored",         survivalY = FALSE,          Yform = "GLM", Lform = "GLM", Aform = "GLM", Cform = "GLM",         calc.support = FALSE, B = 0, ret = FALSE, ncores = 1,          verbose = TRUE, seed = NULL, prog = NULL, ...)

Arguments

Details

By default, expected counterfactual outcomes (specified underYnodes) under the interventionabar are calculated. Other estimands can be specified viacustom.measure.

Ifabar is a vector, then each vector component is used as the intervention value at each time point; that is, interventions which are constant over time are defined. Ifabar is a matrix (of size 'number interventions' x 'time points'), then each row of the length ofAnodes refers to a particular time-varying intervention strategy. The natural intervention can be picked by settingabar='natural'.

The fitted outcome and confounder models are based on generalized additive models (GAMs) as implemented in themgcv package. Model families are picked automatically and reported in the output ifverbose=TRUE (see manual for modifications, though they hardly ever make sense). The model formulas are standard GLMs or GAMs (with penalized splines for continuous covariates), conditional on the past, unless specific formulae are given. It is recommended to use customized formulae to reduce the risk of model mis-specification and to ensure that the models make sense (e.g., not too many splines are used when this is computationally not meaningful). This can be best facilitated by using objects generated throughmake.model.formulas, followed bymodel.formulas.update and/ormodel.update (see examples for those functions).

For survival settings, it is required that i)survivalY=TRUE and ii) after a Cnode/Ynode is 1, every variable thereafter is set toNA. See manual for an example. By default, the package intervenes on Cnodes, i.e. calculates counterfactual outcomes under no censoring.

Ifcalc.support=TRUE, conditional and crude support measures (i.e., diagnostics) are calculated as described in Section 3.4 of Schomaker et al. (2023). Another useful diagnostic for multiple time points is the natural course scenario, which can be evaluated underabar='natural' andcbar='natural'.

To parallelize computations automatically, it is sufficient to setncores>1, as appropriate. To make estimates under parallelization reproducible, use theseed argument. To watch the progress of parallelized computations, set a path in theprog argument: then, a text file reports on the progress, which is particularly useful if lengthy bootstrapping computations are required.

Value

Returns an object of ofclass ‘gformula’:

Author(s)

Michael Schomaker

See Also

plot.gformula for plotting results as (causal) dose response curves,custom.measure for evaluating custom estimands andmi.boot for usinggformula on multiply imputed data.

Examples

## Not run: data(EFV)est <- gformula(X=EFV,                      Lnodes  = c("adherence.1","weight.1",                                  "adherence.2","weight.2",                                  "adherence.3","weight.3",                                  "adherence.4","weight.4"                        ),                        Ynodes  = c("VL.0","VL.1","VL.2","VL.3","VL.4"),                        Anodes  = c("efv.0","efv.1","efv.2","efv.3","efv.4"),                        abar=seq(0,10,1))est## End(Not run)

Compose appropriate model formulas

Description

Function that generates generic model formulas for Y-/L-/A- and Cnodes, according to time ordering and to be used ingformula ormodel.formulas.update.

Usage

make.model.formulas(X, Ynodes = NULL, Lnodes = NULL, Cnodes = NULL, Anodes = NULL,                    survival = FALSE, evaluate = FALSE)

Arguments

Details

This is a helper function to generate model formulas for Y-/L-/A- and Cnodes, according to the time ordering: i.e. to generate GLM/GAM model formulas for the respective nodes given allpast variables. In survival settings, past censoring and outcome nodes are omitted from the formulae. If censoring is present without a survival setting (e.g. Cnodes describe drop-outs and Y is a continuous outcome), then survival should be set as FALSE.

Value

Returns a named list:

See Also

The generated generic model formulas can be updated manually withmodel.update or in an automated manner with screening usingmodel.formulas.update.

Examples

data(EFV)m <- make.model.formulas(X=EFV,                         Lnodes  = c("adherence.1","weight.1",                                     "adherence.2","weight.2",                                     "adherence.3","weight.3",                                     "adherence.4","weight.4"                                    ),                         Ynodes  = c("VL.0","VL.1","VL.2","VL.3","VL.4"),                         Anodes  = c("efv.0","efv.1","efv.2","efv.3","efv.4"),                         evaluate=FALSE) # set TRUE to see fitted models                         m$model.names # all models potentially relevant for gformula(), given full past

Obtaining estimates from multiply imputed data

Description

Combinesgformula estimates obtained from multiple imputed data sets according to theMI Boot andMI Boot pooled methods decribed in Schomaker and Heumann (2018, see reference section below)

Usage

mi.boot(x, fun, cond=NULL, pooled=FALSE, ...)

Arguments

Value

An object of classgformula. Seegformula for details.

Author(s)

Michael Schomaker

References

Schomaker, M., Heumann, C. (2018)Bootstrap inference when using multiple imputation,Statistics in Medicine, 37:2252-2266

Examples

data(EFV)# suppose the following subsets were actually multiply imputed data (M=2)EFV_1 <- EFV[1:2500,]EFV_2 <- EFV[2501:5000,]# first: conduct analysis on each imputed data set. Set ret=T.m1 <- gformula(X=EFV_1,               Lnodes  = c("adherence.1","weight.1",                           "adherence.2","weight.2",                           "adherence.3","weight.3",                           "adherence.4","weight.4"               ),               Ynodes  = c("VL.0","VL.1","VL.2","VL.3","VL.4"),               Anodes  = c("efv.0","efv.1","efv.2","efv.3","efv.4"),               abar=seq(0,5,1), verbose=FALSE, ret=TRUE        )m2 <- gformula(X=EFV_2,               Lnodes  = c("adherence.1","weight.1",                           "adherence.2","weight.2",                           "adherence.3","weight.3",                           "adherence.4","weight.4"               ),               Ynodes  = c("VL.0","VL.1","VL.2","VL.3","VL.4"),               Anodes  = c("efv.0","efv.1","efv.2","efv.3","efv.4"),               abar=seq(0,5,1), verbose=FALSE, ret=TRUE)# second combine resultsm_imp <- mi.boot(list(m1,m2), mean) # uses MI rules & returns 'gformula' objectplot(m_imp)# custom estimand: evaluate probability of suppression (Y=0), among femalesm_imp2 <- mi.boot(list(m1,m2), prop, categ=0, cond="sex==1")plot(m_imp2)

Update model formulas based on variable screening

Description

Wrapper function to facilitate variable screening on all models generated throughmake.model.formulas and return updated formulas in the appropriate format forgformula.

Usage

model.formulas.update(formulas, X, screening = screen.glmnet.cramer,                      with.s = FALSE, by= NA, ...)

Arguments

Details

The default screening algorithm uses LASSO for variable screening (and Cramer's V for the categorized version of all variables if LASSO fails). It is possible to provide user-specific screening algorithms. User-specific algorithms should take the data as first argument,one model formula (i.e. one entry of the list inmodel.formulas) as second argument and return a vector of strings, containing the variable names that remain after screening. Another screening algorithm available in the package isscreen.cramersv, which categorizes all variables, calculates their association with the outcome based on Cramer'sV and selects the 4 variables with strongest associations (can be changed with optionnscreen).The manual provides more information.

The fitted models of the updated models can be evaluated withfit.updated.formulas.

Value

A list of length 4 containing the updated model formulas:

See Also

make.model.formulas,model.update,fit.updated.formulas

Examples

data(EFV)# first: generate generic model formulasm <- make.model.formulas(X=EFV,                         Lnodes  = c("adherence.1","weight.1",                                     "adherence.2","weight.2",                                     "adherence.3","weight.3",                                     "adherence.4","weight.4"                                    ),                         Ynodes  = c("VL.0","VL.1","VL.2","VL.3","VL.4"),                         Anodes  = c("efv.0","efv.1","efv.2","efv.3","efv.4"),                         evaluate=FALSE)                          # second: update these model formulas based on variable screening with LASSOglmnet.formulas <-  model.formulas.update(m$model.names, EFV)glmnet.formulas # third: use these models for estimationest <- gformula(X=EFV,                Lnodes  = c("adherence.1","weight.1",                            "adherence.2","weight.2",                            "adherence.3","weight.3",                            "adherence.4","weight.4"                ),                Ynodes  = c("VL.0","VL.1","VL.2","VL.3","VL.4"),                Anodes  = c("efv.0","efv.1","efv.2","efv.3","efv.4"),                Yform=glmnet.formulas$Ynames, Lform=glmnet.formulas$Lnames,                abar=seq(0,2,1))est

Update GAM models

Description

A wrapper to simplify the update of GAM models

Usage

model.update(gam.object, form)

Arguments

Details

Thegam object needs to be fitted with the option control=list(keepData=T), otherwise the function can not access the data that is needed to update the model fit. Note that bothfit.updated.formulas andmake.model.formulas with optionevaluate=T produce results that are based on this option.

Value

An object ofclass ‘gam’, ‘glm’ and ‘lm’.

Examples

data(EFV)m <- make.model.formulas(X=EFV,                         Lnodes  = c("adherence.1","weight.1",                                     "adherence.2","weight.2",                                     "adherence.3","weight.3",                                     "adherence.4","weight.4"                         ),                         Ynodes  = c("VL.0","VL.1","VL.2","VL.3","VL.4"),                         Anodes  = c("efv.0","efv.1","efv.2","efv.3","efv.4"),                         evaluate=TRUE) # set TRUE for model.update()# update first confounder model of weight manuallymodel.update(m$fitted.models$fitted.L$m_weight.1, .~s(weight.0, by=sex))# manual update of model formulam$model.names$Lnames[2] <- "weight.1 ~ s(weight.0, by=sex)"

Plot dose-response curves

Description

Function to plot dose-response curves based on results returned fromgformula

Usage

## S3 method for class 'gformula'plot(x, msm.method = c("line","loess", "gam", "none"),                         CI = FALSE, time.points = NULL,                         cols = NULL, weight = NULL, xaxis=NULL,                        variable = "psi", difference = FALSE, ...)

Arguments

Details

Time points and variable names should be specified according to the labeling of the results table returned bygformula.

Value

Draws an object ofclass ‘ggplot’.

Examples

data(EFV)est <- gformula(X=EFV,                Lnodes  = c("adherence.1","weight.1",                            "adherence.2","weight.2",                            "adherence.3","weight.3",                            "adherence.4","weight.4"                ),                Ynodes  = c("VL.0","VL.1","VL.2","VL.3","VL.4"),                Anodes  = c("efv.0","efv.1","efv.2","efv.3","efv.4"),                abar=seq(0,10,1))plot(est)plot(est, time.points=c(1,5))

Printing acontrastResult object

Description

Generic S3 function to print the results of acontrastResult object created withcontrast.

Usage

## S3 method for class 'contrastResult'print(x, ...)

Arguments

Value

No return value, called for printing the results of an object of classcontrastResult.

Printing the results of agformula object

Description

Generic S3 function to print the results of an object created withgformula.

Usage

## S3 method for class 'gformula'print(x, ...)

Arguments

Value

No return value, called for printing the results of an object of classclass ‘gformula’.

Helper function to calculate proportions of categorical variables

Description

To be used asfun argument incustom.measure.

Arguments

Value

A scalar numeric value of the proportion ofcateg values invec