*! 2.1.1 NJC 6 November 2004
* 2.1.0 NJC 29 October 2004
* 2.0.0 NJC 16 Feb 2003
* 1.1.2 NJC 27 Sept 2001
program anovaplot, sort
	version 8.0
	syntax [varlist(max=3 numeric default=none)] ///
 	[ , SOrt SCatter(str asis) plot(str asis) * ]

	// initial checking and picking up what -anova- leaves behind
	if "`e(cmd)'" != "anova" {
		di as err "anova estimates not found"
		exit 301
	}

	local covars "`e(varnames)'"
	if `: word count `covars'' > 3 {
		di as err "too many predictors: should be 1, 2 or 3"
		di as txt "(predictors: `covars')"
		exit 103
	}

	if "`varlist'" != "" {
		if !`: list covars === varlist' {
			di as txt "`varlist'" ///
			as err  " not a permutation of " as txt "`covars'"
			exit 498
		}
		local covars "`varlist'"
	}

	local response "`e(depvar)'"
	tokenize `response' `covars'
	args y x1 x2 x3

	quietly {
		// get fit
        	tempvar fit
		predict `fit' if e(sample)

		if "`x2'" != "" {
			tempvar touse obsno group
			gen `touse' = e(sample)
			bysort `touse' `x2' : gen byte `group' = ///
				_n == 1 & `touse'
			gen long `obsno' = _n
			replace `group' = sum(`group')
			forval i = 1 / `= `group'[_N]' {
				tempvar f
				gen `f' = `fit' if `group' == `i'
				su `obsno' if `group' == `i', meanonly
				local value = `x2'[`r(min)']
				local lbl : label (`x2') `value'
				label var `f' `"`lbl'"'
				local fits "`fits'`f' "
			}
			local sub `"`: variable label `x2''"'
			if `"`sub'"' == "" local sub "`x2'"
			local sub "sub(`"Profiles by `sub'"', place(w))"
 			if "`x3'" != "" local sub3 `sub'
			else local sub2 `sub'
		}
		else local fits "`fit'"

		// set up graph defaults
		numlist "2/`= 1 + `: word count `fits'''"
		local order "`r(numlist)'"

		local what : variable label `y'
		if `"`what'"' == "" local what "`y'"
		local ytitle "yti("`what'")"
		if "`x3'" != "" local byby "by(`x3', `sub3')"

		anova_terms
		local cont "`r(continuous)'"
		if !`: list x1 in cont' {
			levels `x1' if e(sample), local(levels)
			local xla "xla(`levels', valuelabels)"
		}
		else local ms "ms(none ..)"
  	}

	// graph
	twoway scatter `y'  `x1' if e(sample), `scatter'                 ///
	|| connected `fits' `x1' if e(sample),                           ///
	`ms' `xla' `ytitle' sort `byby' `sub2' legend(order(`order') pos(5)) ///
	`options' ///
	|| `plot'
end

/* Possibilities for future, suggested by Ken Higbee 7 Sept 2001:

For small datasets plotting the points can be helpful, but even in those cases,
I find that the data labels (or even just plain numbers) fill up the plot to
the point that you can not see the forest because of all the trees.  I think
that an easier to view plot could be produced using more of an approach like
the -serrbar- command.  I am thinking along the lines of using -adjust- which
in addition to giving you the linear predictions can also give the errors.
[...]
I would probably have some option that allowed you to get the data points
plotted and/or the error bars for each adjusted cell prediction.

With more than 3 vars, I believe to get meaningful views you need to set the
extra vars to some value (by default their mean for continuous vars and one of
the levels for a categorical var).  I would have it default to means for
continuous vars and first level for categorical vars with some syntax that
would allow users to specify if they desire something else.

*/