/************************************************************************** Name: {~/Documents/Research/NejmOregonLetter/Programs/sim1T.do} Description: {This program runs simulations calculating power under various assumptions of effect sizes, and calculating the sample size required for 80% power to detect the effect sizes.} Input: {} Output: {} External Macros: {none} Run under: {Stata} Modification History: {11May13; SamR; Created. 13May13; SamR; Changed N to be the effective sample size given weights, recalculated outcome means among those who do not have Medicaid, looked at more effect sizes and increased number of reps. 14May13; SamR; Major update: added calculations of required sample sizes, used global macros to save results. 16May13; SamR; Added additional comments.} **************************************************************************/ version 12 #delimit ; capture log close; clear; set more off; log using ~/Documents/Research/NejmOregonLetter/Programs/sim1T.log, replace; * This sets the list of effect sizes we're testing in percent reductions; global effectsizes = "10 20 30 40 50"; global reps = 10000; global outcomes = "highBP highChol lowHDL highHbA1c depressed"; set seed 12345; /************************************************************************** N is the number of observations and will be the first argument after calling the program. For the actual study, N is 12,229, but the effective sample size given analytic weights is 10,192. **************************************************************************/ cap program drop sim; program define sim, rclass; args N outcome effectsize; drop _all; * Set significance level of p = .05; local alpha = .05; /************************************************************************** Nlose is the number in the sample who lost the lottery. We assume the proportion of losers in the sample is non-stochastic. **************************************************************************/ local Nlose = round(`N' * 4786/10192); local Nlose_plus1 = `Nlose' + 1; /************************************************************************** Here we define macros for the means among those without Medicaid on the following outcomes: highBP = 1 if blood pressure elevated highChol = 1 if high total cholesterol lowHDL = 1 if low HDL cholesterol highHbA1c = 1 if elevated HbA1c depressed = 1 if positive screen for depression We calculate pre-Medicaid means using control group means taken from table 2 in the main paper (page 1716-7), inflated by .185 * effectsize because 18.5 percent of controls had Medicaid (and where here the effect size is the point estimate from table 2 in the paper). These means are only slightly higher than the control group means from the paper; if we did not make this adjustment, the calculated power would be very slightly lower. **************************************************************************/ local highBP_mean = .16546; local highChol_mean = .14550; local lowHDL_mean = .28522; local highHbA1c_mean = .052721; local depressed_mean = .31693; /************************************************************************** z is an indicator for winning the lottery. **************************************************************************/ set obs `N'; gen z = 0 in 1/`Nlose'; replace z = 1 in `Nlose_plus1'/`N'; /************************************************************************** x is an indicator for enrollment in Medicaid. Probabilities conditional on lottery result are from supplemental materials table S9 (page 34) **************************************************************************/ gen x = rbinomial(1, 0.185) in 1/`Nlose'; replace x = rbinomial(1, 0.4264) in `Nlose_plus1'/`N'; gen `outcome' = .; /************************************************************************** This takes random draws for whether an individual has the given outcome. Those without Medicaid have the probability above, and those with Medicaid have the probability reduced by the postulated effect size. **************************************************************************/ replace `outcome' = rbinomial(1,``outcome'_mean' - `effectsize'/100 * ``outcome'_mean' * x); /************************************************************************** This section runs a regression of the outcome on the lottery result and calculates p values. **************************************************************************/ reg `outcome' z; local p_`outcome'_`effectsize' = 2*ttail(e(df_r),abs(_b[z]/_se[z])); cap return scalar p_`outcome'_`effectsize' = `p_`outcome'_`effectsize''; cap return scalar sig_`outcome'_`effectsize' = `p_`outcome'_`effectsize'' < `alpha'; end; /************************************************************************** This loop calculates power for each outcome and each effect size reported in the letter. **************************************************************************/ foreach outcome in $outcomes {; foreach effectsize in $effectsizes {; simulate pw_`outcome'_`effectsize'=r(sig_`outcome'_`effectsize'), reps($reps) nodots: sim 10192 `outcome' `effectsize'; qui: sum pw_`outcome'_`effectsize'; global pw_`outcome'_`effectsize' = r(mean); macro list pw_`outcome'_`effectsize'; }; }; /************************************************************************** This loop does a grid search to find the lowest sample size that has 80% power to detect each effect size for each outcome. Our search is precise to three significant digits, and we run fewer simulations than above to reduce computing time. **************************************************************************/ global reps = $reps / 10; foreach outcome in $outcomes {; foreach effectsize in $effectsizes {; local N = 1; local increment = 5; while `N' < 10^(`increment'+3) {; local N = max(round(`N' - 10^(`increment'+1), 10^`increment'),10); local sig = 0; while `sig' < 0.8 {; local N = round(`N' + 10^`increment', 10^`increment'); simulate pw_`outcome'_`effectsize' = r(sig_`outcome'_`effectsize'), reps($reps) nodots: sim `N' `outcome' `effectsize'; qui: sum pw_`outcome'_`effectsize'; local sig=r(mean); }; local --increment; }; global N_`outcome'_`effectsize' = `N'; macro list N_`outcome'_`effectsize'; }; }; macro list; log close;