Robust inference for the stepped wedge design

Abstract

Stepped wedge designed trials are a type of cluster‐randomized study in which the intervention is introduced to each cluster in a random order over time. This design is often used to assess the effect of a new intervention as it is rolled out across a series of clinics or communities. Based on a permutation argument, we derive a closed‐form expression for an estimate of the intervention effect, along with its standard error, for a stepped wedge design trial. We show that these estimates are robust to misspecification of both the mean and covariance structure of the underlying data‐generating mechanism, thereby providing a robust approach to inference for the intervention effect in stepped wedge designs. We use simulations to evaluate the type 1 error and power of the proposed estimate and to compare the performance of the proposed estimate to the optimal estimate when the correct model specification is known. The limitations, possible extensions, and open problems regarding the method are discussed.

A Simulation Study of Statistical Approaches to Data Analysis in the Stepped Wedge Design

“…# Null condition, 40 clusters, with random treatment effect, time effect, correlation = 0 gm1a = NULL gm2a = NULL gmz1 = NULL gmz2 = NULL gm3a = NULL gm4a = NULL gmz3 = NULL gmz4 = NULL for(i in 1:500) { design <-swDsn(clusters=c(10,10,10,10), extra.time=0, all.ctl.time0=TRUE) # mu1 is the difference between a0 and a1 # tau: random intercept # eta: random treatment swGenData.nScalar <-swSim(design,family=gaussian(link="identity"), n=100, mu0=3, mu1=3, time.effect=c(0,.2,.3,.4,.5), sigma=1, tau=2, eta=0, rho=0, seed=i, retTimeOnTx=FALSE) # time effect, no random treatment effect glmod1 <-glmer(response.var ~ tx.var+as.factor(time.var)+(1|cluster.var), data = swGenData.nScalar, gaussian, nAGQ = 10) storep1, storep1_z, storep2, storep2_z, storep3, storep3_z) zscore1 = summary(mod1)$coefficients [2,5] mod2 <-gee(response.var~tx.var+as.factor(time.var), data = swdata, family = "gaussian",id = cluster.var, corstr = "exchangeable") result2 = mod2$coefficients myfunc <-function(n, seed) { design <-swDsn(clusters=c(10,10,10,10), extra.time=0, all.ctl.time0=TRUE) # first generate the original dataset with a time effect and 0 treatment effect, # rho is the correlation between slope and intercept swGenData.nScalar <-swSim(design,family=gaussian(link="identity"), n=10, mu0=3, myfunc <-function(n, seed) { design <-swDsn(clusters=c(10,10,10,10), extra.time=0, all.ctl.time0=TRUE) # first generate the original dataset with a time effect and 0 treatment effect, # rho is the correlation between slope and intercept swGenData.nScalar <-swSim(design,family=gaussian(link="identity"), n=10, mu0=3, zscore1 = summary(mod1)$coefficients [2,5] mod2 <-gee(response.var~tx.var+as.factor(time.var), data = swdata, family = "gaussian",id = cluster.var, corstr = "exchangeable") result2 = mod2$coefficients[ [2]][1] zscore2 = summary(mod2)$coefficients [2,5] mod3 <-gee(response.var~tx.var, data = swdata, family = "gaussian",id = cluster.var, corstr = "independence") result3 = mod3$coefficients [ storep1, storep1_z, storep2, storep2_z, storep3, storep3_z)) } finalp1 = 0 finalp1_z = 0 finalp2 = 0 finalp2_z = 0 finalp3 = 0 finalp3_z = 0 for(i in 1:500) { resultList = NULL resultList = checkfunc(1000, 40, i) finalp1 = finalp1 + resultList [1] finalp1_z = finalp1_z + resultList [2] finalp2 = finalp2 + resultList [3] finalp2_z = finalp2_z + resultList [4] finalp3 = finalp3 + resultList [5] finalp3_z = finalp3_z + resultList[6] } print(finalp1/500) print(finalp1_z/500) print(finalp2/500) print(finalp2_z/500) print(finalp3/500) print(finalp3_z/500)…”

“…myfunc <-function(n, seed) { design <-swDsn(clusters=c(10,10,10,10), extra.time=0, all.ctl.time0=TRUE) # first generate the original dataset with a time effect and 0 treatment effect, # rho is the correlation between slope and intercept swGenData.nScalar <-swSim(design,family=gaussian(link="identity"), n=10, mu0=3, mu1=3.21, time.effect=c(0,.2,.3,.4,.5), sigma=1, tau=2, eta=0, rho=0, seed = seed, retTimeOnTx=FALSE) data = swSummary(response.var, tx.var, time.var, cluster.var, swGenData.nScalar, type="mean", digits=3)$swDsn newdata = NULL listn=c(1:n) x = sample(listn) newdata=data[x,] return(newdata) } checkfunc = function(numberOfPermutation, nCluster, seed) { design <-swDsn(clusters=c(10,10,10,10), extra.time=0, all.ctl.time0=TRUE) # first generate the original dataset with a time effect and 0 treatment effect, # rho is the correlation between slope and intercept swdata <-swSim(design,family=gaussian(link="identity"), n=10, mu0=3, mu1=3.21, time.effect=c(0,.2,.3,.4,.5), sigma=1, tau=2, eta=0, rho=0, seed = seed, retTimeOnTx=FALSE) storep1 = 0 storep1_z = 0 storep2 = 0 storep2_z = 0 storep3 = 0 storep3_z = 0 matrix = myfunc(nCluster, seed) mod1 <-gee(response.var~tx.var+as.factor(time.var), data = swdata, family = "gaussian",id = cluster.var, corstr = "independence") result1 = mod1$coefficients[ [2]][1] zscore1 = summary(mod1)$coefficients [2,5] mod2 <-gee(response.var~tx.var+as.factor(time.var), data = swdata, family = "gaussian",id = cluster.var, corstr = "exchangeable") result2 = mod2$coefficients storep1, storep1_z, storep2, storep2_z, storep3, storep3_z)) } finalp1 = 0 finalp1_z = 0 finalp2 = 0 finalp2_z = 0 finalp3 = 0 finalp3_z = 0 for(i in 1:500) { resultList = NULL resultList = checkfunc(1000, 40, i) finalp1 = finalp1 + resultList [1] finalp1_z = finalp1_z + resultList [2] finalp2 = finalp2 + resultList [3] finalp2_z = finalp2_z + resultList [4] finalp3 = finalp3 + resultList [5] finalp3_z = finalp3_z + resultList[6] } print(finalp1/500) print(finalp1_z/500) print(finalp2/500) print(finalp2_z/500) print(finalp3/500) print(finalp3_z/500)…”

“…myfunc <-function(n, seed) { design <-swDsn(clusters=c (5,5,5,5), extra.time=0, all.ctl.time0=TRUE) # first generate the original dataset with a time effect and 0 treatment effect, # rho is the correlation between slope and intercept swGenData.nScalar <-swSim(design,family=gaussian(link="identity"), n=10, mu0=3, c(storep1, storep1_z, storep2, storep2_z, storep3, storep3_z, storep4, storep4_z)) } finalp1 = 0 finalp1_z = 0 finalp2 = 0 finalp2_z = 0 finalp3 = 0 finalp3_z = 0 finalp4 = 0 finalp4_z = 0 for(i in 1:500) { resultList = NULL resultList = checkfunc(1000, 20, i) finalp1 = finalp1 + resultList [1] finalp1_z = finalp1_z + resultList [2] finalp2 = finalp2 + resultList [3] finalp2_z = finalp2_z + resultList [4] finalp3 = finalp3 + resultList [5] finalp3_z = finalp3_z + resultList [6] finalp4 = finalp4 + resultList [7] finalp4_z = finalp4_z + resultList[8] } print(finalp1/500) print(finalp1_z/500) print(finalp2/500) print(finalp2_z/500) print(finalp3/500) print(finalp3_z/500) print(finalp4/500) print(finalp4_z/500)…”

“…myfunc <-function(n, seed) { design <-swDsn(clusters=c(10,10,10,10), extra.time=0, all.ctl.time0=TRUE) # first generate the original dataset with a time effect and 0 treatment effect, # rho is the correlation between slope and intercept swGenData.nScalar <-swSim(design,family=gaussian(link="identity"), n=10, mu0=3, mu1=3. gmz1 = summary(glmod1)$coefficients [2,3] # with time effect, no random treatment effect glmod2 <-lmer(response.var ~ tx.var+as.factor(time.var)+(1|cluster.var), data = data) gm2a = fixef(glmod2)[ [2]] [1] gmz2 = summary(glmod2)$coefficients [2,3] # with random treatment and time effect, intercept and slope are not correlated glmod3 <-lmer(response.var ~ tx.var+as.factor(time.var) + (1|cluster.var) + (tx.var-1|cluster.var), data = data) gm3a = fixef(glmod3)[ [2]] [1] gmz3 = summary(glmod3)$coefficients [2,3] # with random treatment and time effect, intercept and slope are correlated glmod4 <-lmer(response.var ~ tx.var+as.factor(time.var) + (tx.var+1|cluster.var), data = data) gm4a = fixef(glmod4)[ [2]] [1] gmz4 = summary(glmod4)$coefficients [2,3] dataPerm = data # copy the original dataset gm1a_1 = NULL gmz1_1 = NULL gm2a_1 = NULL gmz2_1 = NULL gm3a_1 = NULL gmz3_1 = NULL gm4a_1 = NULL gmz4_1 = NULL return(c(storep1, storep1_z, storep2, storep2_z, storep3, storep3_z, storep4, storep4_z)) } finalp1 = 0 finalp1_z = 0 finalp2 = 0 finalp2_z = 0 finalp3 = 0 finalp3_z = 0 finalp4 = 0 finalp4_z = 0 for(i in 1:500) { resultList = NULL resultList = checkfunc(1000, 40, i) finalp1 = finalp1 + resultList [1] finalp1_z = finalp1_z + resultList [2] finalp2 = finalp2 + resultList [3] finalp2_z = finalp2_z + resultList [4] finalp3 = finalp3 + resultList [5] finalp3_z = finalp3_z + resultList [6] finalp4 = finalp4 + resultList [7] finalp4_z = finalp4_z + resultList[8] } print(finalp1/500) print(finalp1_z/500) print(finalp2/500) print(finalp2_z/500) print(finalp3/500) print(finalp3_z/500) print(finalp4/500) print(finalp4_z/500)…”

“…resultList = checkfunc(1000, 20, i) finalp1 = finalp1 + resultList[1] finalp1_z = finalp1_z + resultList[2] finalp2 = finalp2 + resultList[3] finalp2_z = finalp2_z + resultList[4] finalp3 = finalp3 + resultList[5] finalp3_z = finalp3_z + resultList[6] } print(finalp1/500) print(finalp1_z/500) print(finalp2/500) print(finalp2_z/500) print(finalp3/500) print(finalp3_z/500)# S2myfunc <-function(n, seed) { design <-swDsn(clusters=c(5,5,5,5), extra.time=0, all.ctl.time0=TRUE) # first generate the original dataset with a time effect and 0 treatment effect, # rho is the correlation between slope and intercept swGenData.nScalar <-swSim(design,family=gaussian(link="identity"), n=10, mu0=3, mu1=3, time.effect=c(0,.2,.3,.4,.5), sigma=1, tau=2, eta=2, rho=0, seed = seed, retTimeOnTx=FALSE) data = swSummary(response.var, tx.var, time.var, cluster.var, swGenData.nScalar, type="mean", digits=3)$swDsn newdata = NULL listn=c(1:n) x = sample(listn) newdata=data[x,] return(newdata) } checkfunc = function(numberOfPermutation, nCluster, seed) { design <-swDsn(clusters=c(5,5,5,5), extra.time=0, all.ctl.time0=TRUE) # first generate the original dataset with a time effect and 0 treatment effect, # rho is the correlation between slope and intercept swdata <-swSim(design,family=gaussian(link="identity"), n=10, mu0=3, mu1=3, time.effect=c(0,.2,.3,.4,.5), sigma=1, tau=2, eta=2, rho=0, seed = seed, retTimeOnTx=FALSE) (nCluster, seed) mod1 <-gee(response.var~tx.var+as.factor(time.var), data = swdata, family = "gaussian",id = cluster.var, corstr = "independence") result1 = mod1$coefficients[[2]][1] zscore1 = summary(mod1)$coefficients[2,5] mod2 <-gee(response.var~tx.var+as.factor(time.var), data = swdata, family = "gaussian",id = cluster.var, corstr = "exchangeable") result2 = mod2$coefficients[[2]][1] zscore2 = summary(mod2)$coefficients[2,5] mod3 <-gee(response.var~tx.var, data = swdata, family = "gaussian",id = cluster.var, corstr = "independence") result3 = mod3$coefficients[[2]][1] zscore3 = summary(mod3)$coefficients[2,5] dataPerm = swdata # copy the original dataset result1_1 nCluster) x = sample(listn) newdata = matrix[x,] for(a in 1:nrow(newdata)) { for(b in 1:ncol(newdata)){ if(newdata[a,b]==0) { newtx = append(newtx, rep(0,10)) } else{ newtx = append(newtx, rep(1,10)) } } } dataPerm$tx.var = as.vector(t(newtx))mod1_1 <-gee(response.var~tx.var+as.factor(time.var), data = dataPerm, family = "gaussian",id = cluster.var, corstr = "independence")result1_1[i] = mod1_1$coefficients[[2]][1] zscore1_1[i] = summary(mod1_1)$coefficients[2,5]mod2_1 <-gee(response.var~tx.var+as.factor(time.var), data = dataPerm, family = "gaussian",id = cluster.var, corstr = "exchangeable")result2_1[i] = mod2_1$coefficients[[2]][1] zscore2_1[i] = summary(mod2_1)$coefficients[2,5]mod3_1 <-gee(response.var~tx.var, data = dataPerm, family = "gaussian",id = cluster.var, corstr = "independence")result3_1[i] = mod3_1$coefficients[[2]][1] zscore3_1[i] = summary(mod3_1)$coefficients[2,5] } sort1 = sort(result1_1) sort2 = sort(result2_1) sort3 = sort(result3_1) sortz1 = sort(zscore1_1) sortz2 = sort(zscore2_1) sortz3 = sort(zscore3_1) # z score if((zscore1< sortz1[25]) | (zscore1> sortz1[975])) { storep1_z = 1} else {storep1_z = 0} if((zscore2< sortz2[25]) | (zscore2> sortz2[975])) { storep2_z = 1} ...…”