Hello,
I have written the following code for a three-level repeated measures model in a stepped-wedged randomized trial design context. The empty local macros have yet to be filled in with variables.
*Generate multi-level data
capture program drop swcrt
program define swcrt, rclass
version 15.1
preserve
clear
*Stepped-wedge specifications
local num_clus 3 6 9 18 36
local ICC_lvl3
local ICC_lvl2
local clussize 25
*Model specifications
local intercept
local timecoeff
local intcoeff
local sigma_u3
local sigma_u2
local sigma_error
local nrep 100
local alpha 0.05
args num_clus ICC_lvl3 ICC_lvl2 clussize intercept timecoeff intrvcoeff
assert `num_clus' > 0 `ICC_lvl3' > 0 `ICC_lvl2' > 0 `clussize' > 0 `intercept' > 0 `timecoeff' > 0 `intrvcoeff' > 0 `u_3' > 0 `u_2' > 0 `error' > 0
/*Generate simulated multi—level data*/
qui
set seed 12345
clear
set obs `num_clus'
qui gen cluster = _n
//Generate cluster-level errors //
qui gen u_3 = rnormal(0,`sigma_u3')
expand `clussize'
bysort cluster: gen individual = _n
//Generate patient-level errors //
qui gen u_2 = rnormal(0,`sigma_u2')
expand 6
//Set up time//
bysort cluster individual: gen time = _n-1
//Set up intervention variable//
gen intrv = (time>=cluster)
//Generate residual errors
qui gen error = rnormal(0,`sigma_error')
//Generate outcome y
qui gen y = `intercept' + `timecoeff'*time + `intrvcoeff'*intrv + u_3 + u_2 + error
/*Fit multi-level model to simulated dataset*/
mixed y intrv i.time ||cluster: ||individual:, covariance(unstructured) reml dfmethod(kroger)
return scalar p = r(p)
drop _all
end swcrt
*Postfile to store results
tempfile powerresults
capture postutil clear
postfile step int num_clus intrv ICC using `results'
*Loop over number of clusters, ICC
foreach x of local `num_clus'{
display as text "Number of clusters" as result "`c'"
foreach x of local `ICC'{
display as text "Intra-class correlation" as result `i'
forvalue i = 1/`nrep'{
display as text "Iterations" as result `nrep'
quietly swcrt `c' `ICC' `intrvcoeff' `u_3' `u_2' `error'
post step (`c') (`ICC') (`r(p)')
}
}
}
postclose step
use `powerresults', clear
Have I coded this correctly/are there any improvements that can be made? The code as is (I think?) allows me to calculate the minimum number of clusters that would satisfy a power of 80%, but how would I also evaluate a) Type 1 error rate, and b) Bias of the intervention effect estimate? I have looked everywhere, but don't see any sample code that evaluates these measures.
Would welcome everyone's thoughts.
I have written the following code for a three-level repeated measures model in a stepped-wedged randomized trial design context. The empty local macros have yet to be filled in with variables.
*Generate multi-level data
capture program drop swcrt
program define swcrt, rclass
version 15.1
preserve
clear
*Stepped-wedge specifications
local num_clus 3 6 9 18 36
local ICC_lvl3
local ICC_lvl2
local clussize 25
*Model specifications
local intercept
local timecoeff
local intcoeff
local sigma_u3
local sigma_u2
local sigma_error
local nrep 100
local alpha 0.05
args num_clus ICC_lvl3 ICC_lvl2 clussize intercept timecoeff intrvcoeff
assert `num_clus' > 0 `ICC_lvl3' > 0 `ICC_lvl2' > 0 `clussize' > 0 `intercept' > 0 `timecoeff' > 0 `intrvcoeff' > 0 `u_3' > 0 `u_2' > 0 `error' > 0
/*Generate simulated multi—level data*/
qui
set seed 12345
clear
set obs `num_clus'
qui gen cluster = _n
//Generate cluster-level errors //
qui gen u_3 = rnormal(0,`sigma_u3')
expand `clussize'
bysort cluster: gen individual = _n
//Generate patient-level errors //
qui gen u_2 = rnormal(0,`sigma_u2')
expand 6
//Set up time//
bysort cluster individual: gen time = _n-1
//Set up intervention variable//
gen intrv = (time>=cluster)
//Generate residual errors
qui gen error = rnormal(0,`sigma_error')
//Generate outcome y
qui gen y = `intercept' + `timecoeff'*time + `intrvcoeff'*intrv + u_3 + u_2 + error
/*Fit multi-level model to simulated dataset*/
mixed y intrv i.time ||cluster: ||individual:, covariance(unstructured) reml dfmethod(kroger)
return scalar p = r(p)
drop _all
end swcrt
*Postfile to store results
tempfile powerresults
capture postutil clear
postfile step int num_clus intrv ICC using `results'
*Loop over number of clusters, ICC
foreach x of local `num_clus'{
display as text "Number of clusters" as result "`c'"
foreach x of local `ICC'{
display as text "Intra-class correlation" as result `i'
forvalue i = 1/`nrep'{
display as text "Iterations" as result `nrep'
quietly swcrt `c' `ICC' `intrvcoeff' `u_3' `u_2' `error'
post step (`c') (`ICC') (`r(p)')
}
}
}
postclose step
use `powerresults', clear
Have I coded this correctly/are there any improvements that can be made? The code as is (I think?) allows me to calculate the minimum number of clusters that would satisfy a power of 80%, but how would I also evaluate a) Type 1 error rate, and b) Bias of the intervention effect estimate? I have looked everywhere, but don't see any sample code that evaluates these measures.
Would welcome everyone's thoughts.
Comment