*! version 2.0.4 JMGarrett 14May98  STB-43 sg33.1
/* Program to calculate adjusted probabilities for nominal variables       */
/* Form:  adjprop y, by(x1 x2) adj(cov1 cov2 ...) options                  */
/* Options required:  by (x1 [x2])   <-- nominal variables only            */
/* Options allowed:  adjust, model, level, graph, bar, graph_options       */

program define adjprop
  version 4.0
  #delimit ;
    local options "BY(string) Adjust(string) Model Graph L2title(string)
      L1title(string) T2title(string) Bar Level(real 95) *" ;
  #delimit cr
  local varlist "req ex min(1) max(1)"
  local if "opt"
  parse "`*'"
  parse "`varlist'", parse(" ")
  preserve
  capture keep `if'
  local yvar="`1'"
  quietly drop if `yvar'==.

* If there are covariates, drop missing values, calculate means
  parse "`adjust'", parse(" ")  
  local numcov=0
  local i=1
  while "`1'"~="" {
    local equal=index("`1'","=")
    if `equal'==0  {
       local cov`i'="`1'"
       local mcov`i'="mean"
       }
    if `equal'~=0  {
       local cov`i'=substr("`1'",1,`equal'-1)
       local mcov`i'=substr("`1'",`equal'+1,length("`1'"))
       }
    quietly drop if `cov`i''==.
    local covlist `covlist' `cov`i''
    local covdisp `covdisp' `1'
    local i=`i'+1
    macro shift
    local numcov=`i'-1
    }
  local i=1
  while `i'<=`numcov' {
    if "`mcov`i''"=="mean" {
      quietly sum `cov`i''
      local mcov`i'=_result(3)
      }
    local i=`i'+1
    }
  keep `yvar' `by' `covlist'

* Read in X variables and create dummy variables
  parse "`by'", parse(" ")
  local xvar1="`1'"
  local vlblx1 : variable label `xvar1'
  quietly drop if `xvar1'==.
  quietly tab `xvar1', gen(X)
  local numcat1=_result(2)
  local i=2
  while `i'<=`numcat1'  {
    local xlist1 `xlist1' X`i'
    local i=`i'+1
    }
  macro shift
  if "`1'"==""  {
    local x2=0
    local numcat=`numcat1'
    }
  if "`1'"~=""  {
    local xlist1 ""
    local i=2
    while `i'<=`numcat1'  {
      rename X`i' X1`i'
      local xlist1 `xlist1' X1`i'
      local i=`i'+1
      }
    local x2=1
    local xvar2="`1'"
    quietly drop if `xvar2'==.
    local vlblx2 : variable label `xvar2'
    quietly tab `xvar2', gen(X2)
    local numcat2=_result(2)
      local i=2
      while `i'<=`numcat2'  {
        local xlist2 `xlist2' X2`i'
        local i=`i'+1
        }
    local numcat=`numcat1'*`numcat2'
  macro shift
  if "`1'"~=""  {
    disp in red "Only two X variables allowed in the by( ) option"
    exit
    }
    
    * create interaction terms
    local i=2
    while `i'<=`numcat1' {
      local j=2
      while `j'<=`numcat2' {
        quietly gen I`i'`j'=X1`i'*X2`j'
        local intlist `intlist' I`i'`j'
        local j=`j'+1
        }
      local i=`i'+1
      }
    }

* Run logistic regression to get parameter estimates
   if "`model'"~="model"  {
     quietly logistic `yvar' `xlist1' `xlist2' `intlist' `covlist'
     }
   if "`model'"=="model"  {
     logistic `yvar' `xlist1' `xlist2' `intlist' `covlist'
     }
  local varlbly : variable label $S_E_depv
  estimates hold logest
  
* LR test for overall association, and interaction if present
  local l1=-2*_result(2)
  quietly logistic `yvar' `covlist'
  local l2=-2*_result(2)
  local df=`numcat'-1
  local chisq=`l2'-`l1'
  local probchi=chiprob(`df',`chisq')
  if `x2'==1  {
    quietly logistic `yvar' `xlist1' `xlist2' `covlist'
    local l2=-2*_result(2)
    local dfi=(`numcat1'-1)*(`numcat2'-1)
    local chisqi=`l2'-`l1'
    local pchisqi=chiprob(`dfi', `chisqi')
    }
  estimates unhold logest
   
* Collapse to 1 obs. per category, dummy variables, and covariates
  tempvar count
  quietly gen `count'=1
  sort `by'
  #delimit ;
    collapse `count' `yvar' `xlist1' `xlist2' `intlist', by(`by')
      max(. `yvar' `xlist1' `xlist2' `intlist') sum(numobs) ;
  #delimit cr
  quietly replace $S_E_depv=.

* Replace covariates with their means (or specified values)
  local i=1
  while `i'<=`numcov'  {
    quietly gen `cov`i''=`mcov`i''
    local i=`i'+1
    }

* Calculate the adjusted probabilities (risks) and confidence intervals
  tempvar linpred
  predict adjprob
  predict se, stdp
  predict `linpred', xb
  local z=invnorm((1-`level'/100)/2)
  gen upper=1/(1+exp(-`linpred'+`z'*se))
  gen lower=1/(1+exp(-`linpred'-`z'*se))

* Print results
  display "   "
  if `numcov'>0  {
    #delimit ;
      display in yel "*" in green "Adjusted" in yel "*" in green
        " Probabilities and `level'% Confidence Intervals" ;
    #delimit cr
    local probtyp="adjprob"
    }
  if `numcov'==0  {
    #delimit ;
      display in yel "*" in green "Unadjusted" in yel "*" in green
        " Probabilities and `level'% Confidence Intervals";
    #delimit cr
    local probtyp="prob"
    quietly gen prob=adjprob
    }
  display "  "
  display in gr "  Outcome:" in yel "      `varlbly' -- $S_E_depv"
  if `x2'==0 {
    display in gr "  Nominal X:" in yel "    `vlblx1' -- `xvar1'"
    }
  if `x2'==1 {
    display in gr "  Nominal X1:" in yel "   `vlblx1' -- `xvar1'"
    display in gr "  Nominal X2:" in yel "   `vlblx2' -- `xvar2'"
    display in gr "  Interaction:" in yel "  `xvar1' * `xvar2'"
    }
  if `numcov'~=0 {display in gr "  Covariates:" in yel "   `covdisp'"}
  if `numcov'==0 {display in gr "  Covariates:" in yel "   (none)"}
  list `by' numobs `probtyp' se lower upper, noob nod
  disp "  "
  disp in gr "  Likelihood ratio test for difference of `numcat' probabilities:"
  disp "  "
  disp in green "    LR Chi2(`df')    =  " in yellow %6.2f `chisq'
  if `probchi'>=.0001 {
    disp in green "    Prob > Chi2   =   " in yellow %7.4f `probchi'
    }
  if `probchi'<.0001 {
    disp in green "    Prob > Chi2   <   " in yellow "0.0001"
    }
  if `x2'==1  {
    disp "  "
    #delimit ;
      disp in gr "  Likelihood ratio test of interaction for" in yellow 
      " `xvar1' * `xvar2'" in green ":";
    #delimit cr
    disp "  "
    disp in green "    LR Chi2(`dfi')    =  " in yellow %6.2f `chisqi'
    if `pchisqi'>=.0001 {
      disp in green "    Prob > Chi2   =   " in yellow %7.4f `pchisqi'
      }
    if `pchisqi'<.0001 {
      disp in green "    Prob > Chi2   <   " in yellow "0.0001"
      }
    }

* Graph the results, if requested
  if "`graph'"=="graph"  {
    more
    if "`l2title'"=="" {local l2title "`varlbly' -- $S_E_depv"}
    if "`l1title'"=="" & `numcov'==0 {
      if "`bar'"=="" {
        local l1title "Unadjusted Probabilities and `level'% CI"
        }
      if "`bar'"=="bar" { 
        local l1title "Unadjusted Probabilities"
        } 
      }
    if "`l1title'"=="" & `numcov'>0 {
     if "`bar'"=="" {
        local l1title "Adjusted Probabilities and `level'% CI"
        }
      if "`bar'"=="bar" { 
        local l1title "Adjusted Probabilities"
        } 
      }
    if `x2'==0 {
      if "`bar'"=="" {
      #delimit ;
        graph adjprob lower upper `xvar1', c(.II) s(Oii) `options'
           l2("`l2title'") l1("`l1title'") ;
      #delimit cr
      }
      if "`bar'"=="bar" {
       if "`t2title'"=="" & "`vlblx1'"=="" {local t2title "`xvar1'"}
       if "`t2title'"=="" & "`vlblx1'"~="" {
         local t2title "`vlblx1' -- `xvar1'"
         }
       sort `xvar1'
       #delimit ;
         graph adjprob, by(`xvar1') bar `options' l2("`l2title'") 
             l1("`l1title'") t2("`t2title'") ;
       #delimit cr
       }
      }
    if `x2'==1 {
      local x2val : value label `xvar2'
      sort `xvar2' `xvar1'
      local i=1
      local n=1
      while `n'<=_N  {
        if `xvar2'[`n']~=`xvar2'[`n'-1] {
          local catv`i'=`xvar2'[`n']
          local catv `catv' `catv`i''
          local i=`i'+1
          }
        local n=`n'+1
        }
       if "`x2val'"~="" {
         local i=1 
           while `i'<=`numcat2' {
             local x2lbl`i' : label `x2val' `catv`i''
             local i=`i'+1
             }
          }
       rename adjprob _M
       reshape groups `xvar2' `catv'
       reshape vars _M
       reshape con `xvar1'
       reshape wide
       local i=1
       while `i'<=`numcat2'  {
         if "`x2val'"~="" {label var _M`catv`i'' "`xvar2' = `x2lbl`i''"}
           else           {label var _M`catv`i'' "`xvar2' = `catv`i''"}
         local i=`i'+1
         }
       if "`bar'"=="bar" {
        if "`l1title'"=="" & `numcov'==0 {
          local l1title "Unadjusted Probabilities" 
          }
        if "`l1title'"=="" & `numcov'>0 {
          local l1title "Adjusted Probabilities" 
          }
        sort `xvar1'
        graph _M*, by(`xvar1') bar `options' l2("`l2title'") l1("`l1title'")
        }
      if "`bar'"~="bar" {
        graph _M* `xvar1', `options' l2("`l2title'") l1("`l1title'")
        }
      }
    }
end