I am a big fan of scientific collaboration. Feel free to contact me to discuss your causal inference related projects for potential collaboration.
- This package simulates survival data suitable for fitting Marginal Structural Model.
library(devtools)
install_github("ehsanx/simMSM")require(simMSM)?simmsmsetwd("C:/data") # change working dirsimmsm(subjects=2500, tpoints=10, psi=0.3, n=1000)
# This code generates 1000 datasets (takes time!)# 2500 subjects in each datasets# Each subject followed upto 10 time-points (say, months)# Causal effect (log-odds) is 0.3| Parameter | Description |
|---|---|
| subjects | Number of Subjects in each simulated dataset |
| tpoints | Maximum number of time-points each subjects are followed |
| psi | Causal effect parameter for Marginal Structural Model |
| n | Number of simulated datasets an user wants to generate |
Below is an example where we generate 1 large data with n=10,000 subjects, each with up to 10 followup times. Log-odds is 0.3.
require(simMSM)
simmsm(subjects=10000, tpoints=10, psi=0.3, n=1)
aggregate.data<- read.csv("r1.csv")
head(aggregate.data)
# id tpoint tpoint2 T0 IT0 chk Y Ym A Am1 L Lm1 Am1L pAt T maxT pL psi#1 1 1 0 75.51818 0 1.000000 0 0 0 0 0 0 0 0.2222222 NA 10 0.3000000 0.3#2 1 2 1 75.51818 0 2.000000 0 0 0 0 1 0 0 0.3202196 NA 10 0.3000000 0.3#3 1 3 2 75.51818 0 3.349859 0 0 1 0 1 1 0 0.4371436 NA 10 0.3913043 0.3#4 1 4 3 75.51818 0 4.699718 0 0 1 1 0 1 0 0.6532898 NA 10 0.2432432 0.3#5 1 5 4 75.51818 0 6.049576 0 0 1 1 0 0 0 0.5333333 NA 10 0.1764706 0.3#6 1 6 5 75.51818 0 7.049576 0 0 0 1 0 0 0 0.5333333 NA 10 0.1764706 0.3# select.seed.valx#1 1#2 1#3 1#4 1#5 1#6 1Here is a helper function to extract results
ext.cox<-function(fit){
x<- summary(fit)
b.se<- ifelse(x$used.robust==TRUE, x$coef["A","robust se"], x$coef["A","se(coef)"])
if (is.null(fit$weights)) {
res<- as.numeric(c(x$n, x$nevent, x$coef["A","coef"], b.se, confint(fit)["A",],x$coef["A", "Pr(>|z|)"], x$used.robust,
rep(NA,3) ))
} else {
res<- as.numeric(c(x$n, x$nevent, x$coef["A","coef"], b.se, confint(fit)["A",],x$coef["A", "Pr(>|z|)"], x$used.robust,
mean(fit$weights), range(fit$weights) ))
}
names(res) <- c("n", "events", "coef", "se", "lowerci", "upperci", "pval", "robust", "meanW", "minW", "maxW")
return(res)
}Below are codes of estimating weights
ww<- glm(A~tpoint+L+Am1+Lm1, family= binomial(logit), data=aggregate.data)
# Step 2: Weight numerator modelww0<- glm(A~tpoint+Am1, family= binomial(logit), data=aggregate.data)
# Step 3: Obtain fir from the weight modelsaggregate.data$wwp<- with(aggregate.data, ifelse(A==0, 1- fitted(ww), fitted(ww)))
aggregate.data$wwp0<- with(aggregate.data, ifelse(A==0, 1- fitted(ww0),fitted(ww0)))
# Step 4: time-dependent weightsaggregate.data$w<- unlist(tapply(1/aggregate.data$wwp, aggregate.data$id, cumprod))
aggregate.data$sw<- unlist(tapply(aggregate.data$wwp0/aggregate.data$wwp, aggregate.data$id, cumprod))
summary(aggregate.data$sw)
# Min. 1st Qu. Median Mean 3rd Qu. Max. # 0.1396 0.7329 0.9091 1.0013 1.1788 6.2040 Below are codes of estimating treatment effect from the weighted outcome model
require(survival)
# Step 5: Weighted outcome modelfit.msm.w<- coxph(Surv(tpoint2, tpoint, Y) ~A+ cluster(id), data=aggregate.data, weight=w, robust=TRUE)
ext.cox(fit.msm.w)
fit.msm<- coxph(Surv(tpoint2, tpoint, Y) ~A+ cluster(id), data=aggregate.data, weight=sw, robust=TRUE)
ext.cox(fit.msm)
# n events coef se lowerci upperci pval robust #9.481300e+04 1.147000e+03 3.341396e-01 6.882629e-02 1.992425e-01 4.690367e-01 1.204931e-06 1.000000e+00 # meanW minW maxW #1.001264e+00 1.395881e-01 6.203989e+00 Note that log-odds (psi) was 0.3 when we generated the data, and you are getting back an estimate of 0.3341396.
Below are additional codes for unweighted models
# Comparison with unweighted modelsfit.cox<- coxph(Surv(tpoint2, tpoint, Y) ~A+ cluster(id), data=aggregate.data, robust=TRUE)
ext.cox(fit.cox)
fit.cox.adj<- coxph(Surv(tpoint2, tpoint, Y) ~A+L+ cluster(id), data=aggregate.data, robust=TRUE)
ext.cox(fit.cox.adj)- Ehsan Karim
(only R porting from the SAS code). I wrote them in R basically to understand the mechanism, but the SAS / SAS IML / Stata codes (I have them as well, available upon request) are faster than this. Feel free to report any errors / update suggestions.
- Young J.G., Hernan M.A., Picciotto S., and Robins J.M. Relation between three classes of structural models for the effect of a time-varying exposure on survival. Lifetime Data Analysis, 16(1):71-84, 2010.
- Young, Jessica G., et al. Simulation from structural survival models under complex time-varying data structures. JSM proceedings, section on statistics in epidemiology. American Statistical Association, Denver, CO (2008)
- Ali R.A., Ali M.A., and Wei Z. On computing standard errors for marginal structural cox models. Lifetime data analysis, pages 1–26, 2013. doi: 10.1007/s10985-013-9255-7.
- Xiao Y., Abrahamowicz M., and Moodie E.E.M. Accuracy of conventional and marginal structural Cox model estimators: A simulation study. The International Journal of Biostatistics, 6(2):1–28, 2010.
- Karim, M. E.; Petkau, J.; Gustafson, P.; Platt, R.; Tremlett, H. and BeAMS study group. Comparison of Statistical Approaches Dealing with Time-dependent Confounding in Drug Effectiveness Studies. Statistical Methods in Medical Research. First published online: September 21. doi: 10.1177/0962280216668554
- Karim, M. E.; Petkau, J.; Gustafson, P.; Tremlett, H. and BeAMS study group. On the application of statistical learning approaches to construct inverse probability weights in marginal structural Cox models: hedging against weight-model misspecification. Communications in Statistics - Simulation and Computation (In Press).
- You can get additional descriptions of MSM analysis in this GitHub repo.
- Stata users can take a look at this GitHub repo