*! version 7.1.7  04nov2009
* Modified from "sts.ado" and "logrank.ado"

* Originally written by David Fisher, 12th November 2009
* Modified (by David Fisher) as follows:

* May/June 2012, to incorporate aspects of meta-analysis
* ...and tidied up September 2012
* Various bug fixes October 2012

* Re-written to be a Peto meta-analysis command July 2013
* Includes all relevant code from "sts.ado", before the modified code from "logrank.ado"

* Note: This program deals with hazard ratios and meta-analysis
* and is therefore designed to handle TWO treatment groups only.
* For more general log-rank tests, use "sts test" as documented in the Stata manuals.

*! version 1.01  David Fisher  31jan2014

*! version 1.01.1 (beta)  David Fisher, November 2014
* Added glab label values
* Various refinements

* To do:
* Make -petometan- call -admetan- to draw tables/forestplot?? -- NO, TOO DIFFERENT
*   also display p-value and HRs plus O-E/V??

* Syntax:
* petometan trt_var [if] [in], [study(trial_id) by(subgroup) strata(other_strata)]
* where
* trt_var = treatment arm variable (coded such that trt=1, control=0)
* trial_id = trial identifier
* subgroup = optional trial-level subgroup identifier

* Data must be -stset-
* Data can be saved using MATsave option.

program define petometan, rclass sortpreserve
	
	version 8
	
	st_is 2 analysis
	local wt : char _dta[st_wt]
	if "`wt'"=="pweight" {
		disp as err `"Cannot specify pweights"'
		exit 198
	}
	
	syntax varname [if] [in] ///
		[, STUdy(varname) BY(varname) STRata(varlist) ///
		MATsave(name) noGRaph noHET noSUbgroup noOVerall OVSTAT(string) ///
		noTItle noSHow Level(cilevel) ]

	local arm "`varlist'"		// treatment arm
		
	st_show `show'
	tempvar touse
	st_smpl `touse' `"`if'"' "`in'"
	
	local w : char _dta[st_w]
	if `"`_dta[st_id]'"' != "" {
		local id `"id(`_dta[st_id]')"'
	}
	local t0 "_t0"
	local t1 "_t"
	local dead "_d"

	tempvar touse
	mark `touse' `if' `in' [`weight'`exp']
	markout `touse' `t1' `dead'
	markout `touse' `arm' `study' `by' `strata', strok
	
	* start of added section Nov 2014 (from logrank.ado)
	if `"`t0'"'!=`""' & `"`id'"'!=`""' {
		local id
	}
	if `"`t0'"'==`""' & `"`id'"'==`""' {
		tempvar t0
		qui gen byte `t0' = 0
	}
	else if `"`t0'"' != `""' { 
		markout `touse' `t0'
	}
	else if `"`id'"'!=`""' {
		markout `touse' `id'
		quietly {
			sort `touse' `id' `t1'
			local ty : type `t1'
			by `touse' `id': gen `ty' `t0' = cond(_n==1, 0, `t1'[_n-1])
		}
		capture assert `t1'>`t0'
		if _rc {
			di in red `"repeated records at same `t1' within `id'"'
			exit 498
		}
		* N.B. `id' is not needed any more
	}

	capture assert `t1'>0 if `touse'
	if _rc { 
		di in red `"survival time `t1' <= 0"'
		exit 498
	}
	capture assert `t0'>=0 if `touse'
	if _rc { 
		di in red `"entry time `t0' < 0"'
		exit 498
	}
	capture assert `t1'>`t0' if `touse'
	if _rc {
		di in red `"entry time `t0' >= survival time `t1'"'
		exit 498
	}
	capture assert `dead'==0 if `touse'
	if _rc==0 {
		di in red `"no test possible because there are no failures"'
		exit 2000
	}
	if `"`strata'"' != "" {
		tempvar isdead
		sort `strata'
		qui by `strata': gen long `isdead' = sum(`dead')
		qui by `strata': replace `isdead' = . if _n<_N
		qui count if `isdead' == 0 & `touse'
		local n_omit = r(N)
		if `n_omit' > 0 {
			if `n_omit' == 1 {
				local endng um
			}
			else {
				local endng a
			}
			local note `"Note: `n_omit' strat`endng' omitted because of no failures"'
		}
	}
	* end of added section
	
	* markout `touse' `t0'
	
	qui count if `touse'
	if !r(N) {
		disp as err "no observations"
		exit 2000
	}
	
	tempvar obs
	qui gen long `obs'=_n
	
	* Treatment arm variable should be coded such that "`arm'==1" denotes treatment and "`arm'==0" control
	summ `arm' if `touse', meanonly
	local armok = `r(min)'==0 & `r(max)'==1
	qui tab `arm' if `touse'
	local armok = `armok' * (`r(r)'==2)
	if !`armok' {
		di as err _n `"Treatment arm should be coded 0 = control, 1 = research"'
		exit 498
	}
	local g `arm'
	local glab : value label `g'
	
	if "`glab'"==`""' {
		tempname glab
		label define `glab' 0 "Control" 1 "Treatment"
	}
	
	/*
	tempvar g
	qui bysort `arm': gen long `g' = (_n==1)
	qui replace `g' = sum(`g')
	if `g'[_N]!=2 {
		di in red `"`by' variable is not binary"'
		exit 498
	}
	
	* Transfer labels, values or strings to new treatment arm variable
	tempname glab
	capture confirm string variable `arm'
	forvalues i=1/2 {
		if !_rc {
			summ `obs' if `g'==`i', meanonly
			local glabi = `arm'[`r(min)']
		}
		else {
			summ `arm' if `g'==`i', meanonly
			local glabi : label (`arm') `r(min)'		// N.B. will return `r(min)' if no label
			local glabi = cond("`glabi'"=="`r(min)'", cond(`i'==1, "Treatment", "Control"), "`glabi'")
		}
		label define `glab' `i' `"`glabi'"', add
	}
	*/
	* Store variable label
	local gvarlab : variable label `arm'
	if `"`gvarlab'"'==`""' local gvarlab `"`arm'"'
	
	* Create new study ID based on order of first occurrence
	local ns = 1					// default
	local nby = 0					// default
	if `"`study'"'!=`""' {
		tempvar s sobs
		qui bysort `touse' `study' (`obs') : gen long `sobs' = `obs'[1]
		qui bysort `touse' `by' `sobs' : gen long `s' = (_n==1) * `touse'
		qui replace `s' = sum(`s')
		local ns = `s'[_N]			// number of `study'*`by' groups -- should be equal to no. of studies!

		* Transfer labels, values or strings to new treatment arm variable
		tempname studylab
		capture confirm string variable `study'
		forvalues i=1/`ns' {
			if !_rc {
				summ `obs' if `s'==`i', meanonly
				local si = `study'[`r(min)']
			}
			else {
				summ `study' if `s'==`i', meanonly
				local si : label (`study') `r(min)'		// N.B. will return `r(min)' if no label
			}
			label define `studylab' `i' `"`si'"', add
		}
		* Store variable label
		local svarlab : variable label `study'
		if `"`svarlab'"'==`""' local svarlab `"`study'"'
	
		* Create new subgroup ID (BY) based on "natural ordering"
		if `"`by'"'!=`""' {
			
			* Check that `by' is trial-level
			qui tab `s' if `touse'
			if `r(r)' != `ns' {			// N.B. `ns' is already stratified by `by'
				disp as err _n "Data is not suitable for meta-analysis" _c
				disp as err " as subgroup variable (in option 'by') is not constant within trials."
				exit 198
			}
			else {
				tempvar byobs bygroup
				qui bysort `touse' `by' : gen int `bygroup' = (_n==1) * `touse'
				qui replace `bygroup' = sum(`bygroup')
				local nby = `bygroup'[_N]				// number of subgroups
				
				* Transfer labels, values or strings to new BY variable
				sort `obs'
				tempname bylab
				capture confirm string variable `by'
				forvalues i=1/`nby' {
					if !_rc {
						summ `obs' if `bygroup'== `i', meanonly
						local bylabi = `by'[`r(min)']
					}
					else {
						summ `by' if `bygroup'== `i', meanonly
						local bylabi : label (`by') `r(min)'		// N.B. will return `r(min)' if no label
					}
					label define `bylab' `i' `"`bylabi'"', add
				}
				
				* Store variable label
				local byvarlab : variable label `by'
				if `"`byvarlab'"'==`""' local byvarlab `"`by'"'
				local by `"`bygroup'"'
			}
		}
	}
	else if `"`by'"'!=`""' {
		disp as err "Cannot specify 'by' without 'study'"
		exit 198
	}

	
	* Begin manipulating data
	preserve 
	
	if `"`weight'"' != `""' { 
		tempvar w 
		qui gen double `w' `exp' if `touse'
		local wv `"`w'"'
		local wntype "double"	// "gen double `n`i''" if weights
	}
	else {
		local w 1
		local wntype "long"		// "gen long `n`i''" if no weights (since in that case must be whole numbers)
	}
	tempvar op n d
	quietly {
		keep if `touse'
		qui count
		return scalar N = r(N)
		keep `s' `g' `wv' `t0' `t1' `by' `strata' `dead'
		
		* Denominators
		* (need to calculate these before limiting to unique times only)
		* Only need to know denominators per study, per subgroup, and overall
		* Strata are irrelevant as main calculations don't use denoms, & strata-specific stats are not presented
		if trim(`"`study' `by'"') != `""' {
			local bystr `"by `by' `s':"'
		}
		* forvalues i=1/2 {
		forvalues i=0/1 {
			sort `by' `s' `g' `t1'
			tempvar NN`i'
			`bystr' gen long `NN`i''=sum(cond(`g'==`i',1,0))
			sort `by' `s' `NN`i''
			`bystr' replace `NN`i''=`NN`i''[_N]
		}
		
		* Now re-define "bystr" for main calculations
		* This time "by" is irrelevant since it must be trial-level
		* but "strata" ARE relevant
		if trim(`"`study' `by' `strata'"') != `""' {
			local bystr `"by `s' `strata':"'
		}		
		local N = _N
		expand 2
		gen byte `op' = 3/*add*/ in 1/`N'
		replace `t1' = `t0' in 1/`N'
		drop `t0'
		local ++N
		replace `op' = cond(`dead'==0,2/*cens*/,1/*death*/) in `N'/l

		sort `s' `strata' `t1' `op' `g'
		`bystr' gen `wntype' `n' = sum(cond(`op'==3,`w',-`w'))
		by `s' `strata' `t1': gen `wntype' `d' = sum(`w'*(`op'==1))

		* Numbers at risk, and observed number of events (failures)
		* forvalues i=1/2 {
		forvalues i=0/1 {
			tempvar ni`i' di`i'
			`bystr' gen `wntype' `ni`i'' = sum(cond(`g'==`i', cond(`op'==3,`w',-`w'), 0))
			by `s' `strata' `t1': gen `wntype' `di`i'' = sum(cond(`g'==`i', `w'*(`op'==1), 0))
			* N.B. `w' is not needed any more
		}
		by `s' `strata' `t1': keep if _n==_N		// keep unique times only

		* Shift `n' up one place so it lines up
		tempvar newn
		`bystr' gen `wntype' `newn' = `n'[_n-1]
		drop `n' 
		rename `newn' `n'
		
		* Shift each of the `ni's up one place so they line up
		* forvalues i=1/2 {
		forvalues i=0/1 {
			`bystr' gen `wntype' `newn' = `ni`i''[_n-1] if _n>1
			drop `ni`i''
			rename `newn' `ni`i''	
		}
		* drop if `d'==0			// keep failure times only - DON'T DO THIS, IN CASE SOME STUDIES HAVE NO FAILURES
		capture drop `strata'		// don't need strata vars anymore (and there may be many of them)

		* Calculate E (expected number of events/failures)
		* forvalues i=1/2 {
		forvalues i=0/1 {
			tempvar ei`i'
			gen double `ei`i'' = `ni`i''*`d'/`n'
		}
		* Calculate V (hypergeometric variance)
		tempvar V
		assert float(`ni0' + `ni1') == float(`n')		// arithmetic check
		gen double `V' = `ni0'*`ni1'*`d'*(`n'-`d')/(`n'*`n'*(`n'-1))

		* Calculate O - E
		tempvar OE OEsq
		gen double `OE'=`di1'-`ei1'						// use treatment arm
		assert float(`OE') == float(`ei0'-`di0')		// arithmetic check
		
		* At this point we have one obs per unique failure time per arm per trial.
		* Now "collapse" to one obs per study (or just one obs), plus (sub)totals
		if `"`study'"'!=`""' {
			sort `s'
			local bys `"by `s':"'
		}
		* foreach x of varlist `OE' `V' `di1' `di2' {
		foreach x of varlist `OE' `V' `di0' `di1' {
			`bys' replace `x' = sum(`x')
		}
		`bys' keep if (_n == _N)
		tempvar obs
		gen long `obs' = 1
		
		* Generate new obs for subgroup/overall statistics
		if `"`study'"'!=`""' {
		
			assert `ns' == _N
			local newn = `ns' + (`"`overall'"'==`""') + (`"`subgroup'"'==`""')*`nby'	// create rows to hold subgroup & overall totals
			set obs `newn'
			tempvar newobs
			gen byte `newobs' = (_n > `ns')		// flag new observations
			
			tempvar Q
			gen double `Q'=.
			label var `Q' "Q"
			
			* Subgroup totals
			if `"`by'"'!=`""' & `"`subgroup'"'==`""' {
				forvalues i=1/`nby' {
					local ii = `ns' + `i'
					replace `by' = `i' in `ii'
				}
				sort `by' `s'
				* foreach x of varlist `V' `OE' `di1' `di2' `NN1' `NN2' `obs' {
				foreach x of varlist `V' `OE' `di0' `di1' `NN0' `NN1' `obs' {
					tempvar bysum`x'
					by `by' : gen double `bysum`x'' = sum(`x')
					by `by' : replace `bysum`x'' = `bysum`x''[_N]
				}
				tempvar qpart
				by `by' : gen double `qpart' = sum(`V'*(((`OE'/`V')-(`bysum`OE''/`bysum`V''))^2))
				by `by' : replace `qpart' = `qpart'[_N]
				replace `Q' = `qpart' if `newobs' & !missing(`by')
				drop `qpart'
			}

			* Overall totals
			if `"`overall'"'==`""' {
				* foreach x of varlist `V' `OE' `di1' `di2' `NN1' `NN2' `obs' {
				foreach x of varlist `V' `OE' `di0' `di1' `NN0' `NN1' `obs' {
					tempvar sum`x'
					gen double `sum`x'' = sum(`x')
					replace `sum`x'' = `sum`x''[_N]
					replace `x' = `sum`x'' if `newobs' & missing(`Q')		// insert totals in new rows
				}
			
				tempvar Vsq qtot 
				gen double `Vsq' = `V'^2
				gen double `qtot'=sum(`V'*(((`OE'/`V')-(`sum`OE''/`sum`V''))^2))
				replace `qtot' = `qtot'[_N]
				replace `Q' = `qtot' if `newobs' & missing(`Q')
				
				local OEtot = `sum`OE''
				local Vtot = `sum`V''
			}
		}
		else {						// if `"`study'"'==`""'
			local OEtot = `OE'
			local Vtot = `V'
			local sum`di0' = `di0'
			local sum`di1' = `di1'
			* local sum`di2' = `di2'
		}

		if `"`overall'"'==`""' {
			local lnHR = `OEtot'/`Vtot'
			local selnHR = 1/sqrt(`Vtot')
			local chi2 = (`OEtot'^2)/`Vtot'			

			* return scalar o = `sum`di1''+`sum`di2''
			return scalar o = `sum`di0''+`sum`di1''
			return scalar OE = `OEtot'
			return scalar V = `Vtot'
			return scalar lnHR =`lnHR'
			return scalar selnHR = `selnHR'
			return scalar chi2 = `chi2'			
		}
		
		if `"`study'"'!=`""' {
			if `"`overall'"'==`""' {
				* drop `sum`di1'' `sum`di2'' `sum`NN1'' `sum`NN2'' `sum`OE'' `sum`V''
				drop `sum`di0'' `sum`di1'' `sum`NN0'' `sum`NN1'' `sum`OE'' `sum`V''
			}
			if `"`by'"'!=`""' & `"`subgroup'"'==`""' {
				* foreach x of varlist `V' `OE' `di1' `di2' `NN1' `NN2' `obs' {
				foreach x of varlist `V' `OE' `di0' `di1' `NN0' `NN1' `obs' {
					replace `x' = `bysum`x'' if `newobs' & !missing(`by')	// insert totals in new rows
				}
				* drop `bysum`di1'' `bysum`di2'' `bysum`NN1'' `bysum`NN2'' `bysum`OE'' `bysum`V''
				drop `bysum`di0'' `bysum`di1'' `bysum`NN0'' `bysum`NN1'' `bysum`OE'' `bysum`V''
			}
		}
	}			// end "quietly"
	
	* Print to screen
	tempvar counts0 counts1 /*counts2*/
	qui gen str `counts0' = string(`di0') + "/" + string(`NN0')		// don't use tempvar here so can send to forestplot
	qui gen str `counts1' = string(`di1') + "/" + string(`NN1')		// don't use tempvar here so can send to forestplot
	* qui gen str `counts2' = string(`di2') + "/" + string(`NN2')
	label var `counts0' "`: label `glab' 0'"						// REVISIT: replace with glab value labels?
	label var `counts1' "`: label `glab' 1'"						// REVISIT: replace with glab value labels?
	* label var `counts2' "`: label `glab' 2'"
	label var `OE' "o-E(o)"
	label var `V' "Var(o)"

	local bystr
	if `"`by'"'!=`""' {
		local bystr `"by(`by')"'
		label values `by' `bylab'
		label var `by' `"`byvarlab'"'
	}
	if `"`study'"'!=`""' {
		label values `s' `studylab'
		label var `s' `"`svarlab'"'
	}
	else {
		tempvar s
		qui gen byte `s'=1
		label var `s' "Study"
	}
	* tabdisp `s', cell(`counts1' `counts2' `OE' `V') format(%04.2f) totals concise `bystr'
	tabdisp `s', cell(`counts0' `counts1' `OE' `V') format(%04.2f) totals concise `bystr'
	
	* Display effect sizes & tests
	if `"`overall'"'==`""' {
		local CI_lo=exp(`lnHR'-(invnorm((100+`level')/200)*`selnHR'))
		local CI_hi=exp(`lnHR'+(invnorm((100+`level')/200)*`selnHR'))
		local Qtot = `Q'[_N]
		local df = `ns'-1
		return scalar Q = `Qtot'
		return scalar k = `ns'
		local Qpval = chi2tail(`df', `Qtot')
		disp _n
		if `"`by'"'!=`""' disp as text "Overall"
		disp as text "Pooled hazard ratio = " as res %5.3f `=exp(`lnHR')'
		disp as text "Two-sided `level'% confidence limit = (" as res %5.3f `CI_lo' as text ", " as res %5.3f `CI_hi' as text ")"
	}
		
	if `"`study'"'!=`""' {
		if `"`overall'"'==`""' {
			disp as text "Q for heterogeneity = " as res %5.3f `Qtot' as text " on " as res `df' as text " d.f., p = " as res %5.3f `Qpval'
		}
	
		* Subgroups
		if `"`by'"'!=`""' & `"`subgroup'"'==`""' {
			local Qsum = 0
			sort `newobs' `by' `s'
			forvalues i=1/`nby' {
				local lnHR`i' = `OE'[`=`ns'+`i'']/`V'[`=`ns'+`i'']
				local selnHR`i' = 1/sqrt(`V'[`=`ns'+`i''])
				local Q`i' = `Q'[`=`ns'+`i'']
				local Qsum = `Qsum' + `Q`i''
				local df`i' = `obs'[`=`ns'+`i''] - 1
				local Qpval`i' = chi2tail(`df`i'', `Q`i'')
				local CI_lo=exp(`lnHR`i''-(invnorm((100+`level')/200)*`selnHR`i''))
				local CI_hi=exp(`lnHR`i''+(invnorm((100+`level')/200)*`selnHR`i''))
				disp as text _n `"`: label `bylab' `i''"'
				disp as text "Pooled hazard ratio = " as res %5.3f `=exp(`lnHR`i'')'
				disp as text "Two-sided `level'% confidence limit = (" as res %5.3f `CI_lo' as text ", " as res %5.3f `CI_hi' as text ")"
				disp as text "Q for heterogeneity = " as res %5.3f `Q`i'' as text " on " as res `df`i'' as text " d.f., p = " as res %5.3f `Qpval`i''
			}
			
			* Between-subgroup heterogeneity
			if `"`overall'"'==`""' {
				local Qb = `Qtot' - `Qsum'
				local dfb = `nby' - 1
				local Qbpval = chi2tail(`dfb', `Qb')
				disp as text _n `"Between-subgroup heterogeneity = "' as res %5.3f `Qb' as text " on " as res `dfb' as text " d.f., p = " as res %5.3f `Qbpval'
			}
		}
	}
	
	* Matrix output
	if `"`matsave'"'!=`""' {
		* qui mkmat `s' `by' `NN1' `NN2' `di1' `di2' `ei1' `OE' `V' if !missing(`s'), matrix(`matsave')
		qui mkmat `s' `by' `NN0' `NN1' `di0' `di1' `ei1' `OE' `V' if !missing(`s'), matrix(`matsave')
		if `"`by'"'!=`""' local byname "by"
		matrix colnames `matsave' = s `byname' N1 N2 d1 d2 e OE V
	}

	* Output data to forestplot program
	if `"`graph'"'==`""' {

		quietly {
			* drop `NN1' `NN2' `di1' `di2' `ei1' `ei2'		
			drop `NN0' `NN1' `di0' `di1' `ei0' `ei1'		
			
			rename `OE' OE
			rename `V' V
			format OE %5.2f
			format V %5.2f

			rename `counts0' counts0
			rename `counts1' counts1
			* rename `counts2' counts2
			rename `s' _STUDY
			if `"`by'"'!=`""' {
				rename `by' _BY
				local _by "_BY"
			}
					
			gen _ES = OE/V
			gen _seES = 1/sqrt(V)
			gen _WT = V/`Vtot'
			if `"`: value label _STUDY'"'!=`""'	decode _STUDY, gen(_LABELS)
			else gen _LABELS = string(_STUDY)
			
			gen _USE=.
			replace _USE=1 if !missing(_STUDY)
			if `"`by'"'!=`""' {
				replace _USE=3 if missing(_STUDY) & !missing(_BY)
				replace _USE=5 if missing(_BY)
			}
			else replace _USE=5 if missing(_STUDY) 
			
			tempvar use5
			qui gen `use5' = (_USE==5)
			
			if `"`by'"'!=`""' {					// subgroup titles
				tempvar expand
				bysort _BY : gen byte `expand' = 1 + 2*(_n==1)*(!`use5')
				expand `expand'
				gsort _BY -`expand' _USE _STUDY
				by _BY : replace _USE=0 if `expand'>1 & _n==2		/* row for labels */
				by _BY : replace _USE=4 if `expand'>1 & _n==3		/* row for blank line */
				drop `expand'
				
				if "`subgroup'"=="" & "`het'"=="" {
					tempvar expand
					bysort _BY : gen byte `expand' = 1 + (_n==_N)*(!`use5')
					expand `expand'
					gsort _BY -`expand' _USE _STUDY
					by _BY : replace _USE=3.5 if `expand'>1 & _n==2 		/* extra row for het */
				}
			}
			
			* Overall heterogeneity - extra row
			if `"`overall'"'==`""' & "`het'"=="" {
				local newN = `=_N+1'
				set obs `newN'
				replace _USE=5.5 in `newN'
			}
			
			* Blank out effect sizes etc. in new rows
			foreach x of varlist _LABELS _ES _seES _WT counts0 counts1 OE V {
				capture confirm numeric variable `x'
				if !_rc replace `x' = . if !inlist(_USE, 1, 3, 5)
				else replace `x' = "" if !inlist(_USE, 1, 3, 5)
			}

			* Add between-group heterogeneity info if appropriate
			if `"`by'"'!=`""' & `"`overall'"'==`""' & `"`het'"'==`""' {
				local newN = _N+1
				set obs `newN'
				replace _USE = 4.5 in `newN'
				replace `use5' = 0 in `newN'
				replace _LABELS = "Heterogeneity between groups: p = " + string(`Qbpval', "%5.3f") in `newN'
			}
			
			replace _LABELS = "Overall" if _USE==5
			
			if "`ovstat'"=="q" { 
				local ovlabel "(Q = " + string(`Qtot', "%5.2f") + " on `df' df, p = " + string(`Qpval', "%5.3f") + ")"
			}
			else {
				local Isq = max(0, (`Qtot' - `df') / `Qtot')
				local ovlabel "(I-squared = " + string(100*`Isq', "%5.1f")+ "%, p = " + string(`Qpval', "%5.3f") + ")"
			}
			replace _LABELS = "`ovlabel'" if _USE==5.5
		
			* Subgroup ("by") labels
			forvalues i=1/`nby' {
				replace _LABELS = "Subtotal" if _USE==3 & _BY==`i'
				
				if "`ovstat'"=="q" {
					local ovlabel "(Q = " + string(`Q`i'', "%5.2f") + " on `dfi' df, p = " + string(`Qpval`i'', "%5.3f") + ")"
				}
				else {
					local Isq`i' = max(0, (`Q`i'' - `df`i'') / `Q`i'')
					local ovlabel "(I-squared = " + string(100*`Isq`i'', "%5.1f")+ "%, p = " + string(`Qpval`i'', "%5.3f") + ")"
				}
				replace _LABELS = "`ovlabel'" if _USE==3.5 & _BY==`i'
					
				local bytext : label `bylab' `i'
				replace _LABELS = "`bytext'" if _USE==0 & _BY==`i'
			}
			
			sort `_by' _USE _STUDY
			replace _USE = 0 if _USE == -1
			replace _USE = 4 if inlist(_USE, -0.5, -1.5, 2.5, 4.5)
			replace _USE = 3 if _USE == 3.5
			replace _USE = 5 if _USE == 5.5
			
			gen _LCI = _ES - invnorm(0.975)*_seES
			gen _UCI = _ES + invnorm(0.975)*_seES

		}	// end quietly
		
		order _USE `_by' _STUDY _LABELS _ES _seES _LCI _UCI _WT

		forestplot, nopreserve by(`_by') labels(_LABELS) eform ///
				nowt nostats lcols(counts0 counts1) rcols(OE V)
	}

end