--- title: "Time-Varying Causal Excursion Effect Mediation in MRT: Continuous Distal Outcomes" author: "Tianchen Qian (t.qian@uci.edu)" date: "`r Sys.Date()`" # output: # prettydoc::html_pretty: # theme: architect # or choose from: cayman, tactile, architect, leonids, hpstr # highlight: github # syntax highlighting style link-citations: yes bibliography: mhealth-ref.bib csl: biostatistics.csl vignette: > %\VignetteIndexEntry{Time-Varying Causal Excursion Effect Mediation in MRT: Continuous Distal Outcomes} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) options(digits = 4) library(MRTAnalysis) ``` The `MRTAnalysis` package now supports analysis of mediated causal excursion effect of time-varying treatments, time-varying mediators, and a continuous **distal outcome** in micro-randomized trials (MRTs), using the function `mcee()`. Distal outcomes are measured once at the end of the study (e.g., weight loss, cognitive score), in contrast to **proximal outcomes** which are repeatedly measured after each treatment decision point. Time-varying mediators can often be the proximal outcomes in MRTs. ## 1. Introduction **Micro-randomized trials (MRTs)** are experimental designs to evaluate and optimize quickly-adapting interventions, where each participant is randomized among treatment options many times throughout the study. This vignette shows how to use the `mcee` (stands for "mediated causal excursion effect") function to assess **how treatment effects on a distal outcome are mediated through intermediate variables** (e.g., the activity prompt -> near-term activity -> long-term cardiovascular health pathway). This package implements estimation of the **natural direct excursion effect (NDEE)** and the **natural indirect excursion effect (NIEE)** in MRTs and related longitudinal designs: - The **NIEE** captures the pathway from treatment at decision time \(t\) (\(A_t\)) to the *immediate mediator* (\(M_t\)), and then from \(M_t\) to the distal outcome \(Y\). - The **NDEE** lumps together *all other pathways* from \(A_t\) to \(Y\) that do not operate through the immediate mediator. The modeling of the direct and indirect effects is allowed to vary with the **decision point index**, specified through the argument `time_varying_effect_form`. For example, one may model NIEE and NDEE as constants over time (using `~1`) or as quadratic functions of the decision point (using `~ dp + I(dp^2)`). The parameter that parameterizes NIEE is denoted by $\beta$ and the one for NDEE is $\alpha$. At present, the software does *not* support mediation effects conditional on baseline or other time-varying covariates. For full details of the estimands, see @qian2025dynamic. ### Estimation framework For applied users: if analyzing an MRT, you typically know the randomization probability, and need to provide a `control_formula_with_mediator` describing covariates and mediators to include in nuisance models. By default, all nuisance regressions are estimated with GLMs, but you can substitute other machine learning methods. For statistically inclined readers: the estimator is **multiply robust** and generalizes the semiparametric estimators of [Tchetgen Tchetgen & Shpitser (2012)](https://doi.org/10.1214/12-AOS990) to longitudinal settings. It's a two-stage estimation process, where the first stage fits five nuisance functions and the second stage estimates the parameterized NIEE and NDEE. The five nuisance functions are (see @qian2025dynamic for the precise definition): - $p$: treatment assignment mechanism (in MRTs this is *known* and should be specified as such --- see later). - $q$: treatment assignment further conditional on the mediator. - $\eta$: outcome regression. - $\mu$: outcome regression further conditional on the mediator. - $\nu$: a "cross-world" outcome regression. ### Three user-facing wrappers The package provides three entry points depending on how much control you want: - **`mcee`**: streamlined interface for MRTs with known randomization probabilities. Recommended default for most MRT analyses. - **`mcee_general`**: more flexible, with explicit configuration of nuisance models. Useful for observational studies or MRTs when experimenting with different learners. - **`mcee_userfit_nuisance`**: accepts user-supplied nuisance predictions (from external ML fits). The remainder of this vignette will first focus on `mcee` for a quick introduction, then show how to use the more flexible wrappers. --- ## 2. Installation You can install the **MRTAnalysis** package from CRAN: ```{r install_package} # install.packages("MRTAnalysis") ``` --- ## 3. Quick Start Example We illustrate the most basic usage of `mcee` on a small simulated dataset. Note that the dataset is in **long format**, one row per subject-by-decision point. The distal outcome is repeated on each row as a constant within subject. ```{r minimum_example} set.seed(123) # Simulate a toy dataset n <- 20 T_val <- 5 id <- rep(1:n, each = T_val) dp <- rep(1:T_val, times = n) A <- rbinom(n * T_val, 1, 0.5) M <- rbinom(n * T_val, 1, plogis(-0.2 + 0.3 * A + 0.1 * dp)) Y <- ave(0.5 * A + 0.7 * M + 0.2 * dp + rnorm(n * T_val), id) # constant within id dat <- data.frame(id, dp, A, M, Y) # Minimal mcee call fit <- mcee( data = dat, id = "id", dp = "dp", outcome = "Y", treatment = "A", mediator = "M", time_varying_effect_form = ~1, # constant-over-time NDEE and NIEE control_formula_with_mediator = ~ dp + M, # nuisance adjustment control_reg_method = "glm", # default method rand_prob = 0.5, # known randomization prob verbose = FALSE ) summary(fit) ``` The `summary` output provides estimates of Natural Direct Excursion Effect (NDEE; $\alpha$) and Natural Indirect Excursion Effect (NIEE; $\beta$), with standard errors, 95\% confidence intervals, and p-values. The only row in the output `Natural Direct Excursion Effect (alpha)` and `Natural Indirect Excursion Effect (beta)` is named `Intercept`, indicating that these effects are modeled as constants over time (like intercepts in the effect models). In particular, the estimated NDEE is 0.17 (95\% CI: -0.08, 0.42; p-value: 0.17) and the estimated NIEE is 0.025 (95\% CI: -0.002, 0.054; p-value: 0.06). The confidence intervals and p-values are based on t-quantiles. ## 4. A More Complex Data Example This package ships a small example dataset, `data_time_varying_mediator_distal_outcome`. It is in **long format**, one row per subject-by-decision point. ### Columns (Variables) - `ID`: subject identifier (use the actual column name present in your data) - `t_val`: decision point index (must be strictly increasing within each subject) - `I`: availability (binary in `{0,1}`); if omitted, the analysis assumes all rows are available (which is not the case for this data example) - `p_A`: randomization probability, i.e., probability of $A_{it} = 1$ at a given decision point conditional on the history information - `A`: binary treatment, coded in `{0,1}` - `M`: mediator (continuous or binary are both allowed) - `Y`: **distal** outcome — constant within each subject (the same value repeated across the subject's rows) - Additional time-varying covariates for adjustment (e.g., `X`, `A_prev`, `M_prev`, etc.) ### Peek at the included dataset ```{r data_example} data(data_time_varying_mediator_distal_outcome) dat <- data_time_varying_mediator_distal_outcome dplyr::glimpse(dat) dplyr::count(dat, id, name = "Ti") |> dplyr::summarise(mean_T = mean(Ti), min_T = min(Ti), max_T = max(Ti)) # Delete some decision points for certain individuals to mimic scenarios # where not all individuals have the same number of decision points. dat <- dat[!((dat$id == 1 & dat$dp == 10) | (dat$id == 2 & dat$dp %in% c(9, 10))), ] dplyr::count(dat, id, name = "Ti") |> dplyr::summarise(mean_T = mean(Ti), min_T = min(Ti), max_T = max(Ti)) ``` --- ## 5. Detailed Walkthrough of `mcee` We focus first on the streamlined MRT workflow with known randomization probability. Throughout, we'll use the included dataset. ### 5.1 Basic usage (GLM; constant effects) ```{r basid_usage} set.seed(1) fit1 <- mcee( data = dat, id = "id", dp = "dp", outcome = "Y", treatment = "A", mediator = "M", availability = "I", rand_prob = "p_A", time_varying_effect_form = ~1, # NDEE and NIEE are constant over time control_formula_with_mediator = ~ dp + M + X, # covariate adjustment control_reg_method = "glm", # nuisance learners for q, eta, mu, nu verbose = FALSE ) summary(fit1) ``` - `time_varying_effect_form`: defines the basis f(t) used to model alpha (direct) and beta (indirect). - `control_formula_with_mediator`: adjustment set for nuisance models; **include the mediator** here (the wrapper will automatically remove it from the nuisance function models that should not include it). This can include `s()` terms if `control_reg_method = "gam"`. ### 5.2 Time-varying effects (linear in `dp`) ```{r time_varying_effect} fit2 <- mcee( data = dat, id = "id", dp = "dp", outcome = "Y", treatment = "A", mediator = "M", availability = "I", rand_prob = "p_A", time_varying_effect_form = ~dp, # NDEE, NIEE vary linearly in time control_formula_with_mediator = ~ dp + M + X, control_reg_method = "glm", verbose = FALSE ) summary(fit2) # rows now labeled (Intercept) and dp ``` Interpretation: the rows for `dp` report how the natural direct/indirect excursion effects change per unit increase in decision point. > **Tip.** If your dataset already contains a basis of time (e.g., spline columns), you can reference them in `time_varying_effect_form`. The package will warn when variables other than `dp` appear there, but this is allowed for precomputed time bases. ### 5.3 Other learners for nuisance fitting You can switch the learner used for fitting the nuisance functions via `control_reg_method`: - Generalized Additive Models: `"gam"` - Random forest (randomForest): `"rf"` - Ranger (fast RF): `"ranger"` - SuperLearner: `"sl"` ```{r other_learners_gam} # Example: GAM (generalized additive model) set.seed(2) fit3 <- mcee( data = dat, id = "id", dp = "dp", outcome = "Y", treatment = "A", mediator = "M", availability = "I", rand_prob = "p_A", time_varying_effect_form = ~dp, control_formula_with_mediator = ~ s(dp) + s(M) + s(X), # spline formula for mgcv::gam control_reg_method = "gam", verbose = FALSE ) summary(fit3) ``` ### 5.4 Estimating effects at specific decision points If interested in estimating direct and indirect effects of intervention at a specific decision point (or a few decision points), you can use the `specific_dp_only` argument. For example, to estimate effects at the first two decision points (1 and 2) only, you can either: - Supply `specific_dp_only = c(1, 2)` in `mcee`, **or** - Provide `weight_per_row` that puts nonzero mass only at those rows (the wrappers normalize weights within subject). ```{r specific_dp_only} fit4 <- mcee( data = dat, id = "id", dp = "dp", outcome = "Y", treatment = "A", mediator = "M", availability = "I", rand_prob = "p_A", time_varying_effect_form = ~1, control_formula_with_mediator = ~ dp + M + X, control_reg_method = "glm", specific_dp_only = c(1, 2), verbose = FALSE ) summary(fit4) ``` The effect estimates are now an average over decision points 1 and 2. In the next section of the vignette, we'll show how to use the more flexible wrappers, `mcee_general` and `mcee_userfit_nuisance`. ## 6. Testing for Linear Combinations of Coefficients Every wrapper in this package returns an object of class `mcee_fit`. This object is a list with several components. ### 6.1 Estimates `$mcee_fit` - `alpha_hat`, `beta_hat`: estimated coefficients for the **natural direct excursion effect** ($\alpha$) and **natural indirect excursion effect** ($\beta$). - `alpha_se`, `beta_se`: standard errors. - `varcov`: estimated covariance matrix for both $\alpha$ and $\beta$ coefficients (ordered as $(\alpha^T, \beta^T)^T$). - `alpha_varcov`, `beta_varcov`: estimated covariance matrices for $\alpha$ and $\beta$ separately. ```{r fit1_mcee_fit} fit1$mcee_fit ``` ### 6.2 Linear combinations The `lincomb_joint`, `lincomb_alpha`, and `lincomb_beta` arguments to `summary()` lets you compute linear combinations of the estimated coefficients. For example: #### Example 1: Difference $\alpha - \beta$ for constant‑effect model (`fit1`) ```{r lincomb_example1} # difference between direct and indirect excursion effects summary(fit1, lincomb_joint = c(1, -1)) ``` This requests the linear combination $\alpha - \beta$, with standard error and t-statistic computed from the covariance matrix. #### Example 2: Effects at decision point 10 for linear‑trend model (`fit2`) Suppose you fit a model with a linear time basis, e.g.: ```{r lincomb_example2} fit2 <- mcee( data = dat, id = "id", dp = "dp", outcome = "Y", treatment = "A", mediator = "M", availability = "I", rand_prob = "p_A", time_varying_effect_form = ~dp, control_formula_with_mediator = ~ dp + M + X, control_reg_method = "glm", verbose = FALSE ) ``` To get the effect at decision point 10 (the last decision point), you need the intercept plus 9 times the slope: ```{r lincomb_example2_continued} summary(fit2, lincomb_alpha = c(1, 9), lincomb_beta = c(1, 9)) ``` If you prefer a joint contrast (e.g., alpha(t=10) − beta(t=10)), stack the two contrasts into a single row over `c(alpha, beta)`: ```{r lincomb_example2_joint} L_joint_t10 <- matrix(c( 1, 9, # alpha part -1, -9 # beta part ), nrow = 1) summary(fit2, lincomb_joint = L_joint_t10) ``` ### 6.3 Inspecting nuisance models For diagnostics, the returned object also includes: - `$nuisance_models`: the actual model objects used for fitting each nuisance parameter (`p`, `q`, `eta`, `mu`, `nu`). - `$nuisance_fitted`: the fitted values for each nuisance function (`p1`, `q1`, `eta1`, `eta0`, `mu1`, `mu0`, `nu1`, `nu0`). One can use `summary(..., show_nuisance = TRUE)` to get a compact summary of the nuisance fits. ```{r inspect_nuisance} summary(fit1, show_nuisance = TRUE) head(fit1$nuisance_fitted$mu1) ``` > **Tip.** Because `varcov` is the joint covariance of `c(alpha, beta)`, you can build any custom Wald test by forming your own contrast matrix `L` and using `L %*% c(alpha, beta)` with variance `L %*% varcov %*% t(L)`. This is equivalent to using `lincomb_joint = L` in `summary()`. ## 7. Other Wrappers: `mcee_general` and `mcee_userfit_nuisance` While `mcee` covers the common MRT use case (known randomization, one control formula for all nuisance functions), two additional wrappers provide more flexibility in controlling how the nuisance functions are fitted: ### 7.1 `mcee_general` This interface lets you configure each nuisance function separately. - **Configuration objects** specify: - `known`: supply a fixed vector/scalar (commonly for treatment mechanism `p` in MRT). - `formula`: regression formula, which can include `s()` terms if using `method = "gam"`. - `method`: learner (`glm`, `gam`, `lm`, `rf`, `ranger`, `sl`). - `family`: optional; defaults inferred from nuisance type (`binomial` for `p` and `q`; `gaussian` for `eta`, `mu`, `nu`). Helper functions can be used to build these objects: `mcee_config_maker`, `mcee_config_known`, `mcee_config_glm`, `mcee_config_gam`, `mcee_config_lm`, `mcee_config_rf`, `mcee_config_ranger`, `mcee_config_sl`, `mcee_config_sl_user`. #### Example: flexible nuisance configs ```{r } # Families (binomial vs Gaussian) are chosen automatically when omitted; for `p` and `q` the default is binomial, for the outcome regressions it is Gaussian. cfg <- list( p = mcee_config_known("p", dat$p_A), # known randomization prob in MRT q = mcee_config_glm("q", ~ dp + X + M), eta = mcee_config_glm("eta", ~ dp + X), mu = mcee_config_glm("mu", ~ dp + X + M), nu = mcee_config_glm("nu", ~ dp + X) ) fit_gen <- mcee_general( data = dat, id = "id", dp = "dp", outcome = "Y", treatment = "A", mediator = "M", availability = "I", time_varying_effect_form = ~dp, config_p = cfg$p, config_q = cfg$q, config_eta = cfg$eta, config_mu = cfg$mu, config_nu = cfg$nu, verbose = FALSE ) summary(fit_gen) ``` > **Tip.** This interface is particularly useful for **observational studies** where the treatment mechanism `p` is unknown and must be estimated. You can set `config_p` to a regression formula and method of your choice. ### 7.2 `mcee_userfit_nuisance` Sometimes you fit nuisance models externally, perhaps with custom ML workflows. In that case, you can bypass Stage‑1 nuisance model fitting and supply predictions directly. - Required inputs: vectors `p1`, `q1`, `eta1`, `eta0`, `mu1`, `mu0`, `nu1`, `nu0`, each aligned with the rows of your data. - `p1`: $P(A_t = I_t | H_t)$ - `q1`: $P(A_t = I_t | H_t, M_t)$ - `eta1`: $E(Y | H_t, A_t = I_t)$ - `eta0`: $E(Y | H_t, A_t = 0)$ - `mu1`: $E(Y | H_t, A_t = I_t, M_t)$ - `mu0`: $E(Y | H_t, A_t = 0, M_t)$ - `nu1`: $E[E(Y | H_t, A_t = I_t, M_t) | H_t, A_t = 0]$ - `nu0`: $E[E(Y | H_t, A_t = 0, M_t) | H_t, A_t = I_t]$ - The wrapper will enforce `p1 = 1` and `q1 = 1` at `I == 0`. #### Example: manual nuisance fits with GLM ```{r } # Fit nuisance regressions manually: p, q, eta, mu, nu p1_hat <- dat$p_A # known randomization prob in MRT p1_hat[dat$I == 0] <- 1 # manually set this to avoid wrapper message q1_hat <- predict(glm(A ~ dp + X + M, family = binomial(), data = dat[dat$I == 1, ]), type = "response", newdata = dat ) q1_hat[dat$I == 0] <- 1 # manually set this to avoid wrapper message eta1_hat <- predict(lm(Y ~ dp + X, data = dat[dat$A == dat$I, ]), newdata = dat) eta0_hat <- predict(lm(Y ~ dp + X, data = dat[dat$A == 0, ]), newdata = dat) mu1_hat <- predict(lm(Y ~ dp + X + M, data = dat[dat$A == dat$I, ]), newdata = dat) mu0_hat <- predict(lm(Y ~ dp + X + M, data = dat[dat$A == 0, ]), newdata = dat) nu1_hat <- predict(lm(mu1 ~ dp + X, data = cbind(dat, mu1 = mu1_hat)[dat$A == 0, ]), newdata = dat) nu0_hat <- predict(lm(mu0 ~ dp + X, data = cbind(dat, mu0 = mu0_hat)[dat$A == dat$I, ]), newdata = dat) fit_usr <- mcee_userfit_nuisance( data = dat, id = "id", dp = "dp", outcome = "Y", treatment = "A", mediator = "M", availability = "I", time_varying_effect_form = ~dp, p1 = p1_hat, q1 = q1_hat, eta1 = eta1_hat, eta0 = eta0_hat, mu1 = mu1_hat, mu0 = mu0_hat, nu1 = nu1_hat, nu0 = nu0_hat, verbose = FALSE ) summary(fit_usr) ``` - Particularly useful when applying **custom ML pipelines** not wrapped by `mcee`. - For MRTs, you can still supply `p1` as a known constant vector. ### 7.3 Comparing results Different wrappers will produce the exact same result if the same nuisance models are used. Here we compare `fit2` (from `mcee`) and `fit_gen` (from `mcee_general`) which use the same GLM nuisance models and the same control variables. ```{r check_equal_1} fit2$mcee_fit[c("alpha_hat", "beta_hat", "alpha_se", "beta_se")] fit_gen$mcee_fit[c("alpha_hat", "beta_hat", "alpha_se", "beta_se")] all.equal(fit2$mcee_fit, fit_gen$mcee_fit, tolerance = 1e-6) ``` ```{r check_equal_2} all.equal(fit2$mcee_fit, fit_usr$mcee_fit, tolerance = 1e-6) ``` ## 8. Best Practices - **Start with `mcee`.** For analyzing MRTs with known randomization probabilities, this is the simplest choice. - **Moderator formula:** Use only functions of `dp` in `time_varying_effect_form` (or simply set it to `~1`). Do not include baseline or time-varying covariates here. If you've precomputed spline or polynomial bases of `dp`, referencing those columns is allowed. - **Include the mediator** in `control_formula_with_mediator` (or in `config_mu`/`config_q` for `mcee_general`). - **Check your data structure.** Rows must be grouped by subject, with strictly increasing `dp`. Outcomes must be constant within subjects. No missing values are allowed. - **Check ML learners fits:** While you can specify flexible learners (GAM, RF, Ranger, SuperLearner), always check summary outputs and nuisance fits for plausibility. --- ## 9. Conclusion The `mcee` package provides a multiply robust estimator of **natural direct and indirect excursion effects** in micro-randomized trials and observational longitudinal studies with time-varying treatments, time-varying mediators and a distal outcome. Three wrapper functions serve different needs: - `mcee`: the streamlined MRT workflow with known randomization. - `mcee_general`: full control over nuisance model specifications. - `mcee_userfit_nuisance`: for injecting externally fitted nuisance parameters. Together, these tools allow both applied researchers and methodologists to estimate mediated excursion effects with transparent assumptions and flexible modeling options. ## Appendix: Estimating Equations and Asymptotic Variance This appendix documents how the core calculation happens in `.mcee_core_rows()`, which is a detailed derivation of Eq. (12) and the asymptotic variance of Theorem 4 in @qian2025dynamic. The proposed estimating function is \[ \psi(\gamma; \zeta) := \sum_{t=1}^T \omega(t) \begin{bmatrix} \{ \phi_t^{10}(p_t, q_t, \mu_t, \nu_t) - \phi_t^{00}(p_t, \eta_t) - f(t)^\top \alpha \} f(t) \\ \{ \phi_t^{11}(p_t, \eta_t) - \phi_t^{10}(p_t, q_t, \mu_t, \nu_t) - f(t)^\top \beta \} f(t) \end{bmatrix}. \] The estimating equation (EE), omitting the nuisance parameters and allowing different \(T_i\) for each subject \(i\), is \[ \frac{1}{n}\sum_{i=1}^n \sum_{t=1}^{T_i} \omega(i,t) \begin{bmatrix} \{ \phi_{it}^{10} - \phi_{it}^{00} - f(t)^\top \alpha \} f(t) \\ \{ \phi_{it}^{11} - \phi_{it}^{10} - f(t)^\top \beta \} f(t) \end{bmatrix} = 0. \] The estimators are computed by \[ \hat\alpha = \left[ \frac{1}{n}\sum_{i=1}^n \sum_{t=1}^{T_i} \omega(i,t) f(t)f(t)^\top \right]^{-1} \left[ \frac{1}{n}\sum_{i=1}^n \sum_{t=1}^{T_i} \omega(i,t) \{ \phi_{it}^{10} - \phi_{it}^{00} \} f(t) \right], \] \[ \hat\beta = \left[ \frac{1}{n}\sum_{i=1}^n \sum_{t=1}^{T_i} \omega(i,t) f(t)f(t)^\top \right]^{-1} \left[ \frac{1}{n}\sum_{i=1}^n \sum_{t=1}^{T_i} \omega(i,t) \{ \phi_{it}^{11} - \phi_{it}^{10} \} f(t) \right]. \] The asymptotic variance of \((\alpha^\top, \beta^\top)^\top\) can be estimated by \[ \text{Var}(\hat\gamma) \approx \text{Bread}^{-1} \, \text{Meat} \, (\text{Bread}^{-1})^\top, \] where \[ \text{Bread} = \frac{1}{n}\sum_{i=1}^n \sum_{t=1}^{T_i} \omega(i,t) \begin{bmatrix} f(t)f(t)^\top & 0 \\ 0 & f(t)f(t)^\top \end{bmatrix}, \] \[ \text{Meat} = \frac{1}{n}\sum_{i=1}^n \left( \sum_{t=1}^{T_i} \omega(i,t) \begin{bmatrix} \{ \phi_{it}^{10} - \phi_{it}^{00} - f(t)^\top \alpha \} f(t) \\ \{ \phi_{it}^{11} - \phi_{it}^{10} - f(t)^\top \beta \} f(t) \end{bmatrix} \right)^{\otimes 2}. \] ## Reference