/* Title: A procedure to tabulate and plot results after flexible modeling of a quantitative covariate Authors: Nicola Orsini, Sander Greenland Last update: 13 September 2010 Do-file to reproduce analysis and figures of the paper */ set scheme sj version 11 use http://nicolaorsini.altervista.org/data/pa_luts, clear tabstat tpa, stat(n p25 p50 p75 mean sd) by(ipss2) format(%2.0f) egen cases = total(ipss2==1), by(tpac) egen noncases = total(ipss2==0), by(tpac) generate odds = cases/noncases bysort tpac : generate flag = _n==1 twoway (bar odds tpac, barw(.9) fcolor(white) lc(black) lw(medium)) if flag==1, ylabel(0(.1).9, angle(horiz) format(%2.1fc)) ytitle("Odds (Cases/Noncases) of LUTS") xtitle("Total physical activity, MET-hours/day") xlabel(1(1)10, valuelabel angle(45)) scheme(s1mono) legend(off) name(figure1, replace) table tpac, contents(freq sum ipss2 min tpa median tpa max tpa) sjlog using xblc1, replace generate all = 1 table all, contents(freq p20 tpa p40 tpa p60 tpa p80 tpa) generate tpa2 = tpa^2 generate tpa3 = tpa^3 generate tpap1 = max(0,tpa-37.2)^3 generate tpap2 = max(0,tpa-39.6)^3 generate tpap3 = max(0,tpa-42.3)^3 generate tpap4 = max(0,tpa-45.6)^3 sjlog close, replace sjlog using xblc2, replace logit ipss2 tpa tpa2 tpa3 tpap1 tpap2 tpap3 tpap4 sjlog close, replace testparm tpa2 tpa3 tpap1 tpap2 tpap3 tpap4 sjlog using xblc3, replace xblc tpa tpa2 tpa3 tpap1 tpap2 tpap3 tpap4, covname(tpa) reference(29) eform at(29 32 35 38 40 43 45 48 52 55) generate(pa or lb ub) sjlog close, replace sjlog using xblc4, replace twoway (rcap lb ub pa, sort) (scatter or pa, sort), legend(off) scheme(s1mono) xlabel(29 32 35 38 40 43 45 48 52 55) ylabel(.2(.2)1.2, angle(horiz) format(%2.1fc) ) ytitle("Unadjusted Odds Ratios of LUTS" ) xtitle("Total physical activity, MET-hours/day") name(figure2, replace) sjlog close, replace sjlog using xblc5, replace capture drop pa or lb ub xblc tpa tpa2 tpa3 tpap*, covname(tpa) at(29(1)55) reference(29) eform generate(pa or lb ub) twoway (rcap lb ub pa, sort) (scatter or pa, sort), legend(off) scheme(s1mono) xlabel(29(2)55) xmtick(29(1)55) ylabel(.2(.2)1.2, angle(horiz) format(%2.1fc)) ytitle("Unadjusted Odds Ratios of LUTS") xtitle("Total physical activity, MET-hours/day") name(figure3, replace) sjlog close, replace sjlog using xblc6, replace capture drop pa or lb ub quietly levelsof tpa, local(levels) quietly xblc tpa tpa2 tpa3 tpap*, covname(tpa) at(`r(levels)') reference(29) eform generate(pa or lb ub) twoway (line lb ub pa, sort lc(black black) lp(- -)) (line or pa, sort lc(black) lp(l)) if inrange(pa,29,55), legend(off) scheme(s1mono) xlabel(29(2)55) xmtick(29(1)55) ylabel(.2(.2)1.2, angle(horiz) format(%2.1fc)) ytitle("Unadjusted Odds Ratios of LUTS") xtitle("Total physical activity, MET-hours/day") name(figure4, replace) sjlog close, replace // Left restricted cubic splines sjlog using xblc7, replace testparm tpa2 tpa3 sjlog close, replace sjlog using xblc8, replace logit ipss2 tpa tpap1 tpap2 tpap3 tpap4 sjlog close, replace sjlog using xblc9, replace xblc tpa tpap*, covname(tpa) at(29 32 35 38 40 43 45 48 52 55) reference(29) eform sjlog close, replace // Right restricted cubic splines sjlog using xblc10, replace generate tpan = -tpa generate tpapn1 = max(0,45.6-tpa)^3 generate tpapn2 = max(0,42.3-tpa)^3 generate tpapn3 = max(0,39.6-tpa)^3 generate tpapn4 = max(0,37.2-tpa)^3 logit ipss2 tpan tpapn* sjlog close, replace sjlog using xblc11, replace xblc tpan tpapn*, covname(tpa) at(29 32 35 38 40 43 45 48 52 55) reference(29) eform sjlog close, replace // Both tails restricted cubic splines sjlog using xblc12, replace mkspline tpas = tpa, knots(37.2 39.6 42.3 45.6) cubic sjlog close, replace sjlog using xblc13, replace logit ipss2 tpas1 tpas2 tpas3 sjlog close, replace testparm tpas2 tpas3 sjlog using xblc14, replace xblc tpas*, covname(tpa) at(29 32 35 38 40 43 45 48 52 55) reference(29) eform sjlog close, replace // Figure 5 - Overlay different types of cubic splines * both tails unrestricted logit ipss2 tpa tpa2 tpa3 tpap1 tpap2 tpap3 tpap4 local aic1 : di %6.0f -2*e(ll) + 2*e(k) capture drop pa quietly levelsof tpa, local(levels) quietly xblc tpa tpa2 tpa3 tpap1 tpap2 tpap3 tpap4, covname(tpa) at(`r(levels)') reference(29) eform generate(pa orun lbun ubun) * left restricted logit ipss2 tpa tpap1 tpap2 tpap3 tpap4 local aic2 : di %6.0f -2*e(ll) + 2*e(k) capture drop pa quietly levelsof tpa, local(levels) quietly xblc tpa tpap1 tpap2 tpap3 tpap4, covname(tpa) at(`r(levels)') reference(29) eform generate(pa orlr lblr ublr) * right restricted logit ipss2 tpan tpapn* local aic3 : di %6.0f -2*e(ll) + 2*e(k) capture drop pa quietly levelsof tpa, local(levels) quietly xblc tpan tpapn*, covname(tpa) at(`r(levels)') reference(29) eform generate(pa orrr lbrr ubrr) * both tails restricted logit ipss2 tpas1 tpas2 tpas3 local aic4 : di %6.0f -2*e(ll) + 2*e(k) capture drop pa quietly levelsof tpa, local(levels) quietly xblc tpas1 tpas2 tpas3, covname(tpa) at(`r(levels)') reference(29) eform generate(pa orbr lbbr ubbr) di "`aic1'" di "`aic2'" di "`aic3'" di "`aic4'" tw (line orun orrr orlr orbr pa, sort lc(black black black black) lp(solid shortdash longdash dash_dot )) if inrange(pa,29,55), scheme(s1mono) xlabel(29(2)55) xmtick(29(1)55) ylabel(.2(.2)1.2, angle(horiz) format(%2.1fc)) ytitle("Unadjusted Odds Ratios of LUTS" ) xtitle("Total physical activity, MET-hours/day") legend(label(1 "Both-tail unrestricted AIC = `aic1'") label(2 "Right-tail restricted AIC = `aic3'") label(3 "Left-tail restricted AIC = `aic2'") label(4 "Both-tail restricted AIC = `aic4'") ring(0) pos(1) col(1)) name(figure5, replace) // Adjusting for other covariates sjlog using xblc15, replace quietly summarize age generate agec = age - r(mean) logit ipss2 tpan tpapn* agec sjlog close, replace sjlog using xblc16, replace xblc tpan tpapn*, covname(tpa) at(29 32 35 38 40 43 45 48 52 55) reference(29) eform sjlog close, replace sjlog using xblc17, replace capture drop pa or lb ub quietly levelsof tpa, local(levels) quietly xblc tpan tpapn*, covname(tpa) at(`r(levels)') reference(29) eform generate(pa or lb ub) twoway (line lb ub pa, sort lc(black black) lp(- -)) (line or pa, sort lc(black) lp(l)) if inrange(pa,29,55), legend(off) scheme(s1mono) xlabel(29(2)55) xmtick(29(1)55) ylabel(.2(.2)1.2, angle(horiz) format(%2.1fc)) ytitle("Age-adjusted Odds Ratios of LUTS") xtitle("Total physical activity, MET-hours/day") name(figure6, replace) sjlog close, replace // Uncertainty for the predicted response sjlog using xblc18, replace capture drop pa quietly levelsof tpa, local(levels) quietly xblc tpan tpapn*, covname(tpa) at(`r(levels)') reference(29) eform generate(pa rcc lbo ubo) pr twoway (line lbo ubo pa, sort lc(black black) lp(- -)) (line rcc pa, sort lc(black) lp(l)) if inrange(pa,29,55), legend(off) scheme(s1mono) xlabel(29(2)55) xmtick(29(1)55) ylabel(.2(.1).8, angle(horiz) format(%2.1fc)) ytitle("Age-adjusted Odds (Cases/Noncases) of LUTS") xtitle("Total physical activity, MET-hours/day") name(figure7, replace) sjlog close, replace // Use of xblc command after other modeling approaches (Figure 8) // Categorical model sjlog using xblc19, replace xi:logit ipss2 i.tpac agec, or sjlog close, replace capture drop pa sjlog using xblc20, replace quietly levelsof tpa, local(levels) quietly xblc _Itpac_2- _Itpac_10, covname(tpa) at(`r(levels)') eform generate(pa oddsc lboc uboc) pr sjlog close, replace // Linear splines: a simple alternative sjlog using xblc21, replace generate tpa38p = max(tpa-38, 0) logit ipss2 tpa tpa38p agec xblc tpa tpa38p, covname(tpa) at(29 32 35 38 40 43 45 48 52 55) reference(29) eform sjlog close, replace capture drop pa quietly levelsof tpa, local(levels) quietly xblc tpa tpa38p, covname(tpa) at(`r(levels)') pr eform generate(pa oddls lbls ubls) // Fractional polynomials sjlog using xblc23, replace mfp logit ipss2 tpa agec, df(agec:1) sjlog close, replace sjlog using xblc24, replace xblc Itpa__1 Itpa__2, covname(tpa) at(29 32 35 38 40 43 45 48 52 55) reference(29) eform sjlog close, replace capture drop pa quietly levelsof tpa, local(levels) quietly xblc Itpa__1 Itpa__2, covname(tpa) at(`r(levels)') eform generate(pa oddfp lbfp ubfp) pr // Right restricted cubic spline logit ipss2 tpan tpapn* agec capture drop pa capture drop oddrr capture drop lbrr ubrr quietly levelsof tpa, local(levels) quietly xblc tpan tpapn*, covname(tpa) at(`r(levels)') pr eform generate(pa oddrr lbrr ubrr) // Overlay the different strategies twoway (line oddsc pa, sort lc(black) lp(longdash_dot)) (line oddls pa, sort lc(black) lp(dash)) (line oddfp pa, sort lc(black) lp(longdash)) (line oddrr pa, sort lc(black) lp(l)) if inrange(pa,29,55), scheme(s1mono) xlabel(29(2)55) xmtick(29(1)55) ylabel(.2(.1).8, angle(horiz) format(%2.1fc)) ytitle("Age-adjusted Odds (Cases/Noncases) of LUTS" ) xtitle("Total physical activity, MET-hours/day") legend( label(1 "Categorical model") label(2 "Linear spline model") label(3 "Fractional polynomial model") label(4 "Right-restricted cubic-spline model") ring(0) col(1) pos(1)) name(figure8, replace)