*! ml_benchmark v1.5.0 metaLong for Stata 14.1 *! Benchmark Calibration of Longitudinal ITCV Against Observed Covariates *! Subir Hait, Michigan State University | haitsubi@msu.edu *! *! FIX v1.5.0: replaced preserve/restore inside loop with tempfile approach *! because Stata 14.1 does not allow nested preserve (r(621) already preserved) program define ml_benchmark, rclass version 14.1 syntax varlist(min=2 max=2 numeric) [if] [in] , /// STUDy(varname) /// TIME(varname numeric) /// METAfile(string asis) /// SENSfile(string asis) /// COVariates(varlist numeric) /// [ /// ALpha(real 0.05) /// MINk(integer 3) /// SAVing(string asis) /// REPLACE /// ] tokenize `varlist' local yi `1' local vi `2' /* ---- Mark analytic sample ----------------------------------------- */ marksample touse markout `touse' `time' `yi' `vi' `covariates' capture confirm string variable `study' if !_rc { quietly replace `touse' = 0 if trim(`study') == "" } else { quietly replace `touse' = 0 if missing(`study') } local cov_list `covariates' local n_covs : word count `cov_list' if `n_covs' == 0 { di as error "covariates() must contain at least one numeric covariate." exit 198 } /* ---- Load pooled meta results -------------------------------------- */ preserve quietly use `"`metafile'"', clear capture confirm variable time k theta se df t_stat p_val ci_lb ci_ub tau2 if _rc { di as error "metafile() invalid: expected variables not found." restore exit 198 } tempname meta_mat mkmat time k theta se df t_stat p_val ci_lb ci_ub tau2, matrix(`meta_mat') local n_times = rowsof(`meta_mat') restore /* ---- Load sensitivity results -------------------------------------- */ preserve quietly use `"`sensfile'"', clear capture confirm variable time itcv_alpha if _rc { di as error "sensfile() invalid: expected variables not found." restore exit 198 } tempname sens_mat mkmat time itcv_alpha, matrix(`sens_mat') restore /* ---- Save working data to tempfile --------------------------------- */ /* ONE save/use cycle; NO preserve inside the main loop */ quietly keep if `touse' quietly keep `study' `time' `yi' `vi' `cov_list' tempfile workdata quietly save `workdata' /* ---- Result matrix ------------------------------------------------- */ /* cols: time cov_idx k r_partial t_stat df p_val itcv_alpha beats */ local n_rows = `n_times' * `n_covs' tempname bench_mat matrix `bench_mat' = J(`n_rows', 9, .) matrix colnames `bench_mat' = /// time cov_idx k r_partial t_stat df p_val itcv_alpha beats local mat_row = 0 /* ================================================================== */ /* Main loop: time points x covariates */ /* ================================================================== */ forvalues row = 1/`n_times' { local t_val = `meta_mat'[`row', 1] local tau2_t = max(0, `meta_mat'[`row', 10]) /* ITCV_adj threshold for this time point */ local itcv_t = . forvalues i = 1/`=rowsof(`sens_mat')' { if `sens_mat'[`i', 1] == `t_val' { local itcv_t = `sens_mat'[`i', 2] } } local cov_idx = 0 foreach cov of local cov_list { local cov_idx = `cov_idx' + 1 local mat_row = `mat_row' + 1 /* ---- Load subset from tempfile — no preserve needed ------ */ quietly use `workdata', clear quietly keep if `time' == `t_val' quietly keep `study' `yi' `vi' `cov' quietly drop if missing(`yi') | missing(`vi') | missing(`cov') quietly count local n_use = r(N) /* Populate defaults (missing = skipped) */ matrix `bench_mat'[`mat_row', 1] = `t_val' matrix `bench_mat'[`mat_row', 2] = `cov_idx' matrix `bench_mat'[`mat_row', 3] = `n_use' matrix `bench_mat'[`mat_row', 8] = `itcv_t' if `n_use' < `mink' continue /* ---- Check covariate has non-zero variance --------------- */ /* summarize WITHOUT meanonly so that r(Var) is returned */ quietly summarize `cov' if missing(r(Var)) | r(Var) < 1e-15 continue local cov_mean = r(mean) /* ---- Build RE weights and centred covariate -------------- */ tempvar cov_c w quietly gen double `cov_c' = `cov' - `cov_mean' quietly gen double `w' = 1 / (`vi' + `tau2_t') quietly count if missing(`w') | `w' <= 0 if r(N) > 0 continue /* ---- Weighted meta-regression ---------------------------- */ /* One observation per study within this time slice: */ /* ordinary WLS suffices — no clustering needed */ capture quietly regress `yi' `cov_c' [aw=`w'] if _rc continue capture local t_cov = _b[`cov_c'] / _se[`cov_c'] if _rc continue local df_cov = `n_use' - 2 if `df_cov' < 1 local df_cov = 1 local p_cov = 2 * ttail(`df_cov', abs(`t_cov')) local r_par = `t_cov' / sqrt(`t_cov'^2 + `df_cov') local beats = . if !missing(`itcv_t') & !missing(`r_par') { local beats = (abs(`r_par') >= `itcv_t') } matrix `bench_mat'[`mat_row', 4] = `r_par' matrix `bench_mat'[`mat_row', 5] = `t_cov' matrix `bench_mat'[`mat_row', 6] = `df_cov' matrix `bench_mat'[`mat_row', 7] = `p_cov' matrix `bench_mat'[`mat_row', 9] = `beats' } /* end covariate loop */ } /* end time loop */ /* ================================================================== */ /* Build output dataset */ /* ================================================================== */ clear svmat double `bench_mat', names(col) gen str32 covariate = "" local i = 0 foreach cov of local cov_list { local i = `i' + 1 quietly replace covariate = "`cov'" if cov_idx == `i' } drop cov_idx order time covariate k r_partial t_stat df p_val itcv_alpha beats label var time "Follow-up time" label var covariate "Covariate name" label var k "Number of studies" label var r_partial "Partial correlation with effect" label var t_stat "t-statistic (covariate slope)" label var df "Degrees of freedom" label var p_val "p-value" label var itcv_alpha "ITCV_adj threshold at this time" label var beats "1 if |r_partial| >= ITCV_adj" char _dta[ml_type] "ml_benchmark" char _dta[ml_alpha] "`alpha'" char _dta[ml_mink] "`mink'" /* ================================================================== */ /* Display */ /* ================================================================== */ di _newline di as txt " {hline 72}" di as txt " metaLong: Benchmark Calibration of ITCV" di as txt " Covariates: `cov_list'" di as txt " {hline 72}" di as txt _col(4) "Time" _col(10) "Covariate" /// _col(26) "k" _col(30) "|r_partial|" /// _col(44) "ITCV_adj" _col(55) "Beats?" /// _col(63) "p-val" di as txt " {hline 72}" quietly count local nout = r(N) forvalues j = 1/`nout' { local t_r = time[`j'] local c_r = covariate[`j'] local k_r = k[`j'] local rp_r = r_partial[`j'] local ia_r = itcv_alpha[`j'] local bt_r = beats[`j'] local pv_r = p_val[`j'] if missing(`rp_r') { di as txt %7.1f `t_r' " " %-15s "`c_r'" %4.0f `k_r' /// " (skipped)" } else { local bts_lbl = cond(`bt_r'==1, "YES", "no ") di as txt %7.1f `t_r' " " %-15s "`c_r'" as res %4.0f `k_r' /// %12.4f abs(`rp_r') %12.4f `ia_r' " `bts_lbl'" /// %10.4f `pv_r' } } di as txt " {hline 72}" di as txt " 'Beats': |r_partial| >= ITCV_adj means observed confounding" di as txt " of that magnitude would be sufficient to nullify the effect." /* ================================================================== */ /* Save */ /* ================================================================== */ if `"`saving'"' != "" { if "`replace'" != "" quietly save `saving', replace else quietly save `saving' di as txt _newline " Results saved to: " as res `"`saving'"' } return matrix bench = `bench_mat' return scalar alpha = `alpha' return scalar mink = `mink' end