/* Program to calculate adjusted proportions for nominal variables         */
/* 01/24/95  JMGarrett (version 2)                                         */
/* Form:  adjprop y, by(x) adj(cov1 cov2 ...) options                      */
/* Options requires:  by                                                   */
/* Options allowed:  model plot ylabel xlabel titles                       */
/* Note:  For nominal x variable only, logistic regression model           */
/* Note:  Does not need to follow logistic procedure                       */
/*      STB-24: sg33                                                       */
program define adjprop
  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 logistic regression to get parameter estimates
  if "`model'"=="model" {
     logistic `varlst' `covlst'
     more
     }
  if "`model'"~="model"  {
     quietly logistic `varlst' `covlst'
     }
  local varlbly : variable label $S_E_depv
   
* Collapse to 1 obs. per category, dummy variables, covariates, counts
  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 proportions and 95% confidence intervals
  tempvar se linpred
  predict adjprop
  predict `se', stdp
  predict `linpred', xb
  gen lower=1/(1+exp(-`linpred'+1.96*`se'))
  gen upper=1/(1+exp(-`linpred'-1.96*`se'))
  
* Plot the results if requested and print the results  
  if "`plot'"=="plot"  {
     if "`l2title'"==""  {local l2title "`varlbly' -- $S_E_depv"}
     label var adjprop "Estimated Proportion and 95% CI"
     #delimit ;
       graph adjprop lower upper `xvar', c(.II) s(Oii) `options'
         l2("`l2title'") ;
     #delimit cr
     }
  display "   "
  if `numcov'>0  {
      display "Adjusted Proportions and 95% Confidence Intervals" 
      local proptyp="adjprop"
      }
  if `numcov'==0  {
      display "Unadjusted Proportions and 95% Confidence Intervals"
      local proptyp="prop"
      rename adjprop prop
      }
  display "  "
  display "  Outcome variable:      `varlbly' -- $S_E_depv"
  display "  Categorical variable:  `varlblx' -- `xvar'"
  display "  Covariates:  `covlst'"
  list `xvar' numobs `proptyp' lower upper
end