Hi everyone,
So I have this code put together below and it runs perfectly fine, but I noticed an unexpected finding. Now, of course, sometimes unexpected findings are warranted and constitute a new insight in research, but I am inclined to think that this isn't one of those times. This simulation looks at how sample size relates to statistical power generated from models analyzed with generalized estimating equations. The expectation is that power should increase positively with greater sample size, and while that holds true for the first few sample sizes (3 and 6, I notice that the power goes to zero at 9 and this is completely bizarre). I don't see anything wrong with the code, but wondering what else might give rise to this break in trend. Or is it a real finding?
For clarity, here is where I am manipulating the sample size in the code below:
So I have this code put together below and it runs perfectly fine, but I noticed an unexpected finding. Now, of course, sometimes unexpected findings are warranted and constitute a new insight in research, but I am inclined to think that this isn't one of those times. This simulation looks at how sample size relates to statistical power generated from models analyzed with generalized estimating equations. The expectation is that power should increase positively with greater sample size, and while that holds true for the first few sample sizes (3 and 6, I notice that the power goes to zero at 9 and this is completely bizarre). I don't see anything wrong with the code, but wondering what else might give rise to this break in trend. Or is it a real finding?
For clarity, here is where I am manipulating the sample size in the code below:
Code:
local num_clus 3 6 9
Code:
capture program drop swcrt program define swcrt, rclass args num_clus clussize intercept intrvcoeff timecoeff1 timecoeff2 timecoeff3 timecoeff4 timecoeff5 sigma_u3 sigma_error alpha assert `num_clus' > 0 & `clussize' > 0 & `intercept' > 0 & `intrvcoeff' > 0 & `timecoeff1' < 0 & `timecoeff2' < 0 & `timecoeff3' < 0 & `timecoeff4' < 0 & `timecoeff5' < 0 & `sigma_u3' > 0 & `sigma_error' > 0 & `alpha' > 0 *Generate simulated multi—level data qui clear set obs `num_clus' qui gen cluster = _n qui gen group = 1+mod(_n-1,4) *Generate cluster-level errors qui gen u_3 = rnormal(0,`sigma_u3') expand `clussize' bysort cluster: gen individual = _n *Set up time expand 6 bysort cluster individual: gen time = _n-1 *Set up intervention variable gen intrv = (time>=group) *Generate residual errors qui gen error = rnormal(0,`sigma_error') *Generate outcome y qui gen y = `intercept' + `intrvcoeff'*intrv + `timecoeff1'*1.time + `timecoeff2'*2.time + `timecoeff3'*3.time + `timecoeff4'*4.time + `timecoeff5'*5.time + u_3 + error *Fit multi-level model to simulated dataset xtgeebcv y intrv i.time, cluster(cluster) family(gaussian) link(identity) corr(exchangeable) stderr(kc) *Return estimated effect size, absolute and relative bias, p-value, and significance dichotomy tempname M matrix `M' = r(table) return scalar b = _b[intrv] return scalar abias = _b[intrv] - `intrvcoeff' return scalar rbias = [(_b[intrv] - `intrvcoeff')/`intrvcoeff']*100 return scalar p = `M'[4,1] return scalar p_= (`M'[4,1] < `alpha') exit end swcrt *Stepped-wedge specifications local num_clus 3 6 9 local clussize 5 10 15 20 25 *Model specifications local intercept 17.87 local timecoeff1 -5.42 local timecoeff2 -5.72 local timecoeff3 -7.03 local timecoeff4 -6.13 local timecoeff5 -9.13 local intrvcoeff 5.00 local sigma_u3 25.77 local sigma_u2 120.62 local sigma_error 38.35 local alpha 0.05 local nrep 500 *Postfile to store results tempname step tempfile powerresults capture postutil clear postfile `step' num_clus clussize intrvcoeff b p p_ abias rbias using `powerresults.dta', replace *Loop over number of clusters foreach c of local num_clus{ display as text "Number of clusters" as result "`c'" foreach s of local clussize { display as text "Cluster size" as result "`s'" forvalue i = 1/`nrep'{ display as text "Iterations" as result `nrep' quietly swcrt `c' `s' `intercept' `intrvcoeff' `timecoeff1' `timecoeff2' `timecoeff3' `timecoeff4' `timecoeff5' `sigma_u3' `sigma_error' `alpha' post `step' (`c') (`s') (`intrvcoeff') (`r(b)') (`r(p)') (`r(p_)') (`r(abias)') (`r(rbias)') } } } postclose `step' *Open results, calculate power use `powerresults', clear levelsof num_clus, local(num_clus) levelsof clussize, local(clussize) matrix drop _all *Loop over combinations of clusters *Add power results to matrix foreach c of local num_clus{ foreach s of local clussize{ quietly ci proportions p_ if num_clus == `c' & clussize == `s' local power `r(proportion)' quietly ci mean abias if num_clus == `c' & clussize == `s' local abias `r(mean)' quietly ci mean rbias if num_clus == `c' & clussize == `s' local rbias `r(mean)' matrix M = nullmat(M) \ (`c', `s', `intrvcoeff', `power', `abias',`rbias') } } *Display the matrix matrix colnames M = c s intrvcoeff power abias rbias matrix list M, noheader format(%3.2f)
Comment