/* Program to calculate adjusted means for nominal variables in regression */
/* 01/25/94  JMGarrett  (version 2)                                        */
/* Form:  adjmean y, by(x) adj(cov1 cov2 ...) options                      */
/* Options required:  by                                                   */
/* Options allowed:  adjust, model, plot, ylabel, xlabel, titles           */
/* Note:  For nominal variables only                                       */
/* Note:  Does not need to follow regression procedure                     */
/*      STB-24: sg33                                                       */
program define adjmean
  version 4.0
  local options "BY(string) ADJust(string) Model Plot L2title(string) *"
  local varlist "req ex min(1) max(1)"
  local if "opt"
  parse "`*'"
  parse "`varlist'", parse(" ")
  preserve
  capture keep `if'
  keep `varlist' `by' `adjust'
  tempvar count
  quietly gen `count'=1
  local yvar="`1'"
  local xvar="`by'"
  local varlblx : variable label `by'
  quietly drop if `xvar'==. | `yvar'==.

* create dummy variables from x and create variable list
  quietly tab `xvar', gen(x)
  local numcat=_result(2)
  local i=1
  local varlst="`1'"
  while `i'<`numcat'  {
    local varlst="`varlst'" + " " + "x`i'"
    local i=`i'+1
    }

* If there are covariates, drop missing values, set up list, calculate means
  parse "`adjust'", parse(" ")  
  local numcov=0
  local i=1
  while "`1'"~="" {
     local cov`i'="`1'"
     quietly drop if `cov`i''==.
     local i=`i'+1
     local covlst "`covlst' `1'"
     macro shift
     local numcov=`i'-1
     }
  local i=1
  while `i'<=`numcov'  {
     quietly sum `cov`i''
     local mcov`i'=_result(3)
     local i=`i'+1
     }

* Run regression to get parameter estimates
  if "`model'"=="model"  {
     reg `varlst' `covlst'
     more
     }
  if "`model'"~="model"  {
     quietly reg `varlst' `covlst'
     }
  local varlbly : variable label $S_E_depv
   
* Collapse to 1 obs. per category, dummy variables, and covariates
  sort `xvar'
  collapse "`count'" "`varlst'", by(`xvar') max(. "`varlst'") sum(numobs)
  quietly replace $S_E_depv=.
  quietly drop if `xvar'==.
   
* Replace covariates with their means
  local i=1
  while `i'<=`numcov'  {
    quietly gen `cov`i''=`mcov`i''
    local i=`i'+1
    }

* Calculate the adjusted means and 95% confidence intervals
  tempvar se linpred
  predict adjmean
  predict `se', stdp
  predict `linpred', xb
  gen lower=`linpred'-1.96*`se'
  gen upper=`linpred'+1.96*`se'
  
* Plot the results if requested and print results
  if "`plot'"=="plot"  {
     if "`l2title'"=="" {local l2title "`varlbly' -- $S_E_depv"} 
     label var adjmean "Estimated Mean and 95% CI"
     #delimit ;
       graph adjmean lower upper `xvar', c(.II) s(Oii) `options'
          l2("`l2title'")  ;
     #delimit cr
     }
  display "   "
  if `numcov'>0  {
     display "Adjusted Means and 95% Confidence Intervals" 
     local meantyp="adjmean"
     }
  if `numcov'==0  {
     display "Unadjusted Means and 95% Confidence Intervals"
     local meantyp="mean"
     rename adjmean mean
     }
  display "  "
  display "  Outcome variable:      `varlbly' -- $S_E_depv"
  display "  Categorical variable:  `varlblx' -- `xvar'"
  display "  Covariates:  `covlst'"
  list `xvar' numobs `meantyp' lower upper
end