*! 3.0.0 NJC 5 December 2023 
*! 2.0.0 NJC 11 December 2022
*! 1.0.0 NJC 9 March 2011
* qqplot 3.3.3  07mar2005
program qqplotg, sort
	version 8.2      
	
	local opts flip dvm dvp DIFFvsmean GENerate(str) TRANSform(str) 
	local opts `opts' by(str asis) a(str) MISSing lpolyopts(str asis)
	
	// support 2022 syntax group() diffvsmean 
	capture syntax varname(numeric) [if] [in], group(varname) [`opts' *] 

	if _rc capture syntax varname(numeric) [if] [in], over(varname) [`opts' *] 

	if _rc {
		syntax varlist(min=2 max=2 numeric) [if] [in] [, `opts' *]
	} 
		
	// a() expects a string to support input such as 1/3, but the result must be numeric 
	if "`a'" == "" local a = 0.5 
	else { 
		capture local A = `a'
		capture confirm number `A'
		if _rc { 
			di as err "a(`a') should define a number"
			exit 109 
		}
	}

	// supporting 2022 syntax 
	if "`group'" != "" local over `group'
	if "`diffvsmean'" != "" local dvm dvm
	
	if "`dvm'" != "" & "`dvp'" != "" {
		di as error "choose between difference vs mean and difference vs plotting position" 
		exit 198 
	}  

	tempvar touse 
	mark `touse' `if' `in'		/* but do not markout varlist */
	
	if "`by'" != "" { 
		gettoken byvar byopts : by, parse(,)
		Parsebyopts `byopts'
		local byby by(`byvar', `byopts')
			
		if "`missing'" == "" markout `touse' `byvar', strok 
	}
	else { 
		tempvar byvar 
		gen byte `byvar' = 1 
	}

	qui if "`over'" != "" { 
		tempname stub 
		separate `varlist' if `touse', by(`over') veryshortlabel ///
		gen(`stub')
		local varlist `r(varlist)' 

		local nvars : word count `varlist' 
		if `nvars' != 2 { 
			di as err "`nvars' categories in `over': should be 2"
			exit 498
		}
	} 
	
	if "`flip'" != "" { 
		tokenize "`varlist'" 
		local varlist "`2' `1'" 
	} 

	_get_gropts , graphopts(`options') getallowed(RLOPts plot addplot)
	local options `"`s(graphopts)'"'
	local rlopts `"`s(rlopts)'"'
	local plot `"`s(plot)'"'
	local addplot `"`s(addplot)'"'
	_check4gropts rlopts, opt(`rlopts')

	tokenize `varlist'
	tempvar VARY VARX CNT NEWX2 BYGROUP X  thisuse 
	quietly {
		gen `VARY' = `1' if `touse'
		gen `VARX' = `2' if `touse'
		
		if "`transform'" != "" { 
			Trans "`transform'" "`VARY'" "`touse'" "`1'" "`over'"
			Trans "`transform'" "`VARX'" "`touse'" "`2'" "`over'"
		}
		
		egen `BYGROUP' = group(`byvar') if `touse'
		su `BYGROUP', meanonly 
		gen byte `thisuse' = 0 
		gen `X' = . 
		
		forval g = 1/`r(max)' { 
			replace `thisuse' = `BYGROUP' == `g'
			bysort `thisuse' : gen long `CNT' = sum(`VARY'<.) if `thisuse'
			local cnty = `CNT'[_N]
			replace `CNT' = sum(`VARX'<.) if `thisuse'
			local cntx = `CNT'[_N]
			drop `CNT'
			
			if `cntx' == 0 | `cnty' == 0 error 2000 
			if `cnty' > `cntx' {
				QQp2 `VARY' `cnty' `cntx'
			}
			if `cnty' < `cntx' {
				QQp2 `VARX' `cntx' `cnty'  
			}

			QQp1 `VARY' `VARX' `NEWX2'
			replace `X' = `NEWX2' if `thisuse'
			drop `NEWX2'
		}

		_crcslbl `VARY' `1'
		_crcslbl `X' `2'
	}
	
	local yttl : var label `VARY'
	local xttl : var label `X'
	if `"`plot'`addplot'`by'"' == "" {
		local legend legend(nodraw)
	}

	if "`dvm'" != "" { 

	quietly { 
		tempvar diff mean 
		gen `diff' = `VARY' - `X' 
		gen `mean' = (`VARY' + `X') / 2
		local YTTL `"`yttl'"'  
		local yttl `"`yttl' - `xttl'"'
		local xttl `"(`YTTL' + `xttl') / 2"'
	} 
	

	version 8: graph twoway		///
	(scatter `diff' `mean',		///
		sort					///
		ytitle(`"`yttl'"')		///
		xtitle(`"`xttl'"')		///
		ms(Oh) mc(blue)         ///
		yla(, ang(h)) 			/// 
		`legend'				///
		`byby'                  ///
		`options'				///
	)							///
	(function 0,				///
		range(`mean')			///
		lstyle(refline)			///
		yvarlabel("Reference")	///
		`rlopts'				///
	)							///
	(lpoly `diff' `mean', `lpolyopts') ///
	|| `plot' || `addplot'		///
	// blank

	} 
	
	else if "`dvp'" != "" {
		
	quietly {
		tempvar N pp diff 
		
		gen `diff' = `VARY' - `X' 
		bysort `BYGROUP' (`X' `VARY') : gen `N' = sum(`X' < .)
		by `BYGROUP' : replace `N' = `N'[_N]
		by `BYGROUP' : gen `pp' = (_n - `a') / (`N' + 1 - 2 * (`a')) if `X' < . 
	
	    local YTTL `"`yttl'"'  
		local yttl `"`yttl' - `xttl'"'
		local xttl "Fraction of the data"
	} 
	
	
	version 8: graph twoway		///
	(scatter `diff' `pp',		///
		sort					///
		ytitle(`"`yttl'"')		///
		xtitle(`"`xttl'"')		///
		ms(Oh) mc(blue)         ///
		yla(, ang(h)) 			/// 
		xla(0 0.25 "0.25" 0.5 "0.5" 0.75 "0.75" 1)  /// 
		`legend'				///
		`byby'                  ///
		`options'				///
	)							///
	(function 0,				///
		range(`mean')			///
		lstyle(refline)			///
		yvarlabel("Reference")	///
		`rlopts'				///
	)							///
	(lpoly `diff' `pp', `lpolyopts') /// 
	|| `plot' || `addplot'		///
	// blank

			
	}

	else { 

	version 8: graph twoway		///
	(scatter `VARY' `X',	    ///
		sort					///
		ytitle(`"`yttl'"')		///
		xtitle(`"`xttl'"')		///
		ms(Oh) mc(blue)         ///
		yla(, ang(h)) 			/// 
		`legend'				///
		`byby'                  /// 
		`options'				///
	)							///
	(function y=x,				///
		range(`X')  			///
		lstyle(refline)			///
		yvarlabel("Reference")	///
		`rlopts'				///
	)							///
	|| `plot' || `addplot'		///
	// blank
	} 

	if "`generate'" != "" { 
		if "`dvm'" != "" { 
			confirm new variable `generate'd `generate'm  
			gen `generate'd = `diff'  
			label var `generate'd `"`yttl'"' 
			gen `generate'm = `mean' 
			label var `generate'm `"`xttl'"' 
		} 
		
		else if "`dvp'" != "" {
			confirm new variable `generate'd `generate'p 
			gen `generate'd = `diff'  
			label var `generate'd `"`yttl'"' 
			gen `generate'p = `pp' 
			label var `generate'p `"`xttl'"' 
		} 
	
		else {
			confirm new variable `generate'y `generate'x 
		
			gen `generate'y = `DQ' 
			label var `generate'y `"`yttl'"' 
			gen `generate'x = `NEWX2' 
			label var `generate'x `"`xttl'"'
		} 
	} 
end 

program QQp1 
	version 8.2
	tempvar YORDER
	quietly {
		sort `1'
		gen long `YORDER' = _n
		sort `2'
		gen `3' =`2'[`YORDER']
	}
end

program QQp2 /* varname old# new# */	
	version 8.2
	tempvar INT FRAC TEMP
	quietly {
		sort `1'
		gen long `INT' = (_n-0.5) * `2'/`3' + 0.5
		gen `FRAC' = 0.5 + (_n-0.5) * `2'/`3' - `INT'
		gen `TEMP' = (1-`FRAC') * `1'[`INT'] + `FRAC' * `1'[`INT'+1]
		replace `1' = `TEMP' 
	}
end

program Trans 
	args transform varname touse which over

	if strpos("`transform'", "@") {
		local recipe = subinstr("`transform'", "@", "`varname'", .)
	}
	else local recipe `transform'(`varname') 
		
	count if missing(`varname') & `touse'
	local bad = r(N)
	replace `varname' = `recipe' if `touse'
	count if missing(`varname') & `touse'
	local BAD = r(N)
			
	if `BAD' > `bad' {
		if "`over'" != "" { 
			su `over' if `varname' < ., meanonly 
			local where "for `over' == `r(min)'"
		}
		else local where "for `which'"
				
		local nbad = `BAD' - `bad'
		di as err "warning: `transform' leads to `nbad' missing " plural(`nbad', "value") " `where'"
	}
end 

program Parsebyopts 
	syntax , [ LEGend(str asis) note(str asis) * ]
	if `"`legend'"' == "" local legend off 
	c_local byopts legend(`legend') note(`"`note'"') `options'
end