*! version 1.1.0 25feb2022 capture program drop srtfreq program define srtfreq, rclass version 15.1 syntax varlist(numeric fv) [if] [in], INTervention(varlist numeric fv max=1) [, NPerm(integer 0) NBoot(integer 0) SEED(integer 1020252) noDOT PERCentile BASIC noIsily PASTE *] quietly { preserve if "`paste'"!="" { if "`nperm'" != "0" { cap drop PermC_I* PermUnc_I* } if "`nboot'" != "0" { cap drop BootC_I* BootUnc_I* } } if "`nboot'" == "0" & "`basic'"!="" { noi disp as error "Please specify number of bootstraps" error 198 } if "`percentile'"!="" & "`basic'"!="" | "`nperm'" != "0" & "`nboot'" != "0" { noi disp as error "you have included resampling options that cannot be specified at the same time" error 198 } local ci percentile if "`basic'"!="" local ci `basic' tempfile orig save `orig' local intervraw: copy local intervention fvrevar `intervention', list local intervention `r(varlist)' cap { local newvarlist foreach var of local varlist { if "`var'" != "`intervention'" & !regexm("`var'","i[^\.]*\.`intervention'$") & !regexm("`var'","\((#[0-9])\)*\.`intervention'$") | regexm("`var'","#") { local newvarlist `newvarlist' `var' } else { noi disp as txt "Note: Inclusion of the intervention variable in the variable list is redundant." } } local varlist: copy local newvarlist } fvrevar `varlist', list local varlist_clean `r(varlist)' marksample touse markout `touse' `intervention' keep if `touse' tempvar r brkn_fctor tempname test1 test0 Beta Beta1 Beta2 b Sigma1 Sigma2 Sigma x max nt nc N_total colnumber coef gettoken depvar indepvar: varlist tab `intervention' scalar `max' = r(r) baseset, max(`max') intervention(`intervraw') local refcat `r(refcat)' levelsof `intervention', local(levels) tokenize `levels' tempfile srt save `srt' tempfile srt2 save `srt2' reg `depvar', `options' scalar `Sigma2' = e(rmse) `isily' reg `depvar' i.`intervention' `indepvar' , `options' scalar `Sigma1' = e(rmse) matrix `Beta' = e(b) matrix rownames `Beta' = "Beta" matrix `coef' = `Beta'[1,1..`=`max''] mata funczero("`coef'") matrix `coef'=`coef'' matrix `Sigma' = round(`Sigma1'^2,.01),round(`Sigma2'^2,.01) matrix rownames `Sigma' = "Sigma2" matrix colnames `Sigma' = "Conditional" "Unconditional" scalar `colnumber' = colnumb(r(table),"`=`refcat'+0'b.`intervention'") if "`=`colnumber''"=="1" { matrix `test0' = r(table)[1...,2...] } else { matrix `test0' = r(table)[1...,1..`=`colnumber'-1'],r(table)[1...,`=`colnumber'+1'...] } matrix `test1' = `test0' forvalues i = 1/`=rowsof(`test0')' { forvalues j = 1/`=colsof(`test0')' { matrix `test1'[`i',`j']= round(`test0'[`i',`j'],.01) } } matrix `Beta1' = (`test1'["b", .] \ `test1'["ll".."ul", . ]) matrix `Beta2'=`Beta1'' matrix colnames `Beta2' = "Estimate" "95% LB" "95% UB" forvalues i = 1/`=`max''{ if "`=`refcat'+0'" != "``i''" { local rowname `rowname' "`intervention'``i''" local fnc `fnc' `i' } else { local m `i' } } tokenize `fnc' qui tab `intervention', generate(`brkn_fctor') matcell(`x') /* extract matrix x of broken Intervention to use in calculations below*/ drop `brkn_fctor'1 clear forvalues s = 1/2 { forvalues j = 1/`=`max'-1' { tempname cd`j'`s' varcd`j'`s' secd`j'`s' cdlb`j'`s' cdub`j'`s' jdf g`j'`s' varg`j'`s' seg`j'`s' glb`j'`s' gub`j'`s' A`j'`s' g`j'`s' scalar `cd`j'`s'' = (`coef'[`j',1]/ `Sigma`s'' ) scalar `nt' = `x'[`m',1] scalar `nc' = `x'[``j'',1] scalar `varcd`j'`s''= ((`nt'+`nc')/(`nt'*`nc')+`cd`j'`s''^2/(2*(`nt'+`nc'))) scalar `secd`j'`s'' = sqrt(`varcd`j'`s'') scalar `cdlb`j'`s'' = (`cd`j'`s'' - 1.96*`secd`j'`s'') scalar `cdub`j'`s'' = (`cd`j'`s'' + 1.96*`secd`j'`s'') scalar `jdf' = (1 - (3/(4*(`nt'+`nc'-2)-1))) scalar `g`j'`s'' = (`jdf'*`cd`j'`s'') scalar `varg`j'`s'' = (`jdf'^2 * `varcd`j'`s'') scalar `seg`j'`s'' = sqrt(`varg`j'`s'') scalar `glb`j'`s'' = (`g`j'`s'' - 1.96*`seg`j'`s'') scalar `gub`j'`s'' = (`g`j'`s'' + 1.96*`seg`j'`s'') matrix `A`j'`s'' = round(`g`j'`s'',.01),round(`glb`j'`s'',.01),round(`gub`j'`s'',.01) } } forvalues s = 1/2 { tempname G`s' P`s' matrix `G`s'' = `A1`s'' matrix `P`s'' = round(`g1`s'',.01) } forvalues s = 1/2 { if `=`max'-1'>1 { matrix `G`s'' = `A1`s'' foreach i of numlist 2/`=`max'-1'{ /*extract g's (coefficients of the long calculation to be used later)*/ matrix `G`s'' = `G`s'' \ `A`i'`s'' } matrix `P`s'' = round(`g1`s'',.01) foreach i of numlist 2/`=`max'-1'{ matrix `P`s'' = (`P`s'' \ round(`g`i'`s'',.01)) } } } matrix CondES = `G1' matrix rownames CondES = `rowname' matrix colnames CondES = "Estimate" "95% LB" "95% UB" matrix UncondES = `G2' matrix rownames UncondES = `rowname' matrix colnames UncondES = "Estimate" "95% LB" "95% UB" use `srt' //====================================================// //=================== =================// //================== PERMUTATIONS ==================// //================= ===================// //====================================================// if "`nperm'" != "0" { tempvar shuffle noisily di as txt " Running Permutations..." set seed `seed' if `nperm'<1000 { display as error "error: nPerm must be greater than 1000" error 7 } forvalues k = 1/`nperm' { if "`dot'" == "" { if !mod(`k', 100) { noi di _c "`k'" } else { if !mod(`k', 10) { noi di _c "." } } } keep `intervention' gen `shuffle' = runiform() mata funcsh("`intervention'","`shuffle'") drop `shuffle' tempfile permute save `permute' use `srt2' merge 1:1 _n using `permute', update replace nogenerate tempfile srt2 save `srt2' gettoken depvar indepvar: varlist reg `depvar', `options' scalar `Sigma2' = e(rmse) `isily' reg `depvar' i.`intervention' `indepvar' , `options' scalar `Sigma1' = e(rmse) matrix `b' = e(b) capture drop `coef' matrix `coef' = `b'[1,1..`=`max''] mata funczero("`coef'") matrix `coef'=`coef'' qui tab `intervention', generate(`brkn_fctor') matcell(`x') drop `brkn_fctor'1 clear forvalues s = 1/2 { forvalues j = 1/`=`max'-1' { tempname c`k'd`j'`s' var`k'cd`j'`s' g`k'`j'`s' scalar `c`k'd`j'`s'' = (`coef'[`j',1]/ `Sigma`s'' ) scalar `nt' = `x'[`m',1] scalar `nc' = `x'[``j'',1] scalar `var`k'cd`j'`s''= ((`nt'+`nc')/(`nt'*`nc')+`c`k'd`j'`s''^2/(2*(`nt'+`nc'))) scalar `jdf' = (1 - (3/(4*(`nt'+`nc'-2)-1))) scalar `g`k'`j'`s'' = (`jdf'*`c`k'd`j'`s'') } } use `srt2' } clear use `srt' count scalar `N_total' = `r(N)' if `nperm'>`=`N_total'' { set obs `nperm' } forvalues k = 1/`nperm' { forvalues s = 1/2 { forvalues i = 1/`=`max'-1' { capture gen Treat_PEst`s'_`i' = . replace Treat_PEst`s'_`i' = `g`k'`i'`s'' in `k' } } } /* Permutation Test*/ tempname pval matrix `pval'=J(2,`=`max'-1',.) matrix rownames `pval' = "Conditional ES" "Unconditional ES" matrix colnames `pval' = `rowname' forvalues s=1/2 { forvalues i=1/`=`max'-1' { tempvar p`s'_`i' gen `p`s'_`i'' = abs(Treat_PEst`s'_`i')>=abs(`cd`i'`s'') summarize `p`s'_`i'', meanonly matrix `pval'[`s',`i'] =round(r(mean),.01) } } return matrix Pv = `pval' tokenize `levels' if "`paste'"!="" { local m forvalues i = 1/`=`max'' { if "`=`refcat'+0'" != "``i''" { local m = `m' + 1 rename (Treat_PEst1_`m' Treat_PEst2_`m') (PermC_I``i'' PermUnc_I``i'') } } keep PermC_I* PermUnc_I* tempfile perMES save `perMES' use `orig' merge 1:1 _n using `perMES', nogenerate tempfile orig save `orig' } if "`dot'" == "" { noi di as txt "" } noi di as txt " Permutations completed." } /*if nperm*/ //====================================================// //==================== ====================// //==================== BOOTSTRAPS ====================// //==================== ====================// //===================================================// if "`nboot'" != "0" { clear use `srt' count cap scalar `N_total' = `r(N)' noisily di as txt " Running Bootstraps..." if `nboot'<1000 { display as error "error: nBoot must be greater than 1000" error 7 } keep `varlist_clean' `intervention' tempfile touseit save `touseit' set seed `seed' forvalues i = 1/`nboot' { if "`dot'" == "" { if !mod(`i', 100) { noi di _c "`i'" } else { if !mod(`i', 10) { noi di _c "." } } } bsample keep `varlist_clean' `intervention' gettoken depvar indepvar: varlist reg `depvar', `options' scalar `Sigma2' = e(rmse) `isily' reg `depvar' i.`intervention' `indepvar', `options' scalar `Sigma1' = e(rmse) matrix `b' = e(b) capture drop `coef' matrix `coef' = `b'[1,1..`=`max''] mata funczero("`coef'") matrix `coef'=`coef'' forvalues s = 1/2 { forvalues j = 1/`=`max'-1' { tempvar Treat`j'_RM`i'SE`s' scalar `Treat`j'_RM`i'SE`s'' = (`coef'[`j',1]/ `Sigma`s'' ) } } clear use `touseit' } if `nboot'>`=`N_total'' { set obs `nboot' } forvalues i = 1/`nboot' { forvalues s = 1/2 { forvalues j = 1/`=`max'-1' { cap gen Treat`j'_RMSE`s' = . replace Treat`j'_RMSE`s' = `Treat`j'_RM`i'SE`s'' in `i' } } } forvalues s=1/2 { if "`ci'"=="percentile" { forvalues i=1/`=`max'-1' { centile Treat`i'_RMSE`s', centile(2.5) tempname x_`i'`s' y_`i'`s' scalar `x_`i'`s''=r(c_1) centile Treat`i'_RMSE`s', centile(97.5) scalar `y_`i'`s''=r(c_1) } } if "`ci'"=="basic" { forvalues i=1/`=`max'-1' { centile Treat`i'_RMSE`s', centile(2.5) tempname x_`i'`s' y_`i'`s' scalar `y_`i'`s''=2*`cd`i'`s''-r(c_1) centile Treat`i'_RMSE`s', centile(97.5) scalar `x_`i'`s''=2*`cd`i'`s''-r(c_1) } } tempname L1`s' L2`s' L`s' matrix `L1`s''=round(`x_1`s'',.01) matrix `L2`s''=round(`y_1`s'',.01) if `=`max'-1'>1 { foreach i of numlist 2/`=`max'-1' { matrix `L1`s''=`L1`s'' \ round(`x_`i'`s'',.01) } foreach i of numlist 2/`=`max'-1' { matrix `L2`s''=`L2`s'' \ round(`y_`i'`s'',.01) } } matrix `L`s'' = (`P`s'',`L1`s'',`L2`s'') matrix rownames `L`s'' = `rowname' matrix colnames `L`s'' = "Estimate" "95% LB (BT)" "95% UB (BT)" } tokenize `levels' if "`paste'"!="" { local m forvalues i = 1/`=`max'' { if "`=`refcat'+0'" != "``i''" { local m = `m' + 1 rename (Treat`m'_RMSE1 Treat`m'_RMSE2) (BootC_I``i'' BootUnc_I``i'') } } keep BootC_I* BootUnc_I* tempfile origBoot save `origBoot' clear use `orig' merge 1:1 _n using `origBoot', nogenerate tempfile orig save `orig' } if "`dot'" == "" { noi di as txt "" } noi di as txt " Bootstraps completed." matrix CondES = `L1' matrix UncondES = `L2' }/*if nboot*/ cap noi { noisily { return matrix Beta = `Beta2' return matrix Sigma2 = `Sigma' matrix list CondES return matrix CondES = CondES matrix list UncondES return matrix UncondES = UncondES } } clear use `orig' restore, not } end capture program drop baseset program define baseset, rclass syntax, max(name) INTervention(varlist fv) if regexm("`intervention'", "bn\.") | regexm("`intervention'", "^i\(?([0-9] ?)+\)?\.") { noi disp as error "i(numlist) not allowed; you must specify a base level" error 198 } local refcat if regexm("`intervention'", "([0-9]?)[ ]*\.") local refcat = regexs(1) local allow opt1 if "`refcat'" == "" { if regexm("`intervention'", "\(\#*([0-9]?)\)[ ]*\.") local refcat = regexs(1) if "`refcat'"!="" local allow opt2 } if "`refcat'" == "" { if regexm("`intervention'", "\(([a-zA-Z]+)*\)\.") local refcat = regexs(1) if "`refcat'"!="" local allow opt3 } fvrevar `intervention', list local intervention `r(varlist)' levelsof `intervention', local(levels) tokenize `levels' tempname min Max scalar `min'=`1' scalar `Max' = ``=`max''' if "`allow'" == "opt1" { forvalues i=1/`=`max'' { if "`refcat'" != "``i''" local s = `s'+1 /*checking cases were intervention is irregular (i.e. 1 4 9) and user has specified a number of baseline that is not 1,4 or 9 and not below 1 a > nd not above 9*/ } if "`refcat'" != "" { if "`refcat'">"`=`Max''" | "`refcat'"<"`=`min''" | "`s'" == "`=`max''" { noi disp as error "{bf:Warning:} selected baseline level `refcat' is out of bounds; level `=`Max'' chosen instead." } } else { local refcat = `=`min'' } if "`s'" == "`=`max''" & "`refcat'" != "`=`min''" { local refcat = `=`Max'' } if "`refcat'" != "" { fvset base `refcat' `intervention' } } else if "`allow'"=="opt2" { fvset base ``refcat'' `intervention' local refcat = ``refcat'' } else if "`allow'"=="opt3" { fvset base `refcat' `intervention' if strpos("`refcat'","first") >0 { local refcat = `=`min'' } if strpos("`refcat'","last") >0 { local refcat = `=`Max'' } if strpos("`refcat'","freq")>0 { tempname maximum z tab `intervention', matcell(`maximum') mata: st_local("matr", strofreal(max(st_matrix("`maximum'")))) forvalues i = 1/`=`max''{ scalar `z' = `maximum'[`i',1] if "`matr'"== "`=`z''" local refcat = ``i'' } } } return local refcat = `refcat' end