R Implementation of the 2-Stage Regression Estimator of the Structural Nested Mean Model for Examining Time-Varying Causal Effect Moderation

(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

         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):

 

snmm1 <- twostage.snmm(
  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 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:

 

L <- rbind( c(0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
            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

 

Linear Contrasts of Interest
                                             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:

 

snmm2 <- twostage.snmm(
           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

 

Linear Contrasts of Interest
                                            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 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.

Follow Us