#delim ; prog def regpar, rclass; version 14.0; /* Estimate logit scenario proportions and z-transformed population attributable risk from existing estimation results assumed to contain parameters of a model whose predicted values are proportions bounded between 0 and 1, and calculate confidence intervals for the scenario proportions and the untransformed population attributable risk. *|Author: Roger Newson *!Date: 03 September 2015. */ syntax [if] [in] [pweight aweight fweight iweight], [ , Level(cilevel) post * ]; if "`post'"=="" {; tempname oldest; cap estimates store `oldest'; }; * Create estimation results *; tempname cimat; _regpar `if' `in' [`weight'`exp'] , `options' level(`level') cimatrix(`cimat'); * Copy e() results to r() *; return scalar level=`level'; local mscalars: e(scalars); local nmscalar: word count `mscalars'; forv i1=`nmscalar'(-1)1 {; tempname ms_`i1'; local ename: word `i1' of `mscalars'; scal `ms_`i1''=e(`ename'); if !missing(`ms_`i1'') {; return scalar `ename'=`ms_`i1''; }; }; tempname btemp Vtemp; matr def `btemp'=e(b); matr def `Vtemp'=e(V); return matrix V=`Vtemp'; return matrix b=`btemp'; return matrix cimat=`cimat'; return local atspec `"`e(atspec)'"'; return local atzero `"`e(atzero)'"'; if "`post'"=="" {; cap estimates restore `oldest'; }; end; prog def _regpar, eclass; version 14.0; /* Estimate logit scenario proportions and z-transformed population attributable risk from existing estimation results assumed to contain parameters of a model whose predicted values are conditional proportions bounded between 0 and 1, and calculate confidence intervals for the population proportions and the untransformed population attributable risk. */ * Find last estimation command *; local cmd "`e(cmd)'"; if `"`cmd'"'=="" {;error 301;}; local options "Level(cilevel) CImatrix(name)"; if "`cmd'"=="regpar" {; * Replay old estimation results *; syntax [, `options']; }; else {; * Create new estimation results *; syntax [if] [in] [pweight aweight fweight iweight], [ ATspec(string asis) ATZero(string asis) `options' subpop(passthru) PRedict(passthru) vce(passthru) df(passthru) noEsample force ITERate(passthru) ]; * Assign default atspec() and atzero() options and check that each one returns only 1 scenario *; foreach AO in atzero atspec {; cap _ms_at_parse ``AO'', asobserved; if _rc {; disp _n as error `"Invalid at-specification - `AO'(``AO'')"'; error 498; }; else {; cap conf matrix r(at); if !_rc {; local AOrows=rowsof(r(at)); if `AOrows'>1 {; disp _n as error `"Option `AO'(``AO'') returns `AOrows' scenarios"' _n "Only 1 scenario allowed"; error 498; }; if `"``AO''"'=="" {; local `AO' "(asobserved) _all"; }; }; }; }; if `"`atspec'"'=="" & `"`atzero'"'=="" {; disp as error "Scenarios cannot be compared after a constant-only model"; error 498; }; marksample touse; qui margins if `touse' [`weight'`exp'], at(`atzero') at(`atspec') `subpop' `predict' `vce' `df' `esample' `force' post; * Collect scalar e-results from margins *; local mscalars: e(scalars); local i1=0; foreach ename in `mscalars' {; local i1=`i1'+1; tempname ms_`i1'; scal `ms_`i1''=e(`ename'); }; qui nlcom ("Scenario_0":logit(_b[1._at])) ("Scenario_1":logit(_b[2._at])) ("PAR":atanh(_b[1._at]-_b[2._at])), `iterate' `df' post; * Post scalar e-results from margins *; local nmscalar: word count `mscalars'; forv i1=`nmscalar'(-1)1 {; local ename: word `i1' of `mscalars'; if !missing(`ms_`i1'') {; ereturn scalar `ename'=`ms_`i1''; }; }; * Post local e-results *; ereturn local atspec `"`atspec'"'; ereturn local atzero `"`atzero'"'; ereturn local predict "regpar_p"; ereturn local cmdline `""'; ereturn local cmd "regpar"; }; * Display estimation results *; disp as text "Scenario 0: " as result `"`atzero'"' _n as text "Scenario 1: " as result `"`atspec'"' _n as text "Symmetric confidence intervals for the logit proportions" _n "under Scenario 0 and Scenario 1" _n as text "and for the z-transformed population attributable risk (PAR)"; disp as text "Total number of observations used: " as result e(N); ereturn display, level(`level'); * Calculate CI matrix *; if "`cimatrix'"=="" {; tempname cimatrix; }; * Define multiplier for creation of confidence intervals *; tempname mult clfloat; scal `clfloat'=`level'/100; if !missing(e(df_r)) {; * Student's t-distribution *; local dof=e(df_r); scal `mult'=invttail(`dof',0.5*(1-`clfloat')); }; else {; * Normal distribution *; scal `mult'=invnormal(0.5*(1+`clfloat')); }; * Extract estimates amd standard errors for untransformed parameters *; tempname estmat varmat; matr def `estmat'=e(b); matr def `varmat'=e(V); matr def `cimatrix'=J(3,3,.); matr rownames `cimatrix'="Scenario_0" "Scenario_1" "PAR"; matr colnames `cimatrix'="Estimate" "Minimum" "Maximum"; tempname estscal lb ub hwid; * Untransformed scenario proportions *; forv i1=1(1)2 {; scal `estscal'=`estmat'[1,`i1']; scal `hwid'=`varmat'[`i1',`i1']; scal `hwid'=sqrt(`hwid')*`mult'; scal `lb'=`estscal'-`hwid'; scal `ub'=`estscal'+`hwid'; foreach Y in `estscal' `lb' `ub' {; scal `Y'=invlogit(`Y'); }; matr def `cimatrix'[`i1',1]=`estscal'; matr def `cimatrix'[`i1',2]=`lb'; matr def `cimatrix'[`i1',3]=`ub'; }; * Untransformed PAR *; scal `estscal'=`estmat'[1,3]; scal `hwid'=`varmat'[3,3]; scal `hwid'=sqrt(`hwid')*`mult'; scal `lb'=`estscal'-`hwid'; scal `ub'=`estscal'+`hwid'; foreach Y in `estscal' `lb' `ub' {; scal `Y'=tanh(`Y'); }; matr def `cimatrix'[3,1]=`estscal'; matr def `cimatrix'[3,2]=`lb'; matr def `cimatrix'[3,3]=`ub'; * Display CI matrix *; disp _n as text "Asymmetric `level'% CIs for the untransformed proportions" _n "under Scenario 0 and Scenario 1" _n "and for the untransformed population attributable risk (PAR)"; matlist `cimatrix', noheader noblank nohalf lines(none) names(all) format(%9.0g); end;