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