#delim ; prog def punafcc, rclass; version 14.0; /* Estimate log scenario rate ratios and log population unattributable fractions from existing estimation results assumed to contain parameters of a model whose predicted values are case-control or survival rate ratios that can be logged, and calculate a confidence interval for the population attributable fraction. *! Author: Roger Newson *! Date: 05 February 2023 */ syntax [if] [in] [pweight aweight fweight iweight], [ , Level(cilevel) post * ]; if "`post'"=="" {; tempname oldest; cap estimates store `oldest'; }; * Create estimation results *; tempname cimat; _punafcc `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 exp_atspec `"`e(exp_atspec)'"'; return local exp_atzero `"`e(exp_atzero)'"'; return local atspec `"`e(atspec)'"'; return local atzero `"`e(atzero)'"'; if "`post'"=="" {; cap estimates restore `oldest'; }; end; prog def _punafcc, eclass; version 14.0; /* Estimate log scenario rate ratios and log population unattributable fractions from existing estimation results assumed to contain parameters of a model whose predicted values are conditional rate ratios that can be logged, and calculate a confidence interval for the population attributable fraction. */ * Find last estimation command *; local cmd "`e(cmd)'"; if `"`cmd'"'=="" {;error 301;}; local options "EForm Level(cilevel) CImatrix(name)"; if "`cmd'"=="punafcc" {; * Replay old estimation results *; syntax [, `options']; }; else {; * Create new estimation results *; syntax [if] [in] [pweight aweight fweight iweight], [ ATspec(string asis) subpop(string asis) `options' vce(passthru) df(passthru) noEsample force ITERate(passthru) ]; * Assign default atspec() and atzero() options and check that each one returns only 1 scenario and only specifies as-observed or a constant value for a factor or variable *; 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; }; local statlist `"`r(statlist)'"'; foreach STAT in `statlist' {; if !inlist("`STAT'","asobserved","value","values") {; disp _n as error `"Option `AO'(``AO'') specifies disallowed statistic:"' _n as error `"`STAT'"'; 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; * Create subpopulation indicator variable identifying observations in the input subpop() subpopulation which are also cases or deaths *; tempvar subpopind subpopind1 subpopind2; * Generate subpopind1 from case or death variable *; if inlist(`"`e(cmd2)'"',"stcox","stcrreg") | (`"`e(cmd2)'"'=="" & inlist(`"`e(cmd)'"',"stcox","stcrreg")) {; qui gene `subpopind1'=_d if `touse'; }; else {; cap conf numeric var `e(depvar)'; if _rc {; disp _n as error "Valid dependent variable e(depvar) not available"; error 498; }; qui gene `subpopind1'=`e(depvar)' if `touse'; }; qui compress `subpopind1'; * Generate subpopind2 from supplied subpop() option *; local subpopvar: word 1 of `subpop'; cap conf numeric var `subpopvar'; if _rc {; local subpopopt `"1 `subpop'"'; }; else {; local subpopopt `"`subpop'"'; }; cap gene `subpopind2'=`subpopopt'; if _rc {; disp _n as error "Invalid subpop() option"; error 498; }; qui compress `subpopind2'; * Create subpopind from subpopind1 and subpopind2 *; qui gene byte `subpopind'=`subpopind1' & `subpopind2' if !missing(`subpopind1') & !missing(`subpopind2'); qui compress `subpopind'; * Create expression() options for at-specifications for margins *; tempname atmat bscenmat scenval; matr def `bscenmat'=e(b); local BScols=colsof(`bscenmat'); local BSnames: colnames `bscenmat'; foreach AO in atzero atspec {; _ms_at_parse ``AO'', asobserved; matr def `atmat'=r(at); local AOcols=colsof(`atmat'); local AOnames: colnames `atmat'; local exp_`AO' ""; forv i1=1(1)`BScols' {; * Extract variable names and factor levels *; local BScolcur: word `i1' of `BSnames'; _ms_parse_parts `BScolcur'; if inlist("`r(type)'","interaction","product") {; local BScurk_names=r(k_names); local BScurnames ""; local BScurlevels ""; forv i2=1(1)`BScurk_names' {; local BScurnames `"`BScurnames' `r(name`i2')'"'; local BScurlevels `"`BScurlevels' `r(level`i2')'"'; }; }; else {; local BScurk_names=1; local BScurnames `"`r(name)'"'; local BScurlevels `"`r(level)'"'; }; * Evaluate scenario value of factor varlist element *; scal `scenval'=.; forv i2=1(1)`AOcols' {; local AOnamecur: word `i2' of `AOnames'; _ms_parse_parts `AOnamecur'; local namecur2 "`r(name)'"; local levelcur2 "`r(level)'"; forv i3=1(1)`BScurk_names' {; local namecur3: word `i3' of `BScurnames'; local levelcur3: word `i3' of `BScurlevels'; if (`"`namecur3'"'==`"`namecur2'"') & (`"`levelcur3'"'==`"`levelcur2'"') {; if missing(`scenval') {; scal `scenval'=el(`atmat',1,`i2'); }; else if !missing(el(`atmat',1,`i2')) {; scal `scenval'=`scenval'*el(`atmat',1,`i2'); }; }; }; }; * Update expression if appropriate *; if !missing(`scenval') {; local exp_`AO' `"`exp_`AO'' + _b[`BScolcur']*(`=`scenval''-`BScolcur')"'; }; }; * Finalize expression *; local exp_`AO' `"exp(0`exp_`AO'')"'; }; * Create estimation results *; qui margins if `touse' [`weight'`exp'], exp(`exp_atspec') subpop(`subpopind') `vce' `df' `esample' `force' post; * Collect scalar e-results from margins *; local mscalars: e(scalars); local i1=0; foreach rname in `mscalars' {; local i1=`i1'+1; tempname ms_`i1'; scal `ms_`i1''=e(`rname'); }; qui nlcom ("PUF":log(_b[_cons])), `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 exp_atspec `"`exp_atspec'"'; ereturn local exp_atzero `"`exp_atzero'"'; ereturn local atspec `"`atspec'"'; ereturn local atzero `"`atzero'"'; ereturn local predict "punafcc_p"; ereturn local cmdline `""'; ereturn local cmd "punafcc"; }; * Display estimation results *; if "`eform'"!="" {; local eformopt "eform(Ratio)"; disp as text "Scenario 0: " as result `"`atzero'"' _n as text "Scenario 1: " as result `"`atspec'"' _n as text "Confidence interval for the population unattributable faction (PUF)"; }; else {; local eformopt ""; disp as text "Scenario 0: " as result `"`atzero'"' _n as text "Scenario 1: " as result `"`atspec'"' _n as text "Confidence interval for the log population unattributable faction (PUF)"; }; disp as text "Total number of observations used: " as result e(N); ereturn display, `eformopt' 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 estimate amd standard error for log PUF *; tempname estmat varmat; matr def `estmat'=e(b); matr def `varmat'=e(V); tempname estscal lb ub hwid; scal `estscal'=`estmat'[1,1]; scal `hwid'=`varmat'[1,1]; scal `hwid'=sqrt(`hwid')*`mult'; scal `ub'=`estscal'-`hwid'; scal `lb'=`estscal'+`hwid'; foreach Y in `estscal' `lb' `ub' {; scal `Y'=1-exp(`Y'); }; matr def `cimatrix'=J(1,3,.); matr rownames `cimatrix'="PAF"; matr colnames `cimatrix'="Estimate" "Minimum" "Maximum"; matr def `cimatrix'[1,1]=`estscal'; matr def `cimatrix'[1,2]=`lb'; matr def `cimatrix'[1,3]=`ub'; * Display CI matrix *; disp _n as text "`level'% CI for the population attributable fraction (PAF)"; matlist `cimatrix', noheader noblank nohalf lines(none) names(all) format(%9.0g); end;