* 	 Paper: Logistic quantile regression in Stata

*  Authors: N.Orsini, M.Bottai

use http://nicolaorsini.altervista.org/data/pa_luts

//  Example given in the text with output omitted 

// Transform the [0,35] bounded outcome in [0,1] bounded outcome

gen ipssb =  (ipss-0)/(35-0) 



// Transform the [0,35] bounded outcome in (0,1) bounded outcome

gen ipssc = [(ipss-0)/(35-0)*(30377-1)+.5]/30377



// Beta regression

xi: betafit  ipssc, mu(i.tpac)



// Fractional logit regression

xi: glm ipssb  i.tpac, family(binomial) link(logit) vce(robust)



tw histogram ipss,  start(0) width(1) percent ///
ylabel(0(5)20, angle(horiz)) ymtick(0(1)20)   xaxis(1 2) ///
xlabel(1 "p25" 3 "p50" 7 "p75" 19 "p95", axis(2) grid gmax labsize(small) ) ///
xtitle("", axis(2)) xlabel(0(5)35, axis(1)) xmtick(0(1)35) fc(white) gap(2)


graph export "lqreg1.eps", as(eps) preview(on) replace



graph box ipss, over(tpac, gap(*1) label(angle(45))) ///
	             nooutsides   ///
				 ylabel(0(5)35, angle(horiz)) ymtick(0(1)35) ///
				 b1title("Total physical activity, MET-hours/day") /// yline(7 19, lw(thick) lp(dash))
				 ytitle(" International prostate symptom score") note("")



graph export "lqreg2.eps", as(eps) preview(on) replace

   

sjlog using lqreg1, replace
xi: lqreg ipss i.tpac, nodots seed(123)
sjlog close, replace   

sjlog using lqreg2, replace
lqregpred crude
table tpac, c(mean crude50) f(%2.0f)
sjlog close, replace 

sjlog using lqreg3, replace
xi: lqreg ipss i.tpac, nodots seed(123)
sjlog close, replace 

sjlog using lqreg4, replace
gen cage = age-59
xi: lqreg ipss i.tpac cage, nodots seed(123)
sjlog close, replace 

sjlog using lqreg5, replace
testparm _I*
sjlog close, replace 

sjlog using lqreg6, replace	
quietly xi: lqreg ipss i.tpac cage, quantiles(25 50 75 95) seed(123) nodots
lqregpred vpadj, for(_Itpac_2- _Itpac_10) at(cage=0)
sjlog close, replace 

sjlog using lqreg6g, replace
twoway  ///
(line vpadj25 tpa, lp(dash) lc(black) sort c(J)) ///
(line vpadj50 tpa, lp(longdash_dot) lc(black) sort c(J) ) ///
(line vpadj75 tpa, lp(longdash) lc(black) sort c(J)) ///
(line vpadj95 tpa, lp(l) lc(black) sort c(J)) /// 
if inrange(tpa,30,60), ///
ytitle("International prostate symptom score") ///
ylabel(0(5)35, angle(horiz)) ymtick(0(1)35) ///
xtitle("Total physical activity, MET-hours/day")   ///
xlabel(30(5)60)  xmtick(30(1)60) ///
legend(label(1 "0.25") label(2 "0.50") label(3 "0.75")  ///
label(4 "0.95") textfirst ring(0) pos(1) col(1) ///
subtitle("Age-Adjusted" "Quantiles") order(4 3 2 1)) 
sjlog close, replace 
	

graph export "lqreg3.eps", as(eps) preview(on) replace

sjlog using lqreg7, replace	
test [q95]_Itpac_10 = [q50]_Itpac_10
sjlog close, replace 


// Modeling tpa and age assuming linearity 

sjlog using lqreg8, replace	
sort ipss tpa
lqreg ipss tpa cage, quantiles(25 50 75 95) nodots seed(123)
sjlog close, replace



sjlog using lqreg9, replace	
bysort ipss tpa: gen flag = _n==1
lqregpred adjl if flag == 1, for(tpa) plotvs(tpa) 
sjlog close, replace

graph export "lqreg5.eps", as(eps) preview(on) replace


sjlog using lqreg10, replace	
lqregplot tpa
sjlog close, replace
graph export "lqreg6.eps", as(eps) preview(on) replace

// Modeling tpa and age with restricted cubic splines

sjlog using lqreg11, replace	
mkspline tpas = tpa, nknots(3) cubic 
mkspline cages = cage, nknots(3) cubic
lqreg ipss tpas1 tpas2 cages1 cages2, quantiles(25 50 75 95) nodots seed(123)
sjlog close, replace

sjlog using lqreg12, replace	
lqregpred adjs if flag == 1, for(tpas1 tpas2) plotvs(tpa)
sjlog close, replace 

graph export "lqreg7.eps", as(eps) preview(on) replace