Tutorials › Marginal Structural Models › Simulation Workshop On Feedback-Induced Regression Bias

Time-Varying Treatments and Confounder Feedback · Tutorial 27 of 501

Simulation Workshop On Feedback-Induced Regression Bias

Build and interpret a feedback simulation, derive its regression targets, and see why adjusting for a treatment-affected confounder changes the meaning of the treatment coefficients.

Intermediate 12 min read

What You'll Learn

  • Specify a two-visit data-generating process with treatment-confounder feedback.
  • Derive the population coefficients from baseline-adjusted and time-varying-adjusted regressions.
  • Distinguish those coefficients from the outcome contrast under an assigned treatment history.
  • Write a reproducible simulation and check whether its estimates approach the derived targets.
  • Diagnose misleading regression comparisons and identify what they cannot establish.

Why Simulate Treatment-Confounder Feedback?

When treatment is repeated, an earlier treatment can change a later covariate, and that covariate can influence both later treatment and the outcome. This is treatment-confounder feedback. A simulation makes the resulting regression problem concrete: it lets us compare what ordinary regressions estimate with the contrast under a specified treatment history. The comparison is useful only if those quantities are kept distinct.

We will use continuous treatments and a linear data-generating process so that the population regression coefficients can be derived exactly. This is a teaching example, not a claim that clinical treatment is always continuous or that the linear model is generally adequate.

Definition: In a two-visit process, \(A_0\) is treatment at the first visit, \(L_1\) is a covariate measured after \(A_0\) and before the next treatment, \(A_1\) is treatment at the second visit, and \(Y\) is the final outcome. Feedback occurs here because \(A_0\) affects \(L_1\), while \(L_1\) affects both \(A_1\) and \(Y\).

For each individual, generate independent errors \(U_L,U_A,U_Y\), each distributed as \(N(0,1)\), and generate \(A_0\sim N(0,1)\), independently of those errors. Use the structural equations

$$ L_1=0.8A_0+U_L,\qquad A_1=0.6L_1+U_A,\qquad Y=A_0+2A_1+1.5L_1+U_Y. $$

The equations specify the simulation completely. In particular, the coefficient 0.8 represents the effect of earlier treatment on the later covariate; 0.6 represents the effect of that covariate on later treatment. The outcome depends on both treatments and on \(L_1\). The independent, mean-zero errors make the calculations below especially transparent.

Two Regressions, Two Different Questions

Fit these two population models, or their sample analogues:

  • Baseline-adjusted regression: regress \(Y\) on an intercept, \(A_0\), and \(A_1\).
  • Time-varying-adjusted regression: regress \(Y\) on an intercept, \(A_0\), \(A_1\), and \(L_1\).

The second regression conditions on a covariate that is both a consequence of \(A_0\) and a cause of \(A_1\). Its treatment coefficients are conditional regression coefficients. They are not automatically the mean outcome contrast between assigned treatment histories. The first regression does not condition on \(L_1\), but that omission leaves the \(A_1\) comparison confounded by the part of \(L_1\) that predicts \(A_1\).

Worked Example: Deriving the Baseline-Adjusted Targets

First derive the population coefficients for the regression of \(Y\) on \(A_0\) and \(A_1\), omitting \(L_1\). By the outcome equation, the fitted linear projection has the same coefficients as the projection of \(A_0+2A_1+1.5L_1+U_Y\) onto \(A_0,A_1\). Since \(U_Y\) is independent of both regressors, it contributes zero coefficients. The remaining task is to project \(L_1\) onto the two regressors.

Write \(A_1=0.48A_0+0.6U_L+U_A\). After removing its linear prediction from \(A_0\), the residual is \(0.6U_L+U_A\), with variance \(0.6^2+1=1.36\). Its covariance with the residual \(U_L\) in \(L_1=0.8A_0+U_L\) is \(0.6\). Therefore the coefficient of \(A_1\) in the projection of \(L_1\) on \(A_0,A_1\) is \(0.6/1.36=15/34\).

The coefficient of \(A_0\) in that projection is \(0.8-(15/34)(0.48)=10/17\). Thus the omitted \(L_1\) contributes \(1.5(10/17)=15/17\) to the \(A_0\) coefficient and \(1.5(15/34)=45/68\) to the \(A_1\) coefficient. Adding the coefficients already present in the outcome equation gives

$$ \beta_{A_0}^{\text{baseline}}=1+\frac{15}{17}=\frac{32}{17}\approx1.882,\qquad \beta_{A_1}^{\text{baseline}}=2+\frac{45}{68}=\frac{181}{68}\approx2.662. $$

These are population targets, not values that every finite simulated data set must reproduce exactly. The later-treatment coefficient exceeds 2 because the regression leaves out \(L_1\), which predicts both \(A_1\) and \(Y\).

Worked Example: Deriving the Time-Varying-Adjusted Targets

In the second regression, \(L_1\) is included. The outcome equation itself is linear in the included variables, and its error \(U_Y\) has conditional mean zero given \(A_0,A_1,L_1\). Consequently the population regression coefficients are exactly the structural coefficients:

$$ \beta_{A_0}^{\text{time-varying}}=1,\qquad \beta_{A_1}^{\text{time-varying}}=2,\qquad \beta_{L_1}^{\text{time-varying}}=1.5. $$

This argument requires the specified linear conditional mean and the stated independence of \(U_Y\) from the regressors. It does not say that these conditional coefficients are the effects of assigning treatment histories. In particular, conditioning on \(L_1\) holds fixed a variable that changes when \(A_0\) changes.

Compare the Regression Coefficients With a Treatment-History Contrast

To identify the contrast in this simulation, consider setting treatment to fixed values \(a_0\) and \(a_1\), rather than letting the treatment assignment equations choose those values. Under this intervention, the equation for \(L_1\) still operates, so its mean is \(0.8a_0\). Substituting that equation into the outcome equation gives the mean outcome under the assigned history:

$$ E[Y\mid\text{set }A_0=a_0,A_1=a_1] =a_0+2a_1+1.5(0.8a_0) =2.2a_0+2a_1. $$

The treatment-history contrast for increasing \(a_0\) by one unit while holding the assigned \(a_1\) fixed is therefore 2.2. The corresponding contrast for increasing \(a_1\) by one unit, with \(a_0\) fixed, is 2. The first contrast includes the pathway from \(A_0\) through \(L_1\) to \(Y\); the second does not change \(L_1\), which occurs before \(A_1\).

Result: In this data-generating process, the baseline-adjusted coefficients \((32/17,181/68)\) are not the treatment-history contrasts \((2.2,2)\). The time-varying-adjusted coefficients \((1,2)\) also are not the full treatment-history contrasts \((2.2,2)\): adjustment holds \(L_1\) fixed, removing the \(A_0\)-through-\(L_1\) pathway. These conclusions follow for this specified model and should not be generalized to every regression or data-generating process.

The substitution above is also a derivation of the contrast, rather than an appeal to a regression coefficient. Under the intervention, \(U_L\) and \(U_Y\) retain mean zero, so taking expectations of the substituted outcome removes those errors. The coefficient on \(a_0\) is \(1+(1.5)(0.8)=2.2\), and the coefficient on \(a_1\) remains 2. This distinction is the central check when comparing regression results: first state what contrast is wanted, then decide whether a fitted coefficient represents it.

Worked Example: A Reproducible Simulation Check

The following R code simulates \(n=5000\) individuals and fits both regressions. It reports the coefficient estimates; changing the random seed changes the finite-sample values. The targets derived above are the quantities the estimates should approach as the sample size increases under this data-generating process.

$$ \begin{aligned} &\texttt{set.seed(271)}\\ &\texttt{n <- 5000}\\ &\texttt{A0 <- rnorm(n)}\\ &\texttt{UL <- rnorm(n); UA <- rnorm(n); UY <- rnorm(n)}\\ &\texttt{L1 <- 0.8*A0 + UL}\\ &\texttt{A1 <- 0.6*L1 + UA}\\ &\texttt{Y <- A0 + 2*A1 + 1.5*L1 + UY}\\ &\texttt{coef(lm(Y \textasciitilde{} A0 + A1))}\\ &\texttt{coef(lm(Y \textasciitilde{} A0 + A1 + L1))} \end{aligned} $$

For a numerical benchmark, the baseline-adjusted treatment estimates should be near \(1.882\) for \(A_0\) and \(2.662\) for \(A_1\); the time-varying-adjusted estimates should be near \(1\) and \(2\). These are not promised as exact output from the particular finite sample. If a run gives noticeably different values, rerun with a larger \(n\) and check the simulated equations and model formulas before interpreting the difference.

A useful verification is to compare the sample associations implied by the generation steps: \(A_0\) should predict \(L_1\), and \(L_1\) should predict \(A_1\). Also verify that the fitted time-varying regression includes all three terms \(A_0,A_1,L_1\). Omitting a term, changing the order of generation, or reusing one error vector would simulate a different process.

How to Read the Comparison and Diagnose Problems

A simulation is informative because the data-generating equations give a known benchmark. The regression comparison then illustrates how estimates change when adjustment changes. It does not establish that either fitted coefficient is a causal effect in an observational study. In this constructed process, the equations and independent errors justify the calculations; in real data, those assumptions are not verified merely by fitting the models.

Common warning signs include estimates that fail to move toward the derived targets as \(n\) grows, a later-treatment coefficient that is interpreted without considering \(L_1\), or a baseline-treatment coefficient from the adjusted regression that is reported as the total effect despite conditioning on \(L_1\). For this continuous-treatment example, also inspect the range and variation of simulated treatments: an implementation that accidentally creates almost constant treatment values can make regression estimates unstable.

A frequent mistake is to argue that the time-varying-adjusted regression must be better because it includes more measured information. That reasoning confuses prediction or conditional association with a treatment-history contrast. Here its \(A_0\) coefficient is 1, while the assigned-history contrast is 2.2, because changing \(A_0\) changes \(L_1\). Conversely, omitting \(L_1\) does not automatically recover the total effect: the baseline-adjusted regression's \(A_0\) coefficient is about 1.882, not 2.2, and its \(A_1\) coefficient is about 2.662, not 2.

The appropriate comparison depends on the question. A regression coefficient describes a conditional linear relationship under its model specification. A treatment-history contrast compares outcomes under specified assignments. Methods such as the parametric g-formula and marginal structural models are designed to target marginal contrasts under their own assumptions; naming an alternative does not remove the need to define the estimand, justify assumptions, and diagnose model fit. This workshop isolates one reason ordinary adjustment can miss the desired contrast; it does not cover all sources of bias, nonlinearity, or treatment assignment mechanisms.

1
Write the data-generating process.
Specify which variables precede each treatment and which equations create feedback.
2
State the target contrast.
Describe the treatment history being compared before reading a coefficient as an effect.
3
Derive regression targets.
Use the model equations or population projection, rather than relying only on one random simulation run.
4
Simulate and check.
Compare fitted coefficients with their derived targets, then explain any difference in terms of sampling variation or a changed implementation.

Limits of the Workshop

The model is linear, errors are independent and Gaussian, treatments are continuous, and the relevant variables are measured without error. Those choices make the derivation exact but are not necessary features of feedback in general. With nonlinear outcome models, interactions, binary treatments, measurement error, or omitted common causes, the coefficients and their interpretations can differ substantially. A simulation demonstrates consequences of its assumptions; it cannot show that those assumptions describe a particular clinical study.

Check Your Understanding

Answer each question, then compare with the solutions.

  1. Under the specified process, what is the mean of \(L_1\) when \(A_0\) is set to \(a_0\)?
  2. Why is the baseline-adjusted \(A_1\) coefficient greater than the coefficient 2 in the outcome equation?
  3. What is the treatment-history contrast for increasing \(A_0\) by one unit while holding assigned \(A_1\) fixed? Show the calculation.
  4. Which coefficient from the time-varying-adjusted regression includes the pathway from \(A_0\) through \(L_1\) to \(Y\)? Explain.
  5. If a large simulation produces values far from the population targets, name two checks to perform before drawing a substantive conclusion.

Solutions. (1) The \(L_1\) equation gives \(E[L_1\mid\text{set }A_0=a_0]=0.8a_0\), since \(E[U_L]=0\). (2) The omitted \(L_1\) predicts \(A_1\) and affects \(Y\), so the baseline-adjusted regression assigns some of the \(L_1\)-outcome association to \(A_1\); its derived coefficient is \(2+45/68\approx2.662\). (3) Substitution gives \(1+1.5(0.8)=2.2\). (4) Neither treatment coefficient includes that full pathway when \(L_1\) is held fixed; the \(A_0\) conditional coefficient is 1. (5) Check that the generation equations and model formulas match the specification, and increase \(n\) or rerun with a new seed to assess sampling variation.

Clinical Biostats