* REMR version 1.2 - 15 October 2021
* Authors: Luis Furuya-Kanamori (l.furuya@uq.edu.au), Chang Xu, Suhail AR Doi

	*version 1.1 - ref(#) not part of dataset, added warning message ""Graph not rendered - reference value should match one of the dose levels in the dataset"
	*version 1.2 - keep all data after if/in selected


program define remr, rclass
version 14

syntax varlist(min=3 max=4 numeric) [if] [in] ///
, ID(varname) CATegory(varname) [, or rr rd noGraph bplot Table eform] ///
[rcs(numlist max=1)] [Knots(numlist)] [REFerence (numlist max=1)] [CENter (numlist max=1)] [STARTdose (numlist min=2 max=2)]

tokenize `varlist'

preserve
marksample touse, novarlist 
quietly keep if `touse'

quietly{

*Error messages
	*OR and RR
	if "`or'"!="" & "`rr'"!=""{
	noisily di as error "or and rr cannot both be specified"
	exit 198
	}	
	
	*OR and RD
	if "`or'"!="" & "`rd'"!=""{
	noisily di as error "or and rd cannot both be specified"
	exit 198
	}	
	
	*RR and RD
	if "`rr'"!="" & "`rd'"!=""{
	noisily di as error "rr and rd cannot both be specified"
	exit 198
	}
		
	*OR, RR, and RD with data entry as "[1]es [2]se [3]dose"
	if ("`or'"!="" | "`rr'"!="" | "`rd'"!="") & "`4'"=="" {
	noisily di as text "Note: effect size is specified by `1'"
	}	
	
	*Missing ID or category
	cap assert !missing(`id')
	if _rc {
	noisily di as err "`id' cannot contain missing values"
	exit 198
	}
	cap assert !missing(`category')
	if _rc {
	noisily di as err "`category' cannot contain missing values"
	exit 198
	}
	
	*Reference and Center
	if "`reference'"!="" & "`center'"==""{
	noisily di as error "reference option can only be used when center has been specified"
	exit 198
	}
			
	*Category check
	sort `id' `category' 
	bysort `id': gen __category_check = _n-1 
	cap assert `category' == __category_check
	if _rc & ("`if'"=="" & "`in'"==""){
	noisily di as err "category wrongly specified"
	exit 198
	}
	
	if _rc & ("`if'"!="" | "`in'"!=""){
	noisily di as err "verify category is correctly specified"
	}
	

	*Graph
	if "`graph'" == "nograph" & "`bplot'" !="" {
	di as error "nograph and bplot cannot both be specified"
	exit 198	
	}
	if "`bplot'" !="" & "`reference'" !="" {
	di as error "bplot cannot be displayed when a reference has been specified"
	exit 198	
	}	
		
	
*Data entry 
	*[1]cases [2]non-cases [3]total [4]dose
		if "`3'"!="" & "`4'"!=""{
		
		gen __drop =.
		replace __drop = 1 if `1'==0 & `2'==0 
		replace __drop = 1 if `1'==. | `2'==. | `4'==.
		sort `id' `category' 
		bysort `id': replace __drop = 1 if __drop[1] == 1 	
		drop if __drop==1
				
		gen __data_entry_type = "Data input format cases non_cases total dose, id() cat()"

		sort `id' `category' 
		bysort `id': gen __study_id=_n	
		
			*OR (by default)
			if "`or'"!="" | ("`or'"=="" & "`rr'"=="" & "`rd'"==""){
			gen __continuity = 1 if (`1' ==0 | `2' ==0)
				replace `1' = `1'+0.5 if __continuity==1 
				replace `2' = `2'+0.5 if __continuity==1 
			gen __effect_size = "OR"
			gen __or=.
			replace __or=1 if __study_id==1
			bysort `id': replace __or=(`1'*`2'[1]) / (`2'*`1'[1]) if __study_id!=1
			gen __es=ln(__or)
			gen __se=.
			replace __se=0 if __study_id==1
			bysort `id': replace __se=sqrt((1/`1')+(1/`2')+(1/`1'[1])+(1/`2'[1])) if __study_id!=1
			gen __dose=`4'
			}
		
			*RR
			if "`rr'"!=""{
			gen __continuity = 1 if (`1' ==0 | `2' ==0)
				replace `1' = `1'+0.5 if __continuity==1 
				replace `2' = `2'+0.5 if __continuity==1 
			gen __effect_size = "RR"
			gen __rr = .
			replace __rr=1 if __study_id==1
			bysort `id': replace __rr= (`1'/(`1'+`2')) / (`1'[1]/(`1'[1]+`2'[1])) if __study_id!=1
			gen __es=ln(__rr)
			gen __se=.
			replace __se=0 if __study_id==1
			bysort `id': replace __se=sqrt((1/`1')+(1/`1'[1])-(1/(`1'+`2'))-(1/(`1'[1]+`2'[1]))) if __study_id!=1
			gen __dose=`4'
			}
			
			*RD (continuity only for SE)
			if "`rd'"!=""{
			gen __continuity = 1 if (`1' ==0 | `2' ==0)
			gen __effect_size = "RD"
			gen __rd = .
			replace __rd=0 if __study_id==1
			bysort `id': replace __rd= (`1'/(`1'+`2')) - (`1'[1]/(`1'[1]+`2'[1])) if __study_id!=1
			gen __es=__rd
				replace `1' = `1'+0.5 if __continuity==1 
				replace `2' = `2'+0.5 if __continuity==1 
			gen __se=.
			replace __se=0 if __study_id==1
			bysort `id': replace __se=sqrt( ((`1'*`2')/(`1'+`2')^3) + ((`1'[1]*`2'[1])/(`1'[1]+`2'[1])^3) ) if __study_id!=1
			gen __dose=`4'
			}	
		}
		
	*[1]es [2]se [3]dose
		if "`3'"!="" & "`4'"==""{
		
		gen __drop =.
		replace __drop = 1 if `1'==. | `3'==.
		replace __drop = 1 if `2'==. & `category' !=0
		sort `id' `category' 
		bysort `id': replace __drop = 1 if __drop[1] == 1 	
		drop if __drop==1
				
		gen __data_entry_type = "Data input format es se dose, id() cat()"
		
		gen __es = `1'
		gen __se = `2'
		gen __dose=`3'
		
		sort `id' `category'
		bysort `id': gen __study_id=_n	
		}
		

*Center + start dose
	if "`center'"=="" & "`startdose'"!=""{
	noisily di as err "startdose cannot be specified without center"
	exit 198
	}
		
	if "`center'"!="" & "`startdose'"==""{
	local cent `center'
	gen __cent = `center'
	bysort `id': gen __cent_start = __dose[1]
	bysort `id': gen __dose_cent = (__dose-__cent_start)+`cent'
	replace __dose = __dose_cent
	}

	if "`center'"!="" & "`startdose'"!=""{
	local cent `center'
	gen __cent = `center'
	
	tokenize `startdose'
		local startdose1 `1'
		local startdose2 `2'
	
		cap assert `startdose1'<=`cent' | `cent'>=`startdose2'
		if _rc {
		noisily di as err "center(#) has to be within startdose(#1 #2) range"
		exit 198
		}
	
	bysort `id': gen __cent_start = __dose[1]
	bysort `id': gen __dose_cent = (__dose-__cent_start)+`cent'
	replace __dose = __dose_cent
	
	drop if __cent_start < `startdose1' 
	drop if __cent_start > `startdose2'
	
		cap assert _N != 0
		if _rc {
		noisily di as err "studies' lowest dose outside startdose(#1 #2) range"
		exit 198	
		}
	}
	
	
*Chang-Doi method
	gen __wt = 1/(__se^2)
	bysort `id': egen __maxwt = max(__wt)
	replace __wt =. if `category'==0
	replace __wt = __maxwt if __wt==. & `category'==0
		
	*Restricted cubic spline option
	if "`rcs'"=="" & "`knots'"==""{
		gen __model = "Linear"
		gen __doses1 = __dose 
		gen __n_knots = "None"
		gen __knot_val1 = "None"
	}
	
	if "`rcs'"!="" & "`knots'"==""{
	local rcs `rcs'
	mkspline __doses = __dose, cubic nk(`rcs')	dis
		gen __model = "Non-linear"
		gen __n_knots = r(N_knots)
		matrix knots = r(knots)
		svmat knots, names(__knot_val)
	} 
	
 	if  "`rcs'"=="" & "`knots'"!="" {
	local nc 0
	tokenize "`knots'"
	mkspline __doses = __dose, cubic knots(`knots')	dis
		gen __model = "Non-linear"
		gen __n_knots = r(N_knots)
		matrix knots = r(knots)
		svmat knots, names(__knot_val)
	}
	
	if "`rcs'"!="" & "`knots'"!="" {
	local rcs `rcs'
	local nc 0
	tokenize "`knots'"
	mkspline __doses = __dose, cubic knots(`knots')	nk(`rcs') dis
		gen __model = "Non-linear"
		gen __n_knots = r(N_knots)
		matrix knots = r(knots)
		svmat knots, names(__knot_val)
	}
		

	*Calculations + options eform & reference
	sort __doses1
	list __doses* if __doses1 != __doses1[_n-1]
		regress __es __doses* [aweight=__wt], vce(cluster `id') eform (exp beta)
		matrix model = r(table)
		levelsof __dose, local(levels)
			if "`eform'" !="" & "`reference'"!=""{
				local ref `reference'
				gen __ref = `reference'
				xblc __doses* , covname (__dose) at(`r(levels)') eform gen(__d __xb __lci __uci) ref(`ref')
			}
			if "`eform'" =="" & "`reference'"!=""{
				local ref `reference'
				gen __ref = `reference'
				xblc __doses* , covname (__dose) at(`r(levels)') gen(__d __xb __lci __uci) ref(`ref')
			}
			if "`eform'" !="" & "`reference'"==""{
				xblc __doses* , covname (__dose) at(`r(levels)') eform gen(__d __xb __lci __uci)
			}
			if "`eform'" =="" & "`reference'"==""{
				xblc __doses* , covname (__dose) at(`r(levels)') gen(__d __xb __lci __uci)
			}
			
		
*Display results
noisily di ""
noisily di in smcl as txt "{hline 19}{c TT}{hline 65}"
	
	noisily di in smcl as text "Data entry type" as text "{col 20}{c |}" ///
	"{col 25}" (__data_entry_type[1])
	
	if __data_entry_type=="Data input format cases non_cases total dose, id() cat()"{
	noisily di in smcl as text "Effect size" as text "{col 20}{c |}" ///
	"{col 25}" (__effect_size[1])
	}
	
	noisily di in smcl as text "Model" as text "{col 20}{c |}" ///
	"{col 25}" (__model[1])
	
	if "`center'"!=""{
	noisily di in smcl as text "Center value" as text "{col 20}{c |}" ///
	"{col 25}" (__cent[1])
	}
	
	if "`reference'"!=""{
	noisily di in smcl as text "Reference value" as text "{col 20}{c |}" ///
	"{col 25}" (__ref[1])
	}
	
	noisily di in smcl as text "Num knots" as text "{col 20}{c |}" ///
	"{col 25}" (__n_knots[1])

	if "`rcs'"=="" & "`knots'"==""{
	noisily di in smcl as text "Knot values" as text "{col 20}{c |}" ///
	"{col 25}" (__knot_val1[1])
	}

	if "`rcs'"!="" | "`knots'"!=""{
	sort __knot_val1
		if __n_knots[1]==2{
		noisily di in smcl as text "Knot values" as text "{col 20}{c |}" ///
		"{col 25}" int(__knot_val1[1]) ///
		"{col 28}" int(__knot_val2[1])
		}
		if __n_knots[1]==3{
		noisily di in smcl as text "Knot values" as text "{col 20}{c |}" ///
		"{col 25}" int(__knot_val1[1]) ///
		"{col 28}" int(__knot_val2[1]) ///
		"{col 31}" int(__knot_val3[1])
		}
		if __n_knots[1]==4{
		noisily di in smcl as text "Knot values" as text "{col 20}{c |}" ///
		"{col 25}" int(__knot_val1[1]) ///
		"{col 28}" int(__knot_val2[1]) ///
		"{col 31}" int(__knot_val3[1]) ///
		"{col 34}" int(__knot_val4[1])
		}
		if __n_knots[1]==5{
		noisily di in smcl as text "Knot values" as text "{col 20}{c |}" ///
		"{col 25}" int(__knot_val1[1]) ///
		"{col 28}" int(__knot_val2[1]) ///
		"{col 31}" int(__knot_val3[1]) ///
		"{col 34}" int(__knot_val4[1]) ///
		"{col 37}" int(__knot_val5[1]) 
		}
		if __n_knots[1]==6{
		noisily di in smcl as text "Knot values" as text "{col 20}{c |}" ///
		"{col 25}" int(__knot_val1[1]) ///
		"{col 28}" int(__knot_val2[1]) ///
		"{col 31}" int(__knot_val3[1]) ///
		"{col 34}" int(__knot_val4[1]) ///
		"{col 37}" int(__knot_val5[1]) ///
		"{col 40}" int(__knot_val6[1]) 
		}
		if __n_knots[1]==7{
		noisily di in smcl as text "Knot values" as text "{col 20}{c |}" ///
		"{col 25}" int(__knot_val1[1]) ///
		"{col 28}" int(__knot_val2[1]) ///
		"{col 31}" int(__knot_val3[1]) ///
		"{col 34}" int(__knot_val4[1]) ///
		"{col 37}" int(__knot_val5[1]) ///
		"{col 40}" int(__knot_val6[1]) ///
		"{col 43}" int(__knot_val7[1]) 
		}
		if __n_knots[1]==8{
		noisily di in smcl as text "Knot values" as text "{col 20}{c |}" ///
		"{col 25}" int(__knot_val1[1]) ///
		"{col 28}" int(__knot_val2[1]) ///
		"{col 31}" int(__knot_val3[1]) ///
		"{col 34}" int(__knot_val4[1]) ///
		"{col 37}" int(__knot_val5[1]) ///
		"{col 40}" int(__knot_val6[1]) ///
		"{col 43}" int(__knot_val7[1]) ///
		"{col 46}" int(__knot_val8[1]) 
		}
		if __n_knots[1]==9{
		noisily di in smcl as text "Knot values" as text "{col 20}{c |}" ///
		"{col 25}" int(__knot_val1[1]) ///
		"{col 28}" int(__knot_val2[1]) ///
		"{col 31}" int(__knot_val3[1]) ///
		"{col 34}" int(__knot_val4[1]) ///
		"{col 37}" int(__knot_val5[1]) ///
		"{col 40}" int(__knot_val6[1]) ///
		"{col 43}" int(__knot_val7[1]) ///
		"{col 46}" int(__knot_val8[1]) ///
		"{col 49}" int(__knot_val9[1]) 
		}
	}
	
	gen __num_obs = e(N)
	noisily di in smcl as text "Num obs" as text "{col 20}{c |}" ///
	"{col 25}" (__num_obs[1])
	
	gen __num_studies = .
	bysort `id': replace __num_studies = 1 if __num_studies[1]==.
	replace __num_studies = sum(__num_studies)
	gsort -__num_studies
	noisily di in smcl as text "Num studies" as text "{col 20}{c |}" ///
	"{col 25}" (__num_studies[1])
	
	gen __r2 = e(r2)
	noisily di in smcl as text "R-squared" as text "{col 20}{c |}" ///
	"{col 25}" round(__r2[1],0.001)
	
	gen __root_mse = e(rmse)
	noisily di in smcl as text "Root MSE" as text "{col 20}{c |}" ///
	"{col 25}" round(__root_mse[1],0.001)
	
noisily di in smcl as txt "{hline 19}{c BT}{hline 65}"

	*Output table
		if "`table'" != ""{
			if "`eform'" !="" & "`reference'"!=""{
			noisily regress __es __doses* [aweight=__wt], vce(cluster `id') eform (exp beta)
			matrix model = r(table)
			levelsof __dose, local(levels)
			local ref `reference'
			if "`rd'"!=""{
			noisily xblc __doses* , covname (__dose) at(`r(levels)') eform ref(`ref') format(%7.3f)
			}
				else{
				noisily xblc __doses* , covname (__dose) at(`r(levels)') eform ref(`ref') format(%7.2f)
				}				
			}
			if "`eform'" =="" & "`reference'"!=""{
			noisily regress __es __doses* [aweight=__wt], vce(cluster `id') 
			matrix model = r(table)
			levelsof __dose, local(levels)
			local ref `reference'
			if "`rd'"!=""{
			noisily xblc __doses* , covname (__dose) at(`r(levels)') ref(`ref') format(%7.3f)
			}
				else{
				noisily xblc __doses* , covname (__dose) at(`r(levels)') ref(`ref') format(%7.2f)
				}
			}			
			if "`eform'" !="" & "`reference'"==""{
			noisily regress __es __doses* [aweight=__wt], vce(cluster `id') eform (exp beta)
			matrix model = r(table)
			levelsof __dose, local(levels)
			if "`rd'"!=""{
			noisily xblc __doses* , covname (__dose) at(`r(levels)') eform format(%7.3f)
			}
				else{
				noisily xblc __doses* , covname (__dose) at(`r(levels)') eform format(%7.2f)
				}
			}
			if "`eform'" =="" & "`reference'"==""{
			noisily regress __es __doses* [aweight=__wt], vce(cluster `id')
			matrix model = r(table)
			levelsof __dose, local(levels)
			if "`rd'"!=""{
			noisily xblc __doses* , covname (__dose) at(`r(levels)') format(%7.3f)
			}
				else{
				noisily xblc __doses* , covname (__dose) at(`r(levels)') format(%7.2f)
				}				
			}
		}

		
*Matrices (+ Regression model[above])
sort __d

	cap mkmat __d __xb __lci __uci, matrix(output) nomissing
	cap return matrix output = output

	cap return matrix model = model
	
	cap mkmat __doses*, matrix(doses) nomissing
	cap return matrix doses = doses

		
*Graph
	if "`graph'" != "nograph" & "`bplot'" == "" & (__lci==. & __uci==.){
		noisily di as err "graph not rendered - reference value should match one of the dose levels in the dataset"
		}

	if "`graph'" != "nograph" & "`bplot'" == "" & (__lci!=. & __uci!=.){
		sort __d
		tw (line __xb __d, lcolor(black) lpattern(l) lwidth(thick)) /// 
		(line __lci __d, lcolor(black) lpattern(dash) lwidth(thick)) ///
		(line __uci __d, lcolor(black) lpattern(dash) lwidth(thick)), /// 
		legend(off) sch(s1mono) ///	
		xtitle("dose") ///                    
		plotregion(style(none)) ///
		ylabel(#7, angle(horiz) format(%3.2fc)) ///
		xlabel(#7)
		}
		
	if "`bplot'" != "" {
		bysort `id': egen __dosemin = min(__dose) 
		replace __wt=. if __dose == __dosemin	
		
			if "`eform'" =="eform" {
			replace __es = exp(__es)
			}
		
		sort __d
		tw (scatter __es __dose [aw=__wt] if __se!=., msymbol(oh)) ///
		(line __xb __d, lcolor(black) lpattern(l) lwidth(thick)) /// 
		(line __lci __d, lcolor(black) lpattern(dash) lwidth(thick)) ///
		(line __uci __d, lcolor(black) lpattern(dash) lwidth(thick)), /// 
		legend(off) sch(s1mono) ///	
		xtitle("dose") ///                    
		plotregion(style(none)) ///
		ylabel(#7, angle(horiz) format(%3.2fc)) ///
		xlabel(#7)
		}
}		
*

*Restore data and exit
restore 
end
exit