*(Software available free to members. Register.)*

### Accompanying Articles

This R code accompanies the following journal articles:

Almirall, D., Ten Have, T., & Murphy, S. A. (2010). Structural nested mean models for assessing time-varying effect moderation.* Biometrics, 66*(1), 131-139. PMCID: PMC2875310 View article

Almirall, D., McCaffrey, D. F., Ramchand, R., & Murphy, S. A. (2011). Subgroups analysis when treatment and moderators are time-varying.* Prevention Science. *Advance online publication. * *doi:10.1007/s11121-011-0208-7 PMCID: PMC3135740 View abstract

### Overview

Robins' structural nested mean model (SNMM), as described in Almirall, et al. (2010), is useful when scientists are interested in understanding causal effect moderation in the presence of time-varying treatments (i.e., exposures) and time-varying moderators. This R code implements a simple version of the two-stage regression residuals estimator (Almirall et al., 2010) for the SNMM. It returns an estimate of the asymptotic variance-covariance matrix (and therefore standard errors) for the estimated SNMM parameters (stage 2), taking into account sampling error in the estimation of the stage-1 parameters. Read more about causal inference.

You will need to install R in order to use the twostage.snmm() function. R can be downloaded for free by visiting the R-Project Home Page.

### Download

This file contains the R function twostage.snmm(), which implements the 2-stage regression estimator. The header at the beginning of this file explains how to specify each of the function arguments. These include

- specifying the stage 1 regression models for the moderator (one for each moderator at each time point) given the past (modgivpast argument),
- specifying additional terms (coefnuisfunc argument), which are multiplied with the stage 1 residuals to form models for the (nuisance) error terms corresponding to each moderator, and
- specifying models for the intermediate causal effect functions (intermeff argument).

### A Worked Example

In the following, we describe how we fit the SNMM discussed in Almirall, McCaffrey, Ramchand, and Murphy (2011), using the 2-stage regression with residuals estimator as implemented in twostage.snmm(). We cannot provide you with the data for you to carry out this exact example, but after following the code and description below alongside the paper, you should be prepared to use twostage.snmm() in your own research. In the article, we described how to use a SNMM to examine two scientific questions: (1) "What is the impact of receiving treatment during months 1-3 versus not receiving treatment during months 1-3 (nor during months 7-9) on end-of-study environmental risk outcomes as a function of baseline severity?" and (2) "What is the impact of receiving treatment during months 7-9 versus not receiving treatment during months 7-9 on end-of-study environmental risk outcomes as a function of baseline severity, treatment received between 1-3 months, and severity during months 4-6?" A description of the variables (s0,a1,s1,a2,y) in the data set dat1 used in the illustrative example follows:

s0 binary variable measuring baseline (month 0) severity: 0=low severity, 1=high severity

a1 binary variable measuring receipt of substance use treatment at some point during months 1-3: 0=no, 1=yes

s1 binary variable measuring severity during the 4-6 month interval : 0=low severity, 1=high severity

a2 binary variable measuring receipt of substance use treatment at some point during months 7-9: 0=no, 1=yes

y continuous outcome variable measuring environmental risk for substance use (larger values of y indicates poor outcomes)

The following code excerpt was used to estimate the saturated SNMM presented in the article using the 2-stage estimator (text to the right of the # sign are comments designed to help the reader further connect the R code with Equation 4 in the paper):

modgivpast = list(

list(s0 ~ 1, 0), # used to estimate \delta_{S_0}

list(s1 ~ s0 + a1 + a1:s0, 0) # used to estimate \delta_{S_1}

),

coefnuisfunc = list(

~1, # \eta_1

~s0 + a1 + a1:s0 # \eta_2 + \eta_3 s0 + \eta_4 s0 a1

),

intermeff = list(

list(~a1, ~s0 ), # model for \mu_1: \beta_{10} + \beta_{11} s0

list(~a2, ~s0 + s1 + s0:s1 + # analogous model for \mu_2

a1 + s0:a1 + a1:s1 +

s0:a1:s1)

),

response = ~y, # the outcome

dat = dat1, # the data set

verbose = T,

full.return = T,

nmods.vec = c(1,1) # 1 moderator at time 1; 1 moderator at time 2

)

The code above produced the following output, which includes the results reported in Table 1 of the article:

2-Stage Estimator of the Structural Nested Mean Model

*********************************************************

Outcome...................................... y

Sample size.................................. 2870

Number of Time-points........................ 2

Number of Stage-2 Parameters................. 16

Number of Stage-2 Nuisance Parameters........ 5

Number of Causal Parameters, plus Intercept.. 11

---------------------------------------------------------

Marginal mean response under no treatment

-----------------------------------------

Estimate ASE z value Pr(>|z|)

(SNMM Intercept) 39.715 1.3359 29.7287 0

---------------------------------------------------------

Intermediate causal effect function at time 1

----------------------------------------------

Treatment variable: a1

Estimate ASE z value Pr(>|z|)

(Causal 1 Intercept) -2.6850 1.5831 -1.6960 0.0899

s0 -4.9193 3.0396 -1.6184 0.1056

---------------------------------------------------------

Intermediate causal effect function at time 2

----------------------------------------------

Treatment variable: a2

Estimate ASE z value Pr(>|z|)

(Causal 2 Intercept) -6.6851 4.2766 -1.5632 0.1180

s0 2.7402 5.7166 0.4793 0.6317

s1 22.3874 17.2420 1.2984 0.1941

s0:s1 -28.1620 19.1436 -1.4711 0.1413

a1 1.2861 4.4446 0.2894 0.7723

s0:a1 -0.2535 6.0471 -0.0419 0.9666

s1:a1 -27.0626 17.6041 -1.5373 0.1242

s0:s1:a1 35.7499 19.8079 1.8048 0.0711

---------------------------------------------------------

Non-causal nuisance terms at time 1

------------------------------------

De-meaned moderator: e.s0

Estimate ASE z value Pr(>|z|)

(Nuis (1,1) Intcpt) 7.3954 2.8876 2.5611 0.0104

---------------------------------------------------------

Non-causal nuisance terms at time 2

------------------------------------

De-meaned moderator: e.s1

Estimate ASE z value Pr(>|z|)

(Nuis (2,1) Intcpt) 15.6626 4.1926 3.7357 0.0002

s0 1.5247 6.3846 0.2388 0.8113

a1 -0.4813 4.5343 -0.1061 0.9155

s0:a1 -6.0164 6.7970 -0.8852 0.3761

From the output above, we can see that the estimated effect of a1 on y among the subgroup of individuals who had low severity at baseline (s0=0) is -2.69 (ASE=1.58; p-value=0.09). Similarly, among the subgroup of individuals who had low severity at baseline (s0=0), did not receive treatment between months 1-3 (a1=0), and had low severity between months 4-6 (s1=0), the estimated effect of a2 on y is -6.69 (ASE=4.28; p-value=0.12). The effect of a1 on y for s0=1 versus s0=0 differs by -4.92 (p-value=0.11), suggesting that s0 moderates the effect of a1 on y. Because the data in this example arise from an observational study, we remind the reader that in order to interpret the estimated intermediate effects above as causal, we must satisfy the no direct confounders assumption (see Appendix A of the article);this assumption is likely not met in this illustrative example because there may be additional variables other than s(t-1) that are directly related to both a(t) and y. To make this assumption more plausible, additional baseline and time-varying, continuous and discrete, covariates could be included in the model fitted above but we omitted doing this here for expositional simplicity.

The fitted object snmm1 contains additional information, including the estimates from the stage 1 regression (in snmm1$stage1), the estimated asymptotic variance-covariance matrix of the stage 2 parameters (in snmm1$vcov.theta), and a more succinct summary of the output above (in snmm1$ans).

**Inference concerning linear combinations of the SNMM parameters.** Snmm1$vcov.theta can be used to make inference about linear combinations of the parameters of interest (e.g., What is the effect of a1 on y among the subgroup of individuals who had high severity at baseline (s0=1)?). You will need to know in what order the stage 2 parameters are stored in the variance-covariance matrix (the ordering may be different from the order in which they are presented in the output above) to know how which parameters to pick out for the linear combinations of interest (see L below); you can get the order by looking at snmm1$ans. The following code snippet returns estimates for all possible contrasts (linear combinations) of interest, p-values for the test that the contrasts are zero, 95% confidence intervals, and effect sizes for each of the effects:

c(0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0),

c(0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0),

c(0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0),

c(0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0),

c(0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0),

c(0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0),

c(0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,0),

c(0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0),

c(0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0)

)

Lnames <- c("effect of a1 (no future txt), among s0=0",

"effect of a1 (no future txt), among s0=1",

"effect of a2, among s0=0, a1=0, s1=0",

"effect of a2, among s0=0, a1=1, s1=0",

"effect of a2, among s0=1, a1=0, s1=0",

"effect of a2, among s0=1, a1=1, s1=0",

"effect of a2, among s0=0, a1=0, s1=1",

"effect of a2, among s0=0, a1=1, s1=1",

"effect of a2, among s0=1, a1=0, s1=1",

"effect of a2, among s0=1, a1=1, s1=1"

)

rownames(L) <- Lnames; V <- as.matrix(snmm2.1$vcov.theta)

V.L <- (L%*%V)%*%t(L); se.L <- sqrt( diag(V.L) ) # variance-cov matrix and standard errors

est.L <- L %*% snmm1$ans[,1]

lowCI <- est.L - 1.96*se.L # lower 95% confidence bound

uppCI <- est.L + 1.96*se.L # upper 95% confidence bound

waldval <- est.L/se.L # wald statistic

pval <- 2*pnorm(abs(waldval),lower.tail=F) # p-value for test that linear contrast = 0

sdY <- sd(dat1$y)

effsz <- est.L/sdY # get effect sizes by scaling effects by SD(y)

es.lowCI <- lowCI/sdY

es.uppCI <- uppCI/sdY

res2 <- data.frame(EST=est.L, SE=se.L, CI.LOWER=lowCI, CI.UPPER=uppCI, PVALUE=pval,

EFFSIZE=effsz, ES.CI.LOWER=es.lowCI,ES.CI.UPPER=es.uppCI)

cat("\nLinear Contrasts of Interest\n")

print( round(res2[,1:5],3) ); cat("\n"); print( round(res2[,6:8],3) )

The output is

EST SE CI.LOWER CI.UPPER PVALUE

effect of a1 (no future txt), among s0=0 -2.685 1.583 -5.788 0.418 0.090

effect of a1 (no future txt), among s0=1 -7.604 2.595 -12.690 -2.519 0.003

effect of a2, among s0=0, a1=0, s1=0 -6.685 4.277 -15.067 1.697 0.118

effect of a2, among s0=0, a1=1, s1=0 -5.399 1.210 -7.771 -3.026 0.000

effect of a2, among s0=1, a1=0, s1=0 -3.945 3.793 -11.380 3.490 0.298

effect of a2, among s0=1, a1=1, s1=0 -2.912 1.556 -5.963 0.138 0.061

effect of a2, among s0=0, a1=0, s1=1 15.702 16.703 -17.036 48.441 0.347

effect of a2, among s0=0, a1=1, s1=1 -10.074 3.340 -16.620 -3.528 0.003

effect of a2, among s0=1, a1=0, s1=1 -9.719 7.403 -24.229 4.790 0.189

effect of a2, among s0=1, a1=1, s1=1 0.000 3.292 -6.451 6.452 1.000

EFFSIZE ES.CI.LOWER ES.CI.UPPER

effect of a1 (no future txt), among s0=0 -0.128 -0.276 0.020

effect of a1 (no future txt), among s0=1 -0.362 -0.604 -0.120

effect of a2, among s0=0, a1=0, s1=0 -0.318 -0.717 0.081

effect of a2, among s0=0, a1=1, s1=0 -0.257 -0.370 -0.144

effect of a2, among s0=1, a1=0, s1=0 -0.188 -0.542 0.166

effect of a2, among s0=1, a1=1, s1=0 -0.139 -0.284 0.007

effect of a2, among s0=0, a1=0, s1=1 0.748 -0.811 2.306

effect of a2, among s0=0, a1=1, s1=1 -0.480 -0.791 -0.168

effect of a2, among s0=1, a1=0, s1=1 -0.463 -1.154 0.228

effect of a2, among s0=1, a1=1, s1=1 0.000 -0.307 0.307

The results suggest that treatment reduces environmental risk for substance use except in the following two subgroups in which the impact of treatment during months 7-9 was not statistically significant: in the (s0,a1,s1)=(0,0,1) subgroup, the estimated effect of treatment was harmful and large (p-value=0.35); and in the (s0,a1,s1)=(1,1,1) subgroup, the estimated effect was zero.

**Considering a more parsimonious model.** The saturated SNMM that was fit above (snmm1) (and considered in the article) was chosen for expositional simplicity, i.e., to illustrate the methodology. However, investigators may want to consider estimating a more parsimonious SNMM with fewer parameters. In the model fitted above (snmm1), a variety of parameters were not significant, even at conservative levels (say, using a Type-I error rates of 15%). While formal methods for model selection have not yet been developed and evaluated specifically for the SNMM using the 2-stage regression estimator, we can employ the popular, ad-hoc (but useful) approach of removing terms with small to moderate Wald statistics (say, with |z value| < 1.0). We employ this approach in the following reduced-form model fit:

modgivpast = list(

list(s0 ~ 1, 0),

list(s1 ~ s0 + a1 + a1:s0, 0)

),

coefnuisfunc = list(

~1,

~1 # note reduced-form

),

intermeff = list(

list(~a1, ~s0 ),

list(~a2, ~s1 + s0:s1 + a1:s1 + # note reduced-form

s0:a1:s1)

),

response = ~y,

dat = dat1,

verbose = T,

full.return = T,

nmods.vec = c(1,1)

)

Note that this SNMM is no longer saturated, as we are now making modeling assumptions on the data; that is, there are now fewer parameters than there are subgroup combinations of (s0,a1,s1,a2). The output from this model fit was

2-Stage Estimator of the Structural Nested Mean Model

*********************************************************

Outcome...................................... y

Sample size.................................. 2870

Number of Time-points........................ 2

Number of Stage-2 Parameters................. 10

Number of Stage-2 Nuisance Parameters........ 2

Number of Causal Parameters, plus Intercept.. 8

---------------------------------------------------------

Marginal mean response under no treatment

-----------------------------------------

Estimate ASE z value Pr(>|z|)

(SNMM Intercept) 39.4425 1.2017 32.8219 0

---------------------------------------------------------

Intermediate causal effect function at time 1

----------------------------------------------

Treatment variable: a1

Estimate ASE z value Pr(>|z|)

(Causal 1 Intercept) -2.5743 1.4744 -1.7460 0.0808

s0 -4.4462 2.6738 -1.6629 0.0963

---------------------------------------------------------

Intermediate causal effect function at time 2

----------------------------------------------

Treatment variable: a2

Estimate ASE z value Pr(>|z|)

(Causal 2 Intercept) -4.2998 0.9303 -4.6221 0.0000

s1 22.2760 16.3573 1.3618 0.1732

s1:s0 -25.2739 17.6340 -1.4333 0.1518

s1:a1 -26.3297 16.5613 -1.5898 0.1119

s1:s0:a1 31.4240 18.1437 1.7320 0.0833

---------------------------------------------------------

Non-causal nuisance terms at time 1

------------------------------------

De-meaned moderator: e.s0

Estimate ASE z value Pr(>|z|)

(Nuis (1,1) Intcpt) 7.6182 2.5355 3.0046 0.0027

---------------------------------------------------------

Non-causal nuisance terms at time 2

------------------------------------

De-meaned moderator: e.s1

Estimate ASE z value Pr(>|z|)

(Nuis (2,1) Intcpt) 13.3313 1.0962 12.1618 0

And the subgroup specific effects (linear combinations of interest) are now

EST SE CI.LOWER CI.UPPER PVALUE

effect of a1 (no future txt), among s0=0 -2.574 1.474 -5.464 0.316 0.081

effect of a1 (no future txt), among s0=1 -7.021 2.233 -11.397 -2.644 0.002

effect of a2, among s0=0, a1=0, s1=0 -4.300 0.930 -6.123 -2.476 0.000

effect of a2, among s0=0, a1=1, s1=0 -4.300 0.930 -6.123 -2.476 0.000

effect of a2, among s0=1, a1=0, s1=0 -4.300 0.930 -6.123 -2.476 0.000

effect of a2, among s0=1, a1=1, s1=0 -4.300 0.930 -6.123 -2.476 0.000

effect of a2, among s0=0, a1=0, s1=1 17.976 16.322 -14.016 49.968 0.271

effect of a2, among s0=0, a1=1, s1=1 -8.353 3.118 -14.466 -2.241 0.007

effect of a2, among s0=1, a1=0, s1=1 -7.298 6.777 -20.580 5.985 0.282

effect of a2, among s0=1, a1=1, s1=1 -2.203 3.193 -8.462 4.055 0.490

EFFSIZE ES.CI.LOWER ES.CI.UPPER

effect of a1 (no future txt), among s0=0 -0.123 -0.260 0.015

effect of a1 (no future txt), among s0=1 -0.334 -0.543 -0.126

effect of a2, among s0=0, a1=0, s1=0 -0.205 -0.292 -0.118

effect of a2, among s0=0, a1=1, s1=0 -0.205 -0.292 -0.118

effect of a2, among s0=1, a1=0, s1=0 -0.205 -0.292 -0.118

effect of a2, among s0=1, a1=1, s1=0 -0.205 -0.292 -0.118

effect of a2, among s0=0, a1=0, s1=1 0.856 -0.667 2.379

effect of a2, among s0=0, a1=1, s1=1 -0.398 -0.689 -0.107

effect of a2, among s0=1, a1=0, s1=1 -0.347 -0.980 0.285

effect of a2, among s0=1, a1=1, s1=1 -0.105 -0.403 0.193

The subgroup effects of a1 do not change much in the reduced model: treatment during months 1-3 reduces environmental risk for substance use, and the effect is strongest among subjects exhibiting higher severity at baseline. Concerning the effects of a2, because we have removed terms from the model for the intermediate causal effect at time t=2, we now see that the impact of a2 on y is the same among certain subgroups of subjects. In particular, the results suggest that s1 is the primary moderator of the impact of a2 on y. There is a beneficial effect of treatment during months 7-9, and it is stronger in general among subjects with high intermediate severity (s1=1) (it almost doubles). As with the saturated SNMM, there are two exceptions: subjects who had high severity at baseline, were treated during months 1-3, and continued to have high severity during months 4-6 (for this subgroup, the effect is beneficial, but small and non-significant); and subjects who at intake had low severity, did not receive treatment during months 1-3, but whose severity was high during months 4-6 (for this group, we see a large, harmful effect of treatment, but it is not significant).

### Disclaimer

This R code is provided with no guarantees. It is an alpha version, written originally for methodological research purposes only. Therefore, it may not be user-friendly and may be difficult to use. If you find bugs/errors in this code, please email me (dalmiral@umich.edu): I will fix the code and then acknowledge you on this website.

### Recommended Citation

If you use this code or any of this work in your own research, please cite any one (or both) of the related articles:

Almirall, D., Ten Have, T., & Murphy, S. A. (2010). Structural nested mean models for assessing time-varying effect moderation.* Biometrics, 66*(1), 131-139. PMCID: PMC2875310

Almirall, D., McCaffrey, D. F., Ramchand, R., & Murphy, S. A. (2011). Subgroups analysis when treatment and moderators are time-varying.* Prevention Science. *Advance online publication. * *doi:10.1007/s11121-011-0208-7 PMCID: PMC3135740

### Additional References

The original developer of the structural nested mean model is Professor James M. Robins at Harvard University. You will find additional articles on the SNMM authored by Dr. Robins by visiting his website.

### Acknowledgements

Funding for this research was provided by the following grants: R01-DA-015697 (McCaffrey), R01-DA-017507 (Ramchand), R01-MH-080015 (Murphy), and P50-DA-010075 (Murphy). The authors would like to thank Andrew R. Morral, Beth Ann Griffin, and Scott N. Compton for comments and suggestions, and Cha-Chi Fan for guidance with the data.