*! version 7.0.2 6Apr2001
* Written by Peter Cummings
/* csmatch.ado */
program define csmatch, rclass
	version 7.0
	syntax varlist(min=2 max=2 numeric) [if] [in], /*
	*/ [level(real .95)] PAIRMember(string) /*
	*/ Group(varname) [PERSONVar(varlist)] [PAIRVar(varlist)]
	marksample touse
	tokenize `varlist'
	tempvar outcome
	qui gen byte `outcome' = `1' if `touse'
	tempvar expos
	qui gen byte `expos' = `2' if `touse'
	tempvar drivern
	qui gen byte `drivern'= `pairmember' if `touse'
	tempvar pass
	qui gen byte `pass' = `drivern'==0 if `touse'
	markout `touse' `pass' `personvar' `pairvar'
	preserve
	qui keep if `touse'

	/* check that outcome is coded 0/1 */
	capture assert `outcome'==0 | `outcome'==1 | `outcome'==.
		if _rc~=0 {
			di as error "Outcome " as result "`1'" /*
			*/ as error " not coded 0/1"
			exit 450
		}

	/* check exposure is coded 0/1 */
	capture assert `expos'==0 | `expos'==1 | `expos'==.
		if _rc~=0 {
			di as error "Exposure " as result "`2'" /*
			*/ as error " not coded 0/1"
			exit 450
		}

	/* check that groups are never greater than 2 */
	sort `group'
	tempvar grpct
	qui by `group': gen `grpct' = _N
	capture assert `grpct' <3
		if _rc~=0 {
			di as error "Some groups are greater than 2"
			exit 459
		}

	/* convert 2 record format to 1 record format */
	sort `group' `pass'
	tempvar miss1 miss2
	qui by `group': gen byte `miss1' = /*
		*/ `touse'[2] == 0 | `expos'[2]==. | `outcome'[2]==.
	qui by `group': gen byte `miss2' = /*
		*/ `touse'[1] == 0 | `expos'[1]==. | `outcome'[1]==.
	qui replace `miss1' = . if `miss1' ==1
	qui replace `miss2' = . if `miss2' ==1
	markout `touse' `miss1' `miss2'
	tempvar dout dexp
	qui by `group': gen byte `dout' = `outcome'[1] if `touse'
	qui by `group': gen byte `dexp' = `expos'[1] if `touse'
	tempvar pout pexp
	qui by `group': gen byte `pout' = `outcome'[2] if `touse'
	qui by `group': gen byte `pexp' = `expos'[2] if `touse'

	/* deal with occupant by variables. Create a temporary
	variable for each confounder in occupant list and create driver
	and passenger versions */

	local oct: word count `personvar'
	tokenize `personvar'
	local i 1
	/* generate the new occupant variables and put their names into
	local macro stro */
	local stro = ""		/* collect all the new variable names */
	local strod= ""		/* collect the new driver variable names */
	local strop= ""		/* collect the new passenger var names */
	while `i' <= `oct' {
		tempvar dobyv`i' pobyv`i'
		qui by `group': gen int `dobyv`i'' = ``i''[1] if `touse'
		qui by `group': gen int `pobyv`i'' = ``i''[2] if `touse'
		local stro = "`stro' `dobyv`i'' `pobyv`i'' "
		local strod = "`strod' `dobyv`i''"
		local strop = "`strop' `pobyv`i''"
		local i = `i' + 1
		}

	/* eliminate all records in which the two occupants are not the 
	same regarding the personvar variables */
	local i 1
	tempvar nokeep
	qui gen byte `nokeep' = 1
	while `i' <= `oct' {
		qui replace `nokeep' = . if `dobyv`i'' ~= `pobyv`i''
		local i = `i' + 1
		}

	/* change 'touse' to include only the driver record */
	tempvar dronly
	qui gen byte `dronly' = 1 if `pass' ==0
	/* use only records that have at least one occupant with 
	the outcome */
	tempvar onedead
	qui gen byte `onedead' = 1 if `pout'==1 | `dout'==1
	/* use only records that are discordant on exposure */
	tempvar oneexp
	qui gen byte `oneexp' = 1 if `pexp' ~= `dexp'
	markout `touse' `dronly' `onedead' `nokeep' `oneexp'
	qui keep if `touse'
	tempvar prct
	qui count
	qui gen long `prct' = r(N)	/* count pairs available */
	
	/* deal with pairvar, the vehicle variables */
	local vct: word count `pairvar'

	/* GENERATE OVERALL ESTIMATE */
	/* generate the exposed to unexposed ratio */
	/* this "crude" estimate uses all pairs that are in touse, 
	regardless of whether some strata are missing an estimate */
	tempvar a0 b0 c0
	tempvar aaa0 bbb0 ccc0
	tempvar rr0    /* exposed to unexposed rr */
	tempvar vlrr0  /* variance of ln rr */

	/* a b c d */
	qui gen byte `aaa0' = `dout'==1 & `pout'==1 & /*
		*/ (`dexp'==1 & `pexp'==0 /*
		*/ | `pexp'==1 & `dexp'==0)
	qui gen long `a0' = sum(`aaa0')
	qui replace `a0'  = `a0'[_N]
	qui gen byte `bbb0' = `dexp'==1 & `dout'==0 & `pexp'==0 & `pout'==1 /*
			*/ |  `pexp'==1 & `pout'==0 & `dexp'==0 & `dout'==1
	qui gen long `b0' = sum(`bbb0')
	qui replace `b0'  = `b0'[_N]
	qui gen byte `ccc0' = `dexp'==1 & `dout'==1 & `pexp'==0 & `pout'==0 /*
			*/ |  `pexp'==1 & `pout'==1 & `dexp'==0 & `dout'==0
	qui gen long `c0' = sum(`ccc0')
	qui replace `c0'  = `c0'[_N]

	qui gen double `rr0' = (`a0'+`c0')/(`a0'+`b0')
	qui gen double `vlrr0' = (`b0'+`c0')/[(`a0'+`b0')*(`a0'+`c0')]
	tempvar rru0
	tempvar rrl0
	qui gen double `rru0' = exp(ln(`rr0') /*
		*/ + invnorm(`level'*.5 + 0.5)*sqrt(`vlrr0'))
	qui gen double `rrl0' = exp(ln(`rr0') /*
		*/ - invnorm(`level'*.5 + 0.5)*sqrt(`vlrr0'))
	tempvar serr0 	/* standard error of rr */
	qui gen double `serr0' = exp(sqrt(`vlrr0'))
	
	return scalar rr = `rr0'
	return scalar vlrr = `vlrr0'
	return scalar prct = `prct'

	/* generate results by strata according
	to the personvar and pairvar variables */

	/* start a big if/else loop here. Within the first part
		calculate and display the stratum specific estimates.
		The last part will just display crude estimates. */

	if "`stro'`pairvar'" ~="" {	/* large if/else starts here */
		tempvar strata
		qui egen int `strata' = group(`stro' `pairvar')
		qui summarize `strata'
		tempvar strmax
		qui gen `strmax' = r(max)

	/* generate counts by stratum */
	/* generate the driver belted to passenger unbelted ratio */
	tempvar last
	tempvar a b c
	tempvar aaa bbb ccc
	tempvar rr    /* exposed to unexposed */
	tempvar vlrr  /* variance of ln rr */

	/* a b c d */
	sort `strata'
	tempvar last
	qui gen byte `last'=0
	qui by `strata': replace `last'=1 if _n== _N
	qui by `strata': gen byte `aaa' = `dout'==1 & `pout'==1 & /*
		*/ (`dexp'==1 & `pexp'==0 /*
		*/ | `pexp'==1 & `dexp'==0)
	qui by `strata': gen long `a' = sum(`aaa')
	qui by `strata': replace `a'  = `a'[_N]
	qui by `strata': gen byte `bbb' = `dexp'==1 & `dout'==0 & /*
			*/ `pexp'==0 & `pout'==1 /*
			*/ |  `pexp'==1 & `pout'==0 & `dexp'==0 & `dout'==1
	qui by `strata': gen long `b' = sum(`bbb')
	qui by `strata': replace `b'  = `b'[_N]
	qui by `strata': gen byte `ccc' = `dexp'==1 & `dout'==1 & /*
			*/ `pexp'==0 & `pout'==0 /*
			*/ |  `pexp'==1 & `pout'==1 & `dexp'==0 & `dout'==0
	qui by `strata': gen long `c' = sum(`ccc')
	qui by `strata': replace `c'  = `c'[_N]

	local i 1
	while `i'<=`strmax' {
		qui count if `strata' ==`i'
		tempvar prct`i'
		qui gen long prct`i' = r(N)
		return scalar prct`i' = r(N)
		local i = `i' +1
		}
	
	tempvar rru
	tempvar rrl
	qui keep if `last'==1
	qui gen double `rr' = (`a'+`c')/(`a'+`b')
	qui gen double `vlrr' = (`b'+`c')/[(`a'+`b')*(`a'+`c')]
	qui gen double `rru' = exp(ln(`rr') /*
		*/ + invnorm(`level'*.5 + 0.5)*sqrt(`vlrr'))
	qui gen double `rrl' = exp(ln(`rr') /*
		*/ - invnorm(`level'*.5 + 0.5)*sqrt(`vlrr'))  
		
	/* tests of homogeneity */
	tempvar chi df pchi
	qui gen double `chi' = sum(((ln(`rr')-ln(`rr0'))^2)/`vlrr')
	qui replace `chi' = `chi'[_N]
	qui count if `rr'~=.
	qui gen int `df' = r(N) -1
	qui gen double `pchi' = chi2tail(`df',`chi')

	/* collect stratum specific scalars */
	sort `strata'
	local i 1
	while `i'<=`strmax' {
		return scalar rr`i' = `rr'[`i']
		return scalar vlrr`i' = `vlrr'[`i']
		local i = `i' +1
		}

	/* create macros for value labels for strata display */

	local names `personvar' `personvar' `pairvar'
	tokenize `names'
 	capture local lbl1 : value label `1'
	capture local lbl2 : value label `2'
	capture local lbl3 : value label `3'
	capture local lbl4 : value label `4'
	capture local lbl5 : value label `5'
	capture local lbl6 : value label `6'
		
	/* confine output length to 72 characters */
	di
	local c1 = 73 - (8 + length("Available pairs = "))
	local c2 = 73 - (8 + length("Pairs used = "))
	di as text "Matched-pair cohort method" _col(`c1') /*
		*/ "Available pairs = " as result %8.0g `prct'
	di as text "Stratified by: `personvar' `pairvar'"
	di as text "{hline 8}{c TT}{hline 63}"

	if `oct' ~= 0 & `vct' ~= 0 {
		local colp = 11 + 9*`oct'	/* passenger column */
		local colv = 11 + 2*9*`oct'	/* vehicle column */
		di as text _col(9)"{c |}" _col(11)"Driver" /*
			*/ _col(`colp')"Pass"  _col(`colv')"Vehicle"
	}
	else if `oct'==0 & `vct' ~= 0 {
		di as text _col(9)"{c |}" _col(11)"Vehicle"
	}
	else {
		local colp = 11 + 9*`oct'	/* passenger column */
		di as text _col(9)"{c |}" _col(11)"Driver" /*
			*/ _col(`colp')"Pass"
	}

	/* variable list will be 2*oct + 1*vct */
	tokenize `personvar' `personvar' `pairvar'
	di as text "stratum {c |}" /*
		*/ _col(11)"`1'" _col(20)"`2'" /*
		*/ _col(29)"`3'" _col(38)"`4'" /*
		*/ _col(47)"`5'" _col(56)"`6'" _col(65)"count"
	di as text "{hline 8}{c +}{hline 63}"

	local totct: word count `personvar' `personvar' `pairvar'
	tokenize `strod' `strop' `pairvar'
	sort `strata'
	if `totct' == 6 {
		label values `1' `lbl1'		/* label values */
		label values `2' `lbl2'
		label values `3' `lbl3'
		label values `4' `lbl4'
		label values `5' `lbl5'
		label values `6' `lbl6'
		local i 1
		while `i' <=`strmax' {
		if "`lbl1'" ~= "" {
			local disp1 = `1'[`i']
			local dsplbl1 : label `labl1' `disp1'
			}
		else {
			capture local dsplbl1 = `1'[`i']
		}
		if "`lbl2'" ~= "" {
			local disp2 = `2'[`i']
			local dsplbl2 : label `labl2' `disp2'
			}
		else {
			capture local dsplbl2 = `2'[`i']
		}
		if "`lbl3'" ~= "" {
			local disp3 = `3'[`i']
			local dsplbl3 : label `labl3' `disp3'
			}
		else {
			capture local dsplbl3 = `3'[`i']
		}
		if "`lbl4'" ~= "" {
			local disp4 = `4'[`i']
			local dsplbl4 : label `labl4' `disp4'
			}
		else {
			capture local dsplbl4 = `4'[`i']
		}
		if "`lbl5'" ~= "" {
			local disp5 = `5'[`i']
			local dsplbl5 : label `labl5' `disp5'
			}
		else {
			capture local dsplbl5 = `5'[`i']
		}
		if "`lbl6'" ~= "" {
			local disp6 = `6'[`i']
			local dsplbl6 : label `labl6' `disp6'
			}
		else {
			capture local dsplbl6 = `6'[`i']
		}
		
		di as result `stratum'[`i'] /*
			*/ as text _col(9)"{c |}" as result /*
			*/ _col(11)"`dsplbl1'" /*
			*/ _col(20)"`dsplbl2'" /*
			*/ _col(29)"`dsplbl3'" /*
			*/ _col(38)"`dsplbl4'" /*
			*/ _col(47)"`dsplbl5'" /*
			*/ _col(56)"`dsplbl6'" /*
			*/ _col(65) prct`i'[`i']
		local i = `i'+1
		}
	}
	else if `totct' == 5 {
		label values `1' `lbl1'		/* label values */
		label values `2' `lbl2'
		label values `3' `lbl3'
		label values `4' `lbl4'
		label values `5' `lbl5'
		local i 1
		while `i' <=`strmax' {
		if "`lbl1'" ~= "" {
			local disp1 = `1'[`i']
			local dsplbl1 : label `lbl1' `disp1'
			}
		else {
			capture local dsplbl1 = `1'[`i']
		}
		if "`lbl2'" ~= "" {
			local disp2 = `2'[`i']
			local dsplbl2 : label `lbl2' `disp2'
			}
		else {
			capture local dsplbl2 = `2'[`i']
		}
		if "`lbl3'" ~= "" {
			local disp3 = `3'[`i']
			local dsplbl3 : label `lbl3' `disp3'
			}
		else {
			capture local dsplbl3 = `3'[`i']
		}
		if "`lbl4'" ~= "" {
			local disp4 = `4'[`i']
			local dsplbl4 : label `lbl4' `disp4'
			}
		else {
			capture local dsplbl4 = `4'[`i']
		}
		if "`lbl5'" ~= "" {
			local disp5 = `5'[`i']
			local dsplbl5 : label `lbl5' `disp5'
			}
		else {
			capture local dsplbl5 = `5'[`i']
		}
		
		di as result `stratum'[`i'] /*
			*/ as text _col(9)"{c |}" as result /*
			*/ _col(11)"`dsplbl1'" /*
			*/ _col(20)"`dsplbl2'" /*
			*/ _col(29)"`dsplbl3'" /*
			*/ _col(38)"`dsplbl4'" /*
			*/ _col(47)"`dsplbl5'" /*
			*/ _col(65) prct`i'[`i']
		local i = `i'+1
		}
	}
	else if `totct' == 4 {
		label values `1' `lbl1'		/* label values */
		label values `2' `lbl2'
		label values `3' `lbl3'
		label values `4' `lbl4'
		local i 1
		while `i' <=`strmax' {
		if "`lbl1'" ~= "" {
			local disp1 = `1'[`i']
			local dsplbl1 : label `lbl1' `disp1'
			}
		else {
			capture local dsplbl1 = `1'[`i']
		}
		if "`lbl2'" ~= "" {
			local disp2 = `2'[`i']
			local dsplbl2 : label `lbl2' `disp2'
			}
		else {
			capture local dsplbl2 = `2'[`i']
		}
		if "`lbl3'" ~= "" {
			local disp3 = `3'[`i']
			local dsplbl3 : label `lbl3' `disp3'
			}
		else {
			capture local dsplbl3 = `3'[`i']
		}
		if "`lbl4'" ~= "" {
			local disp4 = `4'[`i']
			local dsplbl4 : label `lbl4' `disp4'
			}
		else {
			capture local dsplbl4 = `4'[`i']
		}
		
		di as result `stratum'[`i'] /*
			*/ as text _col(9)"{c |}" as result /*
			*/ _col(11)"`dsplbl1'" /*
			*/ _col(20)"`dsplbl2'" /*
			*/ _col(29)"`dsplbl3'" /*
			*/ _col(38)"`dsplbl4'" /*
			*/ _col(65) prct`i'[`i']
		local i = `i'+1
		}
	}
	else if `totct' == 3 {
		label values `1' `lbl1'		/* label values */
		label values `2' `lbl2'
		label values `3' `lbl3'
		local i 1
		while `i' <=`strmax' {
		if "`lbl1'" ~= "" {
			local disp1 = `1'[`i']
			local dsplbl1 : label `lbl1' `disp1'
			}
		else {
			capture local dsplbl1 = `1'[`i']
		}
		if "`lbl2'" ~= "" {
			local disp2 = `2'[`i']
			local dsplbl2 : label `lbl2' `disp2'
			}
		else {
			capture local dsplbl2 = `2'[`i']
		}
		if "`lbl3'" ~= "" {
			local disp3 = `3'[`i']
			local dsplbl3 : label `lbl3' `disp3'
			}
		else {
			capture local dsplbl3 = `3'[`i']
		}
		
		di as result `stratum'[`i'] /*
			*/ as text _col(9)"{c |}" as result /*
			*/ _col(11)"`dsplbl1'" /*
			*/ _col(20)"`dsplbl2'" /*
			*/ _col(29)"`dsplbl3'" /*
			*/ _col(65) prct`i'[`i']
		local i = `i'+1
		}
	}
	else if `totct' == 2 {
		label values `1' `lbl1'		/* label values */
		label values `2' `lbl2'
		local i 1
		while `i' <=`strmax' {
		if "`lbl1'" ~= "" {
			local disp1 = `1'[`i']
			local dsplbl1 : label `lbl1' `disp1'
			}
		else {
			capture local dsplbl1 = `1'[`i']
		}
		if "`lbl2'" ~= "" {
			local disp2 = `2'[`i']
			local dsplbl2 : label `lbl2' `disp2'
			}
		else {
			capture local dsplbl2 = `2'[`i']
		}
		
		di as result `stratum'[`i'] /*
			*/ as text _col(9)"{c |}" as result /*
			*/ _col(11)"`dsplbl1'" /*
			*/ _col(20)"`dsplbl2'" /*
			*/ _col(65) prct`i'[`i']
		local i = `i'+1
		}
	}
	else {
		label values `1' `lbl1'		/* label values */
		local i 1
		while `i' <=`strmax' {
		if "`lbl1'" ~= "" {
			local disp1 = `1'[`i']
			local dsplbl1 : label `lbl1' `disp1'
			}
		else {
			capture local dsplbl1 = `1'[`i']
		}		

		di as result `stratum'[`i'] /*
			*/ as text _col(9)"{c |}" as result /*
			*/ _col(11)"`dsplbl1'" /*
			*/ _col(65) prct`i'[`i']
		local i = `i'+1
		}
	}


	di as text "{hline 8}{c BT}{hline 63}"
	local levpc = `level'*100
	di as text "Relative Risk estimates" /*
		*/ _col(39)"Confidence interval level = " /*
		*/ %4.1f in ye `levpc' " %"
	di as text "{hline 8}{c TT}{hline 63}"
	di as text "stratum {c |}" /*
		*/ _col(11) "E+O+/E-O+  E+O-/E-O+  E+O+/E-O- " /*
		*/ "      RR    [Conf.   Interval]"
	di as text "{hline 8}{c +}{hline 63}"
	local i 1
	while `i' <=`strmax' {
		di as result `stratum'[`i'] /*
			*/ as text _col(9)"{c |}" as result /*
			*/ _col(13) %7.0f `a'[`i'] /*
			*/ _col(24) %7.0f `b'[`i'] /*
			*/ _col(35) %7.0f `c'[`i'] /*
			*/ _col(48) %6.4f `rr'[`i'] /*
			*/ _col(55) "[" %6.4f `rrl'[`i'] "," /*
			*/ _col(66) %6.4f `rru'[`i'] "]"
		local i = `i'+1
		}
	di as text "{hline 8}{c +}{hline 63}"
	di as text "combined{c |}" as result /*
			*/ _col(48) %6.4f `rr0' /*
			*/ _col(55) "[" %6.4f `rrl0' "," /*
			*/ _col(66) %6.4f `rru0' "]"
	di as text "{hline 8}{c BT}{hline 63}"
	di as text "Homogeneity test of RR estimates " /*
		*/ as result _col(63) "p = " %5.4f `pchi'
	di as text "{hline 72}

	/* now comes the second part of the big if/else loop, 
		showing results if there is only one stratum */
	}
	else 	{
		di
		local c1 = 73 - (8 + length("Available pairs = "))
		di as text "Matched-pair cohort method" _col(`c1') /*
			*/ "Available pairs = " as result %8.0g `prct'
		local levpc = `level'*100
		di as text "Relative Risk estimates" /*
		*/ _col(39)"Confidence interval level = " as text /*
		*/ %4.1f `levpc' " %"
		di as text "{hline 72}"
		di as text "E+O+/E-O+    E+O-/E-O+    E+O+/E-O- " /*
			*/ _col(40)"RR"  _col(55)"[Conf.  Interval]"
		di as text "{hline 72}"
		di as result _col(3) %7.0f `a0' /*
			*/ _col(16) %7.0f `b0' /*
			*/ _col(27) %7.0f `c0' /*
			*/ _col(40) %6.4f `rr0' /*
			*/ _col(55) as text "[" as result %6.4f `rrl0' /*
			*/ as text ", " /*
			*/ as result _col(64) %6.4f `rru0' "]" 
		di as text "{hline 72}"
	}

end