program define mulogit,eclass *!Multivariate and univariate odds ratios (logistic regression) plots /* Date: 4/02/08 Last modification: 4/2/08 Author: Leif Peterson, TMHRI, Houston Generates plots of multivariate and univariate odds ratios and write results into the data editor Format: mulogit depvar varlist [if] [in] See help. */ version 10 syntax varlist [if] [in] local touse `if' `in' local ncol:word count `varlist' local depvar:word 1 of `varlist' forv i=2(1)`ncol' { local indepvar:word `i' of `varlist' local indepvars `indepvars' `indepvar' } logit `depvar' `indepvars' `touse' quietly { matrix b=e(b) matrix V=e(V) cap drop typeor gen typeor="" cap drop myvarnames gen myvarnames="" cap drop or gen or=. cap drop ll95ci gen ll95ci=. cap drop ul95ci gen ul95ci=. cap drop test gen test=. cap label drop testlbl label values test testlbl cap drop multtest gen multtest=. cap label drop multtestlbl label values multtest multtestlbl local varnames : colnames(V) local p = colsof(V)-1 forv i =1(1)`p' { local w:word `i' of `varnames' replace myvarnames="`w'" in `i' replace or=exp(b[1,`i']) in `i' replace ll95ci=exp(b[1,`i'] - 1.96 * sqrt(V[`i',`i'])) in `i' replace ul95ci=exp(b[1,`i'] + 1.96 * sqrt(V[`i',`i'])) in `i' replace test=`i' in `i' replace multtest=`i' in `i' label define testlbl `i' "`w'" , add local z = abs(b[1,`i']/sqrt(V[`i',`i'])) if `z' <1.96 { label define multtestlbl `i' "`w'" , modify } else if `z' >=1.96 & `z' <2.326 { label define multtestlbl `i' "*`w'", modify } else if `z' >=2.326 & `z' <3.090 { label define multtestlbl `i' "**`w'" , modify } else if `z' >=3.090 { label define multtestlbl `i' "***`w'" , modify } replace typeor="Multivariate" in `i' } } quietly { cap drop univtest gen univtest=. cap label drop univtestlbl label values univtest univtestlbl } forv i =2(1)`ncol' { quietly { local indepvar:word `i' of `varlist' } logit `depvar' `indepvar' `touse' quietly { matrix b=e(b) matrix V=e(V) local row = `p' + `i' - 1 replace myvarnames="`indepvar'" in `row' replace or=exp(b[1,1]) in `row' replace ll95ci=exp(b[1,1]- 1.96 * sqrt(V[1,1])) in `row' replace ul95ci=exp(b[1,1]+ 1.96 * sqrt(V[1,1])) in `row' local newtest=`i' - 1 replace test=`newtest' in `row' replace univtest=`newtest' in `row' label define testlbl `newtest' "`indepvar'" , modify local z = abs(b[1,1]/sqrt(V[1,1])) if `z' < 1.96 { label define univtestlbl `newtest' "`indepvar'" , modify } else if `z' >=1.96 & `z' <2.326 { label define univtestlbl `newtest' "*`indepvar'" , modify } else if `z' >=2.326 & `z' <3.090 { label define univtestlbl `newtest' "**`indepvar'" , modify } else if `z' >=3.090 { label define univtestlbl `newtest' "***`indepvar'" , modify } replace typeor="Univariate" in `row' } } twoway rcap ul95ci ll95ci test,legend(off) || scatter or test,legend(off) by(typeor) xlabel(#"`p'", labels angle(vertical) valuelabel) note("*p<0.05, **p<0.01, ***p<0.001") saving(multuniv,replace) twoway rcap ul95ci ll95ci multtest if typeor=="Multivariate"|| scatter or multtest if typeor=="Multivariate", legend(off) xlabel(#"`p'", labels angle(vertical) valuelabel) xtitle("") l2title("Multivariate") l1title("OR(95% CI)") note("*p<0.05, **p<0.01, ***p<0.001") saving(mult,replace) twoway rcap ul95ci ll95ci univtest if typeor=="Univariate"|| scatter or univtest if typeor=="Univariate", legend(off) xlabel(#"`p'", labels angle(vertical) valuelabel) xtitle("") l2title("Univariate") l1title("OR(95% CI)") note("*p<0.05, **p<0.01, ***p<0.001") saving(univ,replace) graph combine univ.gph mult.gph, ycommon end