*** Jan Ditzen - June 2013 *** Programs for the simulation Meta Regression. Do file contains three programs, "funnel_regressions" and "build_funnel". The latter starts the former, so "build_funnel" can not be *** used without "funnel_regressions". The last program is "statistics" which performes the different regressions. /* Changelog: 17.02.2015 - added PEESE regression in program statistics added output of mean of t (t_bar_bar) and sd of b 24.02.2015 - removed MST and weeped out unnecessary parts */ capture program drop funnel_regressions program define funnel_regressions, rclass syntax anything , x_mean(string) x_sd(string) e_mean(string) e_sd(string) j(string) /* Options: x_mean: mean with variable x x_sd: standard deviation of variable x e_mean and e_sd: Mean and SD of error term obs: number of observation in each regression Description: Program performs the individual regressions for a funnel. If repeatet m times, a funnel with m regressions is created. */ clear local obsn = round($obs) set obs `obsn' global obs = $obs + 1/`j' local beta = `anything' *Draw x, draw the error term. gen x = rnormal(`x_mean',`x_sd') gen e = rnormal(`e_mean',`e_sd') *DGP (fixed, always the same) gen y = `beta' * x + e * Local for Regression local reg y x *Run the regression without a constant and save the estimated beta and its standard error in r(). reg `reg', noconst return scalar x = _b[x] return scalar x_se = _se[x] return scalar df = e(N) - 1 - e(df_m) end capture program drop build_funnel program define build_funnel , rclass syntax anything , x_mean(string) x_sd(string) e_mean(string) e_sd(string) n(string) j(string) ci(string) obs(string) mint(string) [Save] /* anything: contains beta Options: x_mean: Mean for variable x x_sd: SD for Variable x e_mean: Mean for Error Term e_sd: SD for Error Term N: Number of visible regressions in funnel j: 1/J is the percentage of selected regressions ci: confidence interval for H0: coeff == 0 mint: minimum/smalles acceptable t value, determines point A obs: starting value for number of observations Description: Program performs m regression with the specified parameters. It uses the estimated b's to perform estimations with the right, no and the wrong augmentation. Moreover it computes the mean of the estimated b. The results are stored in r(). If censored, the x_b's to which censoring applies are dropped. */ local beta = `anything' local N = `n' local m = `j'*`N' local J = `j' global obs `obs' #delimit ; simulate x_b = r(x) x_se = r(x_se) df = r(df) , reps(`m') : funnel_regressions `beta', x_mean(`x_mean') x_sd(`x_sd') e_mean(`e_mean') e_sd(`e_sd') j(`j') ; #delimit cr gen t = x_b / x_se gen p = 1 / x_se *** Selection **SR0 select randomly N point estimates preserve egen jset = seq() , block(`J') collapse (mean) t p x_b x_se , by(jset) statistics `beta' S0 `ci' `save' return scalar SR0_bbar = r(bbar) return scalar SR0_b_sd = r(b_sd) return scalar SR0_t_mean = r(t_mean) return scalar SR0_fat_pet_bm = r(fat_pet_bm) return scalar SR0_fat_pet_cbm = r(fat_pet_cbm) return scalar SR0_fat_pet_bf = r(fat_pet_bf) return scalar SR0_fat_pet_cbf = r(fat_pet_cbf) return scalar SR0_peese_bf = r(peese_bf) return scalar SR0_peese_cbf = r(peese_cbf) return scalar SR0_peese_bp = r(peese_bp) return scalar SR0_peese_cbp = r(peese_cbp) restore **SR1 select N point estimates with the largest t preserve egen jset = seq() , block(`J') by jset, sort: egen max = max(t) keep if max == t statistics `beta' S1 `ci' `save' return scalar SR1_bbar = r(bbar) return scalar SR1_b_sd = r(b_sd) return scalar SR1_t_mean = r(t_mean) return scalar SR1_fat_pet_bm = r(fat_pet_bm) return scalar SR1_fat_pet_cbm = r(fat_pet_cbm) return scalar SR1_fat_pet_bf = r(fat_pet_bf) return scalar SR1_fat_pet_cbf = r(fat_pet_cbf) return scalar SR1_peese_bf = r(peese_bf) return scalar SR1_peese_cbf = r(peese_cbf) return scalar SR1_peese_bp = r(peese_bp) return scalar SR1_peese_cbp = r(peese_cbp) restore **SR2 select N point estimates with the largest b preserve egen jset = seq() , block(`J') by jset, sort: egen max = max(x_b) keep if max == x_b statistics `beta' S2 `ci' `save' return scalar SR2_bbar = r(bbar) return scalar SR2_b_sd = r(b_sd) return scalar SR2_t_mean = r(t_mean) return scalar SR2_fat_pet_bm = r(fat_pet_bm) return scalar SR2_fat_pet_cbm = r(fat_pet_cbm) return scalar SR2_fat_pet_bf = r(fat_pet_bf) return scalar SR2_fat_pet_cbf = r(fat_pet_cbf) return scalar SR2_peese_bf = r(peese_bf) return scalar SR2_peese_cbf = r(peese_cbf) return scalar SR2_peese_bp = r(peese_bp) return scalar SR2_peese_cbp = r(peese_cbp) restore **SR3 select N which are most far away. tempfile data save `data' ** Calculate point A: egen jset = seq() , block(`J') gen id = _n gen sample = 0 forvalues i=1(1)`n' { preserve keep if jset == `i' tabstat t x_b, save matrix a = r(StatTotal) local t = abs(a[1,1]) local x_b = abs(a[1,2]) gen q = t/`t' if t/`t' < x_b/`x_b' replace q = x_b/`x_b' if x_b/`x_b' < t/`t' egen max = max(q) keep if max == q sum id local id = r(max) restore replace sample = 1 if id == `id' } keep if sample == 1 statistics `beta' S3 `ci' `save' return scalar SR3_bbar = r(bbar) return scalar SR3_b_sd = r(b_sd) return scalar SR3_t_mean = r(t_mean) return scalar SR3_fat_pet_bm = r(fat_pet_bm) return scalar SR3_fat_pet_cbm = r(fat_pet_cbm) return scalar SR3_fat_pet_bf = r(fat_pet_bf) return scalar SR3_fat_pet_cbf = r(fat_pet_cbf) return scalar SR3_peese_bf = r(peese_bf) return scalar SR3_peese_cbf = r(peese_cbf) return scalar SR3_peese_bp = r(peese_bp) return scalar SR3_peese_cbp = r(peese_cbp) **SR4 select the first which is larger than mint use `data', clear egen jset = seq() , block(`J') gen id = _n gen sample = 0 forvalues i=1(1)`n' { preserve keep if jset == `i' tabstat t x_b, save matrix a = r(StatTotal) local t = abs(a[1,1]) local x_b = abs(a[1,2]) gen q = t/`t' if t/`t' < x_b/`x_b' replace q = x_b/`x_b' if x_b/`x_b' < t/`t' local minb = `x_b' / `t' * `mint' gen qmin = `mint'/`t' if `mint'/`t' < `minb' / `x_b' replace qmin = `minb' / `x_b' if `minb'/`x_b' < `mint'/`t' gen sample1 = 1 if qmin < q sum sample1 local s = r(N) if `s' == 0 { egen max = max(q) keep if max == q sum id local id = r(max) } if `s' != 0 { keep if sample1 == 1 keep if _n == 1 sum id local id = r(max) } restore replace sample = 1 if id == `id' } keep if sample == 1 statistics `beta' S4 `ci' `save' return scalar SR4_bbar = r(bbar) return scalar SR4_b_sd = r(b_sd) return scalar SR4_t_mean = r(t_mean) return scalar SR4_fat_pet_bm = r(fat_pet_bm) return scalar SR4_fat_pet_cbm = r(fat_pet_cbm) return scalar SR4_fat_pet_bf = r(fat_pet_bf) return scalar SR4_fat_pet_cbf = r(fat_pet_cbf) return scalar SR4_peese_bf = r(peese_bf) return scalar SR4_peese_cbf = r(peese_cbf) return scalar SR4_peese_bp = r(peese_bp) return scalar SR4_peese_cbp = r(peese_cbp) end capture program drop statistics program statistics, rclass syntax anything tokenize `anything' local beta = `1' local sr = "`2'" local ci = `3' local save = "`4'" **PET: reg t p return scalar fat_pet_bf = _b[_cons] return scalar fat_pet_bm = _b[p] return scalar bias_bm = _b[p] - `beta' test p == 1 return scalar fat_pet_cbm = r(p)>`ci' test _cons == 0 return scalar fat_pet_cbf = r(p)>`ci' **PEESE: /* Comments for PEESE: regression eq: (Eq 9c from paper): t_i = beta_f * s_i + beta_p * p_i + v_i2 regression of PEESE is on the t-value, explanatory variables are standard error (s_i; in code x_se) and p_i (1/standard error). */ reg t x_se p , noconst return scalar peese_bf = _b[x_se] return scalar peese_bp = _b[p] test x_se == 0 return scalar peese_cbf = r(p)>`ci' test p == 1 return scalar peese_cbp = r(p)>`ci' **Summary Stats (b_bar and b_sd) sum x_b return scalar bbar = r(mean) return scalar b_sd = r(sd) sum t return scalar t_mean = r(mean) ** Save Funnel if "`save'" == "save" { local name = subinstr("`c(current_time)' `c(current_date)'",":","_",.) local name = subinstr("`name'"," ","_",.) save "`c(pwd)'\datasets\\funnels_`name'_`sr'.dta", replace } end