*! version 1.0

capture program drop dunnett
program define dunnett, rclass

	syntax varlist(min=2 max=2) [if] [in] , CONtrol(real) /*
	*/ [ci(integer 95) Side(integer 2) NOlabel SEParator(integer 5) MIssing]
       	gettoken resp fact: varlist
       
 	/* Mark sample for computations */
      /* If missing is specified, treat missing as category */
	if "`missing'" == ""{
	   	marksample touse 
         	qui count if `touse'
      }

	else{
		marksample touse, novarlist
		qui count if `touse'
	}

	/* Make sure we do not have zero observations in sample */
	if r(N)==0 {
           	error 2000
      }

	/* Check assumptions */
	
	/* Check that separator() is positive */
	if `separator' <0 {
		di in red "separator() incorrectly specified"
		exit 198
	}

	/* Check that ci() option is only 95 or 99 */
      if `ci' != 95 & `ci' != 99{
          	di in red "ci() must be 95 or 99"
           	exit 198
      }
	
	/* Check that side() option is only 1 or 2 */
      if `side' != 1 & `side' != 2{
      	di in red "side() must be 1 or 2"
      	exit 198
      }

	/* Check that number of treatments are less than 22 for 2 sided interval */
      qui levelsof `fact' if `touse', local(levels) missing
      local k: word count `levels'

	/* k is total number of treatments including control */
      if `k' >21 & `side'==2{
      	di in red "two sided intervals must have 21 or less treatments including control"
      	exit 198
      }
	
	/* Check that number of treatments are less than 11 for 1 sided interval */
      if `k' >10 & `side'==1{
      	di in red "one sided intervals must have 10 or less treatments including control"
      	exit 198
      }

	/* Make sure control() is in factor_var */
      local infact = 0
      foreach factor of local levels{
		if `factor' == `control'{
			local infact = 1
		}
      }
      if `infact' == 0{
   		di in red "control = `control' is not in `fact'"
   		exit 198
      }

	/* Check that sample size is the same for all groups */
	qui count if `fact' == `control' & `touse'
	local n = r(N)
	foreach i of local levels{
		qui count if `fact' == `i'  & `touse'
 		if r(N) != `n'{
			di in red "sample sizes of each group must be equal"
			exit 198
		}
	}

	/* Count observations in sample */
	qui count if `touse'
  	local N = r(N)
  	
	/* Call oneway and save mse and df */
  	qui oneway `resp' `fact' if `touse',missing 
 	local mse = r(rss)/r(df_r)
  	local df = r(df_r)
  
	/* Make sure df is at least 5 */
  	if `df' < 5{
    		di in red "degrees of freedom must be at least five"
   		exit 198
  	}

	/* Call program d_table for table value and store in d */
 	d_table `side' `ci' `df' `k' 
  	local d = `r(d)'
	
	/* Calculate Dunnett criterion value */
  	local criterion = `d' * sqrt(2*`mse'/(`n'))
  
	/* Create local list of treatments excluding control */
	local list : list levels - control

	/* If option nolabel is not specified */
	if "`nolabel'" == ""{
		/* extract value labels if labels exist */
		local factorlab: value label `fact'
		if "`factorlab'" != ""{
			local lab`control': label `factorlab' `control'
			if `"`lab`control''"' == ""{
				local lab`control' = `control'
			}
			foreach group of local list{
				if "`group'" == "."{
					local labmiss = .
				}
				else{
					local lab`group': label `factorlab' `group'
					if `"`lab`group''"' == ""{
						local lab`group' = `group'
					}
				}
			}
		}
		/* if factor_var does not have labels */
		if "`factorlab'" == ""{
			local lab`control' = `control'
			foreach group of local list{
				if "`group'" == "."{
					local labmiss = .
				}
				else{
					local lab`group' = `group'
				}
			}

		}
	}
	/* if option nolabel is specified */
	if "`nolabel'" != ""{
		local lab`control' = `control'
		foreach group of local list{
			if "`group'" == "."{
				local labmiss = .
			}
			else{
				local lab`group' = `group'
			}
		}
	}

	/* calculate control mean */
  	qui sum `resp' if `fact'==`control' & `touse', meanonly
  	local control_mean = r(mean)


	/* Calculations and table output for 2 sided intervals */
	if `side' == 2{  
		/* Calculate differences and confidence interval */
  		foreach group of local list{
			if "`group'" == "."{
				local g = "miss"
			}
			else{
				local g = `group'
			}
			qui sum `resp' if `fact'==`group' &  `touse', meanonly
			local mean`g' = r(mean)
			local diff`g' = `mean`g'' - `control_mean'
			local abs_diff`g' = abs(`diff`g'')
			local lb`g' = `diff`g'' - `criterion'
			local ub`g' = `diff`g'' + `criterion'
			if `abs_diff`g'' > `criterion'{
				local conclusion`g'  "Yes"
			}
			else{
				local conclusion`g'  "No"
			}
  		}
		/* Create table output */
	
		/* Header */
		#delimit ;
		di in smcl in gr "{hline 9}{c TT}{hline 74}";
		di in smcl in gr _col(9) " {c |}" _col(74) "Different";
		di in smcl in gr %8s abbrev("`fact'",8) " {c |}"
		_col(17) "Mean"
		_col(29) "Diff"
		_col(36) "[ 2-Sided  `ci'% SCI ]"
		_col(60) "abs(Diff)"
		_col(72) "from Control?";
		#delimit cr
		di in smcl in gr "{hline 9}{c +}{hline 74}"

		/* Output for control */
		#delimit ;
		di in smcl in gr %8s abbrev("`lab`control'' ",8) " {c |}" in ye
		_col(12) %9.0g `control_mean'
		_col(30) "---"
		_col(42) "---        ---"
		_col(66) "---"
		_col(77) "---";
		#delimit cr

		/* Output for all treatments excluding control */
		local i = 0
		foreach group of local list{
			local i = `i' + 1
			if mod(`i',`separator')==0{
				di in smcl in gr  "{hline 9}{c +}{hline 74}"
			}
			if "`group'" == "."{
				local g = "miss"
				local shortlab = .
			}
			else{
				local g = `group'
				local shortlab = abbrev("`lab`g''",8)
			}
			#delimit ;
			di in smcl in gr  %8s "`shortlab'" " {c |}" in ye
			_col(12) %9.0g `mean`g''
			_col(24) %9.0g `diff`g''
			_col(36) %9.0g `lb`g'' "  " %9.0g `ub`g''  
			_col(60) %9.0g `abs_diff`g''
			_col(77) %3s "`conclusion`g''" ;
			#delimit cr
		}
  		
		/* Bottom line and notes */
 		di in smcl in gr  "{hline 9}{c BT}{hline 74}"
		local shortfact = abbrev("`fact'",8)
 		di in smcl in gr _col(6) "Diff = mean(`shortfact')-mean(control)"
 		di in smcl in gr  _col(6) "Different from mean(control) if abs(Diff) > " %8.0g `criterion'
	}


	/* Calculations and table output for 1 sided interval */
	if `side' == 1{
	
		/* Calculate differences and lower and upper bounds */
		foreach group of local list{
			if "`group'" == "."{
				local g = "miss"
			}
			else{
				local g = `group'
			}
			qui sum `resp' if `fact'==`group' &  `touse', meanonly
			local mean`g' = r(mean)
			local diff`g' = `mean`g'' - `control_mean'
			local lb`g' = `diff`g'' - `criterion'
			local ub`g' = `diff`g'' + `criterion'
			if `diff`g'' > `criterion'{
				local greater`g'  "Yes"
			}
			else{
				local greater`g'  "No"
			}
			if `diff`g'' < -`criterion'{
				local less`g' "Yes"
			}
			else{
				local less`g' "No"
			}
		}

		/* Header */
		#delimit ;
		di in smcl in gr "{hline 9}{c TT}{hline 82}";
		di in smcl in gr _col(9) " {c |}"
		_col(38) "`ci'% SCI"
		_col(52) "`ci'% SCI"
		_col(67) "Greater"
		_col(84) "Less";
		di in smcl in gr %8s abbrev("`fact'",8) " {c |}"
		_col(17) "Mean"
		_col(29) "diff"
		_col(36) "Lower Bound"
		_col(50) "Upper Bound"
		_col(64) "than Control?"
		_col(80) "than Control?";
		#delimit cr
		di in smcl in gr "{hline 9}{c +}{hline 82}"

		/* Output for control */
		#delimit ;
		di as smcl in gr  %8s abbrev("`lab`control'' ",8) " {c |}" in ye
		_col(12) %9.0g `control_mean' 
		_col(30) "---"
		_col(44) "---"
		_col(58) "---"
		_col(69) "---"
		_col(85) "---";
		#delimit cr

		/* Output for all treatments excluding control */
		local i = 0
		foreach group of local list{
			local i = `i' + 1
			if mod(`i',`separator')==0{
				di as smcl in gr  "{hline 9}{c +}{hline 82}"
			}
			if "`group'" == "."{
				local g = "miss"
				local shortlab = "."
			}
			else{
				local g = `group'
				local shortlab = abbrev("`lab`g''",8)
			}
			#delimit ;
			di as smcl in gr  %8s "`shortlab'" " {c |}" in ye
			_col(12) %9.0g `mean`g''  
			_col(24) %9.0g `diff`g''
			_col(38) %9.0g `lb`g''
			_col(52) %9.0g `ub`g''
			_col(69) %3s "`greater`g''"
			_col(85) %3s "`less`g''";
			#delimit cr
		}

		/* Bottom line and notes */
 		di as smcl in gr  "{hline 9}{c BT}{hline 82}"
		local shortfact = abbrev("`fact'",8)
 		di as smcl in gr  _col(6) "Diff = mean(`shortfact')-mean(control)"
 		di as smcl in gr  _col(6) "Greater than mean(control) if Diff > " %9.0g `criterion'
 		di as smcl in gr  _col(6) "Less than mean(control) if Diff < -" %9.0g `criterion'
	}

	/* Values to be returned in r() */
	
  	return clear
  	
  	return local fact_var "`fact'"
 	return local resp_var "`resp'"
	return scalar confidence = `ci'
  	return scalar side = `side'
  	return scalar df = `df'
 	return scalar groups = `k' 
	return scalar reps = `n'
 	return scalar N = `N'
  	return scalar mse = `mse'
  	return scalar d = `d'
  	return scalar criterion = `criterion'
	return scalar control = `control'
  	

end


/* Program d_table 
   Contains four matrices of the Dunnett table values
   Takes arguments side, ci, df, and k
   Returns local value d containing table value
*/
  
capture program drop d_table
program define d_table, rclass
	args side ci df k 
	
	/* assign values i and j for grabbing table value from matrix */

	/* j is number of treatments excluding control */
	local j = `k'-1
	
	/* Taking the conservative value for i (df) */

	if `df'<=20{	
		local i = `df'-4
	}
	if `df'>20 & `df'<24{
		local i = 16
	}
	if `df'>=24 & `df'<30{
		local i = 17
	}
	if `df'>=30 & `df'<40{
		local i = 18
	}
	if `df'>=40 & `df'<60{
		local i = 19
	}
	if `df'>=60 & `df'<120{
		local i = 20
	}
	if `df'==120{
		local i = 21
	}
	if `df'>120{
		local i = 22
	}
	
	/* two sided confidence intervals / alpha=0.05 */

   	if `side' == 2 & `ci' == 95{
		capture matrix clear d_2_05
		#delimit ;
		matrix input d_2_05 = 
		(2.57,3.03,3.29,3.48,3.62,3.73,3.82,3.90,3.97,4.03,4.09,4.14,4.26,4.42\
		2.45,2.86,3.10,3.26,3.39,3.49,3.57,3.64,3.71,3.76,3.81,3.86,3.97,4.11\
		2.36,2.75,2.97,3.12,3.24,3.33,3.41,3.47,3.53,3.58,3.63,3.67,3.78,3.91\
		2.31,2.67,2.88,3.02,3.13,3.22,3.29,3.35,3.41,3.46,3.50,3.54,3.64,3.76\
		2.26,2.61,2.81,2.95,3.05,3.14,3.20,3.26,3.32,3.36,3.40,3.44,3.53,3.65\
		2.23,2.57,2.76,2.89,2.99,3.07,3.14,3.19,3.24,3.29,3.33,3.36,3.45,3.57\
		2.20,2.53,2.72,2.84,2.94,3.02,3.08,3.14,3.19,3.23,3.27,3.30,3.39,3.50\
		2.18,2.50,2.68,2.81,2.90,2.98,3.04,3.09,3.14,3.18,3.22,3.25,3.34,3.45\
		2.16,2.48,2.65,2.78,2.87,2.94,3.00,3.06,3.10,3.14,3.18,3.21,3.29,3.40\
		2.14,2.46,2.63,2.75,2.84,2.91,2.97,3.02,3.07,3.11,3.14,3.18,3.26,3.36\
		2.13,2.44,2.61,2.73,2.82,2.89,2.95,3.00,3.04,3.08,3.12,3.15,3.23,3.33\
		2.12,2.42,2.59,2.71,2.80,2.87,2.92,2.97,3.02,3.06,3.09,3.12,3.20,3.30\
		2.11,2.41,2.58,2.69,2.78,2.85,2.90,2.95,3.00,3.03,3.07,3.10,3.18,3.27\
		2.10,2.40,2.56,2.68,2.76,2.83,2.89,2.94,2.98,3.01,3.05,3.08,3.16,3.25\
		2.09,2.39,2.55,2.66,2.75,2.81,2.87,2.92,2.96,3.00,3.03,3.06,3.14,3.23\
		2.09,2.38,2.54,2.65,2.73,2.80,2.86,2.90,2.95,2.98,3.02,3.05,3.12,3.22\
		2.06,2.35,2.51,2.61,2.70,2.76,2.81,2.86,2.91,2.94,2.97,3.00,3.07,3.16\
		2.04,2.32,2.47,2.58,2.66,2.72,2.77,2.82,2.86,2.89,2.92,2.95,3.02,3.11\
		2.02,2.29,2.44,2.54,2.62,2.68,2.73,2.77,2.81,2.85,2.87,2.90,2.97,3.06\
		2.00,2.27,2.41,2.51,2.58,2.64,2.69,2.73,2.77,2.80,2.83,2.86,2.92,3.00\
		1.98,2.24,2.38,2.47,2.55,2.60,2.65,2.69,2.73,2.76,2.79,2.71,2.87,2.95\
		1.96,2.21,2.35,2.44,2.51,2.57,2.61,2.65,2.69,2.72,2.74,2.77,2.83,2.91);
		#delimit cr
		local d = d_2_05[`i',`j']
		matrix drop d_2_05
     	}

	/* two sided confidence intervals / alpha=0.01 */

     	if `side' == 2 & `ci' == 99{
		capture matrix drop d_2_01
		#delimit ;
		matrix input d_2_01 =
		(4.03,4.63,4.98,5.22,5.41,5.56,5.69,5.80,5.89,5.98,6.05,6.12,6.30,6.52\
		3.71,4.21,4.51,4.71,4.87,5.00,5.10,5.20,5.28,5.35,5.41,5.47,5.62,5.81\
		3.50,3.95,4.21,4.39,4.53,4.64,4.74,4.82,4.89,4.95,5.01,5.06,5.19,5.36\
		3.36,3.77,4.00,4.17,4.29,4.40,4.48,4.56,4.62,4.68,4.73,4.78,4.90,5.05\
		3.25,3.63,3.85,4.01,4.12,4.22,4.30,4.37,4.43,4.48,4.53,4.57,4.68,4.82\
		3.17,3.53,3.74,3.88,3.99,4.08,4.16,4.22,4.28,4.33,4.37,4.42,4.52,4.65\
		3.11,3.45,3.65,3.79,3.89,3.98,4.05,4.11,4.16,4.21,4.25,4.29,4.39,4.52\
		3.05,3.39,3.58,3.71,3.81,3.89,3.96,4.02,4.07,4.12,4.16,4.19,4.29,4.41\
		3.01,3.33,3.52,3.65,3.74,3.82,3.89,3.94,3.99,4.04,4.08,4.11,4.20,4.32\
		2.98,3.29,3.47,3.59,3.69,3.76,3.83,3.88,3.93,3.97,4.01,4.05,4.13,4.24\
		2.95,3.25,3.43,3.55,3.64,3.71,3.78,3.83,3.88,3.92,3.95,3.99,4.07,4.18\
		2.92,3.22,3.39,3.51,3.60,3.67,3.73,3.78,3.83,3.87,3.91,3.94,4.02,4.13\
		2.90,3.19,3.36,3.47,3.56,3.63,3.69,3.74,3.79,3.83,3.86,3.90,3.98,4.08\
		2.88,3.17,3.33,3.44,3.53,3.60,3.66,3.71,3.75,3.79,3.83,3.86,3.94,4.04\
		2.86,3.15,3.31,3.42,3.50,3.57,3.63,3.68,3.72,3.76,3.79,3.83,3.90,4.00\
		2.85,3.13,3.29,3.40,3.48,3.55,3.60,3.65,3.69,3.73,3.77,3.80,3.87,3.97\
		2.80,3.07,3.22,3.32,3.40,3.47,3.52,3.57,3.61,3.64,3.68,3.70,3.78,3.87\
		2.75,3.01,3.15,3.25,3.33,3.39,3.44,3.49,3.52,3.56,3.59,3.62,3.69,3.78\
		2.70,2.95,3.09,3.19,3.26,3.32,3.37,3.41,3.44,3.48,3.51,3.53,3.60,3.68\
		2.66,2.90,3.03,3.12,3.19,3.25,3.29,3.33,3.37,3.40,3.42,3.45,3.51,3.59\
		2.62,2.85,2.97,3.06,3.12,3.18,3.22,3.26,3.29,3.32,3.35,3.37,3.43,3.51\
		2.58,2.79,2.92,3.00,3.06,3.11,3.15,3.19,3.22,3.25,3.27,3.29,3.35,3.42);
		#delimit cr
     		local d = d_2_01[`i',`j']	
		matrix drop d_2_01
	}

	/* one sided confidence intervals / alpha=0.05 */ 
    	if `side' == 1 & `ci' == 95{
		capture matrix drop d_1_05
		#delimit ;
		matrix input d_1_05 =
		(2.02,2.44,2.68,2.85,2.98,3.08,3.16,3.24,3.30\
		1.94,2.34,2.56,2.71,2.83,2.92,3.00,3.07,3.12\
            1.89,2.27,2.48,2.62,2.73,2.82,2.89,2.95,3.01\
		1.86,2.22,2.42,2.55,2.66,2.74,2.81,2.87,2.92\
		1.83,2.18,2.37,2.50,2.60,2.68,2.75,2.81,2.86\
		1.81,2.15,2.34,2.47,2.56,2.64,2.70,2.76,2.81\
		1.80,2.13,2.31,2.44,2.53,2.60,2.67,2.72,2.77\
		1.78,2.11,2.29,2.41,2.50,2.58,2.64,2.69,2.74\
		1.77,2.09,2.27,2.39,2.48,2.55,2.61,2.66,2.71\
		1.76,2.08,2.25,2.37,2.46,2.53,2.59,2.64,2.69\
		1.75,2.07,2.24,2.36,2.44,2.51,2.57,2.62,2.67\
		1.75,2.06,2.23,2.34,2.43,2.50,2.56,2.61,2.65\
		1.74,2.05,2.22,2.33,2.42,2.49,2.54,2.59,2.64\
		1.73,2.04,2.21,2.32,2.41,2.48,2.53,2.58,2.62\
		1.73,2.03,2.20,2.31,2.40,2.47,2.52,2.57,2.61\
		1.72,2.03,2.19,2.30,2.39,2.46,2.51,2.56,2.60\
		1.71,2.01,2.17,2.28,2.36,2.43,2.48,2.53,2.57\
		1.70,1.99,2.15,2.25,2.33,2.40,2.45,2.50,2.54\
		1.68,1.97,2.13,2.23,2.31,2.37,2.42,2.47,2.51\
		1.67,1.95,2.10,2.21,2.28,2.35,2.39,2.44,2.48\
		1.66,1.93,2.08,2.18,2.26,2.32,2.37,2.41,2.45\
		1.64,1.92,2.06,2.16,2.23,2.29,2.34,2.38,2.42);
		#delimit cr
	      local d = d_1_05[`i',`j']	
		matrix drop d_1_05
	}
	
	/* two sided confidence intervals / alpha=0.01 */ 
     	if `side' == 1 & `ci' == 99{
      	capture matrix drop d_1_01
		#delimit ;
		matrix input d_1_01 =
		(3.37,3.90,4.21,4.43,4.60,4.73,4.85,4.94,5.03\
		3.14,3.61,3.88,4.07,4.21,4.33,4.43,4.51,4.59\
		3.00,3.42,3.66,3.83,3.96,4.07,4.15,4.23,4.30\
		2.90,3.29,3.51,3.67,3.79,3.88,3.96,4.03,4.09\
		2.82,3.19,3.40,3.55,3.66,3.75,3.82,3.89,3.94\
		2.76,3.11,3.31,3.45,3.56,3.64,3.71,3.78,3.83\
		2.72,3.06,3.25,3.38,3.48,3.56,3.63,3.69,3.74\
		2.68,3.01,3.19,3.32,3.42,3.50,3.56,3.62,3.67\
		2.65,2.97,3.15,3.27,3.37,3.44,3.51,3.56,3.61\
		2.62,2.94,3.11,3.23,3.32,3.40,3.46,3.51,3.56\
		2.60,2.91,3.08,3.20,3.29,3.36,3.42,3.47,3.52\
		2.58,2.88,3.05,3.17,3.26,3.33,3.39,3.44,3.48\
		2.57,2.86,3.03,3.14,3.23,3.30,3.36,3.41,3.45\
		2.55,2.84,3.01,3.12,3.21,3.27,3.33,3.38,3.42\
		2.54,2.83,2.99,3.10,3.18,3.25,3.31,3.36,3.40\
		2.53,2.81,2.97,3.08,3.17,3.23,3.29,3.34,3.38\
		2.49,2.77,2.92,3.03,3.11,3.17,3.22,3.27,3.31\
		2.46,2.72,2.87,2.97,3.05,3.11,3.16,3.21,3.24\
		2.42,2.68,2.82,2.92,2.99,3.05,3.10,3.14,3.18\
		2.39,2.64,2.78,2.87,2.94,3.00,3.04,3.08,3.12\
		2.36,2.60,2.73,2.82,2.89,2.94,2.99,3.03,3.06\	
		2.33,2.56,2.68,2.77,2.84,2.89,2.93,2.97,3.00);
		#delimit cr
	      local d = d_1_01[`i',`j']
		matrix drop d_1_01
    	}
	
	/* return table value in r(d) */
	return clear
	return scalar d = `d'
   
end