*! version 3.0.2 JMGarrett 18Oct03
/* Graph reliablity plot of predictions from logistic regression model   */
/* Form:  relyplot, groups(#) or fractions(3, 4, 5, or 10)               */
/* Note: Uses estimates from prior run logistic regression model         */
/* Note: Renamed from "rely" to "relyplot"  19Nov17                      */  
program relyplot
  version 8.0
  syntax [if] [in], [, GRoups(int 0) FRActions(int 0) CI * ]
  marksample touse
  preserve
  capture assert `groups'>0 | `fractions'>0
  if _rc~=0 {
     disp as error " Either groups(#) or fractions(#) must be specified"
     exit
     }
  if `groups'>0 & `fractions'>0  {
     disp as error "Choose either groups(#) or fractions(#) -- both"  /*
       */ " should not be specified."
     exit
     }
  if `fractions'~=0 {
     capture assert `fractions'==3 | `fractions'==4 | `fractions'==5 | /*
          */ `fractions'==10
     if _rc~=0 {
        disp as error "Fractions can be 3, 4, 5, or 10" 
        exit 
        }
     }
  local d=e(depvar)
  quietly predict pred_d if `touse', p
  quietly keep if pred_d~=.
  keep `d' pred_d
  quietly gen x=1

* Create deciles of prediced risk; calculate observed proportions
  if `fractions'>0 {
    if `fractions'==10 { 
      quietly gen pred_grp=recode(pred_d,.1,.2,.3,.4,.5,.6,.7,.8,.9,1)
      quietly replace pred_grp=round(pred_grp*1000-50,1)
      local pclab "tenths"
      local pclabgr "Tenths"
      #delimit ;
        label define predlbl 50 "0-.1" 150 ".1-.2" 250 ".2-.3"
           350 ".3-.4" 450 ".4-.5" 550 ".5-.6" 650 ".6-.7"
           750 ".7-.8" 850 ".8-.9" 950 ".9-1.0" 1000 "1.0";
      #delimit cr
      local xlabval "50(100)950"
      } 
    if `fractions'==5 {
      quietly gen pred_grp=recode(pred_d,.2,.4,.6,.8,1)
      quietly replace pred_grp=round(pred_grp*1000-100,1)
      local pclab "fifths"
      local pclabgr "Fifths"
      #delimit ;
        label define predlbl 100 "0-.2" 300 ".2-.4" 500 ".4-.6"
           700 ".6-.8" 900 ".8-1.0" ;
      #delimit cr
      local xlabval "100 300 500 700 900"
        }
    if `fractions'==4 {
      qui gen pred_grp=recode(pred_d,.25,.5,.75,1)
      quietly replace pred_grp=round(pred_grp*1000-125,1)
      local pclab "forths"
      local pclabgr "Forths"
      #delimit ;
        label define predlbl 125 "0-.25" 375 ".25-.5" 625 ".5-.75"
           875 ".75-1.0" ;
      #delimit cr
      local xlabval "125(250)875"
      }
    if `fractions'==3 {
      qui gen pred_grp=recode(pred_d,.33,.67,1)
      quietly replace pred_grp=round(pred_grp*1000-165,1)
      local pclab "thirds"
      local pclabgr "Thirds"
      label define predlbl 165 "0-.33" 505 ".33-.67" 835 ".67-1.0"
      local xlabval "165 505 835"
      }
    sort pred_grp
    collapse (mean) obs_prop="`d'" (sum) n=x obs_d="`d'", by(pred_grp)
    }

* Create percentiles of predicted risk; calculate observed proportions
  if `groups'>0 {
    sort pred_d
    quietly gen pred_grp=group(`groups')
    quietly egen max=max(pred_d), by(pred_grp)
    quietly replace max=max[_n-1] if pred_d==pred_d[_n-1]
    #delimit  ;
      collapse (mean) obs_prop="`d'" (min) min_pred=pred_d
         (max) max_pred=pred_d (count) n=pred_d (sum) obs_d="`d'",
          by(pred_grp) ;
    #delimit cr
    quietly gen diffpred=(max_pred-min_pred)/2
    quietly replace pred_grp=min_pred+diffpred
    quietly egen maxpred=max(max_pred)
    quietly egen minpred=min(min_pred)
    quietly egen maxobs=max(obs_prop)
    quietly egen minobs=min(obs_prop)
    local maxmax=max(maxpred,maxobs)
    local minmin=min(minpred,minobs)
    drop maxpred minpred maxobs minobs
    }

* Generate the binomial 95% CI for observed values
  if "`ci'"=="ci" {
     tempvar cise
     qui gen `cise'=sqrt(obs_prop*(1-obs_prop)/n)
     qui gen upper=obs_prop+1.96*`cise'
     qui gen lower=obs_prop-1.96*`cise'
     }
 
* Add two observations for diagonal line on plot
  tempfile diag
  quietly save `diag'
  drop _all
  quietly set obs 2
  if `fractions'>0 {
    quietly gen pred_grp=0 if _n==1
    quietly replace pred_grp=1000 if _n==2
    quietly gen diagy=0 if _n==1
    quietly replace diagy=1 if _n==2
    }
  if `groups'>0 {
    quietly gen pred_grp=`minmin' if _n==1
    quietly replace pred_grp=`maxmax' if _n==2
    quietly gen diagy=`minmin' if _n==1
    quietly replace diagy=`maxmax' if _n==2
    }
  append using `diag'
  label var obs_prop "Observed Proportions"
  label var diagy "Perfect Fit"
 
* Print prediction values and observed proportions
  format obs_prop %9.3f
  if "`ci'"=="ci" {
    format upper lower %9.3f
    }
  if `fractions'>0 {
    label define predlbl 1000 "1.0", modify 
    label val pred_grp predlbl
    disp " "
    disp as text "Reliability Plots: Predicted risks for " /*
       */ as result "`d'" as text " divided"
    disp as text "                      into " as result /*
       */ "`pclab'" as text " of predicted risk"
    disp as text "
    disp " "
    if "`ci'"=="ci" {
       list pred_grp obs_d n obs_prop lower upper     ///
          if pred_grp~=. & n~=., separator(0) noob
       }
    if "`ci'"=="" {
       list pred_grp obs_d n obs_prop if pred_grp~=. & n~=., ///
       separator(0) noob
       }
    }
  if `groups'>0 {
    format pred_grp max_pred min_pred %9.3f
    disp " "
    disp as text "Reliability Plots: Predicted risks for " /*
       */ as result "`d'" as text " divided"
    disp as text "                       into " as result /* 
       */ "`groups'" as text " equal size groups"
    disp " "
    if "`ci'"~="" {
      list pred_grp min_pred max_pred obs_d n obs_prop lower upper ///
         if n~=., separator(0) noob
      }
    if "`ci'"=="" {
      list pred_grp min_pred max_pred obs_d n obs_prop if n~=.,  ///
         separator(0) noob
      }
    }

* Graph results
  format obs_prop %10.2f
  format pred_grp %10.2f
  local title "Reliability of Predicted Values"

  if `fractions'>0 {
    lab var pred_grp "Predicted Risk Divided into `pclabgr'"
    local xtitle "Predicted Risk Divided into `pclabgr'"
    if "`ci'"~="" {
       twoway (scatter obs_prop pred_grp, sort)                    ///
              (rcap upper lower pred_grp, sort)                    ///            
              (line diagy pred_grp, sort),                         ///
              ytitle(Observed Proportions) yscale(range(. .))      ///
              xtitle("`xtitle'") title("`title'")                  ///
              xlabel(`xlabval',valuelabel)                         ///
              legend(order(1 2 "95% CI" 3) rows(1))                ///
              ylabel(, angle(horizontal)) `options'
       }
    if "`ci'"=="" {
       twoway (scatter obs_prop pred_grp, sort)                    ///
              (line diagy pred_grp, sort),                         ///
              ytitle(Observed Proportions) yscale(range(. .))      ///
              xtitle("`xtitle'") title("`title'")                  ///
              xlabel(`xlabval',valuelabel)                         ///
              ylabel(, angle(horizontal)) `options'
       }

  }
 
 if `groups'>0 {
    if "`ci'"~="" {
      twoway (scatter obs_prop pred_grp, sort)                    ///
             (rcap upper lower pred_grp, sort)                    ///            
             (line diagy pred_grp, sort),                         ///
             ytitle(Observed Proportions) yscale(range(. .))      ///
             xtitle("`xtitle'") title("`title'")                  ///
             legend(order(1 2 "95% CI" 3) rows(1))                ///
             ylabel(, angle(horizontal)) `options'
      }
    if "`ci'"=="" {
      twoway (scatter obs_prop pred_grp, sort)                    ///
             (line diagy pred_grp, sort),                         ///
             ytitle(Observed Proportions) yscale(range(. .))      ///
             xtitle("`xtitle'") title("`title'")                  ///
             ylabel(, angle(horizontal)) `options'
      }
   }
end