*! baplot.ado version 1.21 written by PTS (p.seed@umds.ac.uk) (STB-55: sbe33)
*! Produces Bland-Altman plots for two variables
*! See Bland & Altman Lancet Feb 8 1986, pp 307-310
*!
*! syntax: BAanalysis var1 var2 if in, symbol(symbol) format(%6.3f) avlab("Average") difflab("Difference") yline(str) textsize(#) other graph options




cap prog drop BAanalysis
prog define BAanalysis
version 10
     syntax varlist(min=2 max=2 numeric) [if] [in], ///
     [Symbol(string) format(string) avlab(string) difflab(string) ///
     novars noGRaph ci zero mean yline(string) diag saving(string) ///
     text(real 100) ratio(real 1) * ]"

	parse "`varlist'", parse(" ")
	local m1 "`1'"
	local m2 "`2'"


	if "`symbol'" == ""  local symbol "o"
	if "`format'" == "" local format "%6.3f"
	if "`saving'" ~= ""  local saving "saving(`saving')"
	if "`ci'" ~= "" local mean = "mean"


	preserve
	tempvar touse
	mark `touse' `if' `in'
	markout `touse' `m1' `m2'
	qui keep if `touse'

	tempvar av diff
	if "`avlab'" == ""  local avlab "Average"
	if lower("`avlab'") == "nolab"  local avlab " "

	if "`difflab'" == ""  local difflab "Difference"
	if lower("`difflab'") == "nolab"  local difflab " "
	qui gen `av' = (`m1' + `m2')/2
	qui gen `diff' = `m1' - `m2'
	label var `av' "`avlab'"
	label var `diff' "`difflab'"

	qui summ `diff' if `touse'
	local xbar = _result(3)
	local sd = _result(4)^.5
	local n = _result(2)
	local se = `sd'/`n'^.5
	local sig1 = 2*ttail(`n'-2,((`xbar')^2/(`se')^2)^.5)
	local t = -invttail(_result(2)-1, .975)
	local lrr = `xbar' - 1.96*`sd'
	local urr = `xbar' + 1.96*`sd'
	local min = _result(5)
	local max = _result(6)
	local lcb = `xbar' - `t'*`se'
	local ucb = `xbar' + `t'*`se'
      local Se_limits=`sd'*sqrt((1/`n')+1.96*1.96/(2*(`n'-1)))
	local Lower_sd=`sd'*sqrt((`n'- 1)/invchi2((`n'- 1),0.975))
	local Upper_sd=`sd'*sqrt((`n'- 1)/invchi2((`n'- 1),0.025))


	qui summ `av'  if `touse'
	local var_mean = _result(4)
	local var_between = _result(4)- (`sd'^2)/4
	local var_total = _result(4)+ (`sd'^2)/4
	local icc = `var_between '/`var_total '
	local xmin = _result(5)
	local xmax = _result(6)

	local yline "`lrr' `urr'"
	if "`ci'" ~= "" local yline "`yline' `lcb'   `ucb'"
	if "`mean'" ~= "" | "`zero'" == "" local yline "`yline' `xbar'"
	if "`zero'" ~= "" local yline "`yline' 0"

	qui corr `av' `diff'
	local r = _result(4)
	local n = _result(1)
	local sig2 = 2*ttail(`n'-2,( ((`r')^2)*(`n'-2)/(1-(`r')^2))^.5)

	qui spearman `av' `diff'
	local rho = r(rho)
	local sig_rho = r(p)

	qui corr `m1' `m2', cov
	local tau = _result(4)
	qui summ `m1'
	local xmin1 = _result(5)
	local xmax1 = _result(6)
	local err1 = _result(4) - `tau'
	qui summ `m2'
	local err2 = _result(4) - `tau'
	local xmin2 = _result(5)
	local xmax2 = _result(6)
	local xmin3 = floor(min(`xmin2', `xmin1'))
	local xmax3 = ceil(max(`xmax2', `xmax1'))





	#delim ;
		di in gr _n "Bland-Altman comparison of `m1' and `m2':";

	di in gr _n "Number of differences n = " in ye `format' `n';

	di in gr "Range `m1' : " in ye  `format' `xmin1'
		in gr " to " in ye `format' `xmax1';
	di in gr "Range `m2' : " in ye  `format' `xmin2'
		in gr " to " in ye `format' `xmax2';

	di in gr _n "95% Limits of agreement (Reference Range for difference): " ;
	di in gr "Lower Limit: " in ye  `format' `lrr'
		in gr " (CI " in ye `format' `lrr' - 1.96* `Se_limits'
		in gr " to " in ye `format' `lrr' + 1.96* `Se_limits'
		in gr ") " ;

	di in gr "Upper limit: " in ye `format' `urr'
	in gr " (CI " in ye `format' `urr' - 1.96* `Se_limits'
		in gr " to " in ye `format' `urr' + 1.96* `Se_limits'
		in gr ") " ;

      di in gr _n "Mean difference (bias `m1' - `m2'): " in ye `format' `xbar'
		in gr " (CI " in ye `format' `lcb'
		in gr " to " in ye `format' `ucb'
		in gr "  , p =" in ye `format' `sig1'
		in gr ") ";
	di in gr "Estimated sd on differences: " in ye `format' `sd'
		in gr " (CI " in ye `format' `Lower_sd'
		in gr " to " in ye `format' `Upper_sd'
		in gr ") ";

	di in gr _n "Correlation between difference and average: Pearsons r = " in ye `format' `r'
	   in gr ", p =" in ye `format' `sig2' ;
	di in gr "Correlation between difference and average: Spearmans rho = " in ye `format' `rho'
		in gr ", p =" in ye `format' `sig_rho' ;

	di in gr _n "Evaluation of one method:";

	di in gr "Estimated with-in subjects sd : " in ye `format' `sd'/sqrt(2)
		in gr " (CI " in ye `format' `Lower_sd'/sqrt(2)
		in gr " to " in ye `format' `Upper_sd'/sqrt(2)
		in gr ") ";
	di in gr "Estimated between subjects sd  : " in ye `format' sqrt(`var_between')
	in gr " (variance " in ye `format' `var_between'
		in gr ") ";
	di in gr "Estimated total sd : " in ye `format' sqrt(`var_total' )
	in gr " (variance " in ye `format' `var_total'
		in gr ") ";
	di in gr "Estimated ICC: " in ye `format' (`icc' );




	di in gr _n "Comparison of two methods (strong assumptions!!):";

	di in gr "Estimated between subject sd: " in ye `format' sqrt(`tau')
	in gr " (variance " in ye `format' `tau'
		in gr ") ";
	di in gr "Estimated additional error sd (`1'): " in ye `format' sqrt(`err1')
	in gr " (variance " in ye `format' `err1'
		in gr ") ";

	di in gr "Estimated additional error sd (`2'): " in ye `format' sqrt(`err2')
	in gr " (variance " in ye `format' `err2'
		in gr ") ";

	di in gr "Low or negative error variance indicates that model assumptions are violated.";

	#delim cr
	global S_1 `xbar'
	global S_2 `lrr'
	global S_3 `urr'

	if "`graph'" == "" {
		sort `diff' `av'
		tempvar f n
		qui by `diff' `av': gen `f' = _N if `diff' ~= . & `av' ~= .
		qui by `diff' `av': gen `n' = _n if `diff' ~= . & `av' ~= .
		qui graph twoway scatter `diff' `av' if `touse' & `n' == 1 [fw=`f'], symbol(`symbol') `xlabel' `ylabel' ///
 		yline(0,lco(black)) yline(`xbar' `lrr' `urr' ,lco(gray) ) yline(`lcb' `ucb',lpa(dash) lco(gray)) name(q1,replace) `options'
		if "`diag'" ~= "" {

			qui qnorm `diff' , name(q2,replace)
			qui qnorm `av' ,  name(q3,replace)
			qui twoway  ///
			(scatter `m2' `m1' if `touse' & `n' == 1 [fw=`f']  ,msy(O) symbol(`symbol') )  ///
			(function y=x, range(`xmin3'  `xmax3')), xtit(`m1' ) ytit(`m2' ) aspect(1) legend(off) name(q4,replace)
			graph combine q1 q2 q3 q4
			graph drop q1 q2 q3 q4

			}
		** graph twoway scatter `diff' `av' if `touse' & `n' == 1 [fw=`f'], symbol(`symbol') `xlabel' `ylabel'  yline(`yline') `saving' `options'
		}
end