/************************************************************************
Author: Aluisio J D Barros
Email: abarros@equidade.org
Maintained by: Cauane Blumenberg, Leonardo Ferreira
Email: cblumenberg@equidade.org, lferreira@equidade.org

This program calculates both the slope index of inequality (SII) and the
relative index of inequality (RII) using a logistic regression model.

This approach is suitable for binary outcomes where the probability of
the outcome is estimated, instead of its mean.

The logistic model avoids adjusted values out of the zero-one interval
that sometimes happen when a linear model is used for binary outcomes.

Version history
April 2014 - corrected for the case where SII/RII were calculated even when the rank variable was missing for all obs.

Jan 2015 - adjustments made to situations where the coverage is very high (~1) or very low (~0) in several of the quintiles, results would be unstable.
           We now replace these values with .99 or .01 to improve the results. The SII is probably still an overestimate.

Jun 2017 - Sort stable is now in use to make the results be the same in every run of the routine.

Mar 2018 - A line was added to drop from the dataset missing variables for either the stratifier or outcome.

Jan 2019 - the observed probabilities for each group are now calculated using a logistic model as well, with the stratifier as a categorical variable in the model.
           This allows weights to be taken into account properly. Also, a change in the outcome from 1 to .99 and 0 to .01 was removed. The code was changed to
           include the clusters in the model using a local macro instead of if's.

Jun 2023 - simplified the equation that estimates the SII
           the command now works directly with continuous ranking variables - option continuous is now deprecated
           added an option to compare the SII measure between the groups of a third variable (limited to dichotomous variables)
           replaced all named variables to tempvars to prevent overlap
           rii is no longer calculated by default - it now requires the option rii
           p value is now available in the returned matrix
*************************************************************************/

/************************************************************************
Usage:

siilogit rankvar outcomevar [if] [in] [w=weightvar], cluster(clustervar) [nograph]

rankvar must be ordinal and ordered so that low is poor/worst and high is rich/best

The use of w=weight allows for the use a (sample) weights.
The use of cluster(clustervar) allows for clustered samples.
If the option nograph is used, the graph to check for linearity is not displayed.

Example:

siilogit v190 sba [w=v005], cluster(v021)
*************************************************************************/

*** Small program to calculate scalars with the estimate, se and p value
capture program drop calcscalars
program define calcscalars
  args indic
  mat def `indic' = e(b)
  mat def `indic'se = e(V)
  scalar `indic'_r = `indic'[1,1]
  scalar `indic'_se = sqrt(`indic'se[1,1])
  scalar `indic'_p = 2*normal(-abs(`indic'_r/`indic'_se))
end

program define siilogit, rclass
version 17.0

syntax varlist(max=2 numeric) [if] [in] [fweight/] , [COMPare(varname numeric)] [Cluster(varname numeric)] [NOGRaph] [RII]

*** Values to 1 and 2
tokenize `varlist'
tempvar _freq _rfi _rfcum _rkmidp _rkmidpalt _stratgroups _lom _lpm
preserve
quietly {
  set more off

  *** Apply "if" and "in" restrictions and remove missings for concentration variable
  gen byte _use = 1 `if' `in'
  keep if _use==1                                                               // drop observations not to be used in the temporary dataset
  drop if missing(`1',`2')                                                      // drop observations missing for either rank or outcome variables

  *** Check for variables being missing
  inspect `1'
  if r(N) == 0 {
    di as error "Error! Rank variable [" "`1'" "] is missing for all observations."
    error 416
  }
  inspect `2'
  if r(N) == 0 {
    di as error "Error! Outcome variable [" "`2'" "] is missing for all observations."
    error 416
  }
  keep if `2' ~=.
  replace `2' = .99 if `2' > .99                                                // avoid convergence problems
  replace `2' = .01 if `2' < 0.01

  *** Id rank variable
  local rkv `1'
  inspect `rkv'
  local ngroups = cond(r(N_unique)==.,99,r(N_unique))                           // r(N_unique) is missing if more than 99 groups
  local grouped = r(N)==r(N_unique)                                             // local macro grouped = 1 if no. of observations = no. of groups
  local sampsize = r(N)                                                         // sample size

  *** Create midpoints for ranked groups
  sort `rkv', stable
  if "`exp'" == "" gen `_freq'=1
  else gen `_freq' = `exp'
  sum `_freq'
  gen `_rfi' =`_freq' / r(sum)
  gen `_rfcum' = sum(`_rfi')

  if `grouped' == 1 {                                                           // GROUPED DATA
    gen `_rkmidp' =`_rfi'/2 + `_rfcum'[_n-1]
    replace `_rkmidp' = `_rfi'/2 in 1
  }
  else {                                                                        // MICRODATA
    if (`ngroups' > 20) {
      gen `_rkmidp' = `_rfcum'                                                  // continuous stratifier
      local stratngroups = max(min(100,int(`sampsize'/500)*10),20)              // define the number of groups to divide the continuous stratifier for the regression, with a maximum of 100 and a minimum of 20
      xtile `_stratgroups' = `rkv', n(`stratngroups')
      by `_stratgroups', sort: egen `_rkmidpalt' = mean(`_rfcum')
    }
    else {
      by `rkv', sort: egen `_rkmidp' = mean(`_rfcum')                           // categorical stratifier
      gen `_rkmidpalt' = `_rkmidp'
    }
  }

  *** Logistic regression fitted with weights and clusters
  describe, short                                                               // This block adds two observations
  local nobs2 = r(N)+2                                                          // including values 0 and 1 for
  local nobs1 = r(N)+1                                                          // the midpoint rank
  set obs `nobs2'                                                               // so that the graph includes a line that goes
  replace `_rkmidp'=0 in `nobs2'                                                // from zero to one
  replace `_rkmidp'=1 in `nobs1'                                                // end of block

  *** Define the condition of using clusters or not
  if "`cluster'" == "" {
    local vcecond                                                               // condition is blank
  }
  else {
    local vcecond  vce(cluster `cluster')                                       // condition is vce(cluster var)
  }

  *** Estimating fitted values for each group (as a categorical variable)
  if (`ngroups'<=20) {                                                          // stratifier is grouped
    capture glm `2' i.`rkv' [pw=`_freq'], f(bin) l(logit) `vcecond' search iterate(100)
  }
  else {                                                                        // stratifier is continuous
    capture glm `2' i.`_stratgroups' [pw=`_freq'], f(bin) l(logit) `vcecond' search iterate(100)
  }
  if _rc==0 {
    predict _lom, xb                                                            // observed values (logit probabilities)
  }


  *** Fitting the model that will be used for calculating the SII
  capture glm `2' `_rkmidp' [pw=`_freq'], f(bin) l(logit) `vcecond' search iterate(100)
  if _rc==0 {
    return scalar er = 0
    predict _lpm, xb                                                            // predicted values from the model (logit probabilities)

    *** Estimating the SII
    nlcom (SII_1: 1/(1+exp(-1*(_b[_cons] +_b[`_rkmidp'])))) (SII_0: 1/(1+exp(-1*(_b[_cons])))), post
    noi nlcom SII: _b[SII_1] - _b[SII_0], post noheader
    calcscalars sii
    return scalar sii_p = sii_p
    return scalar sii_se = sii_se
    return scalar sii = sii_r

    if "`rii'"!="" {
      *** Estimating the RII
	    capture glm `2' `_rkmidp' [pw=`_freq'], f(bin) l(logit) `vcecond' search iterate(100)
	    nlcom (RII_1: 1/(1+exp(-1*(_b[_cons] +_b[`_rkmidp'])))) (RII_0: 1/(1+exp(-1*(_b[_cons])))), post
	    noi nlcom RII: _b[RII_1] / _b[RII_0], post noheader
      calcscalars rii
      return scalar rii_p = rii_p
      return scalar rii_se = rii_se
      return scalar rii = rii_r
    }

    *** Graph to check for linearity
    if "`nograph'"=="" {
      noisily twoway sc _lom `_rkmidpalt', s(o) || sc _lpm  `_rkmidp', c(l) s(i) xlabel(0 (.1) 1) legend(label(1 "Observed") label(2 "Fitted")) xtitle(Fractional rank) ytitle(Logit probability)
    }

    *** Compare SII/RII between two subgroups
    if "`compare'"!="" {
      inspect `compare'
      if r(N) == 0 {
        di as error "Error! Comparison variable [" "`compare'" "] is missing for all observations."
        error 416
      }
      if r(N_unique) != 2 {
        di as error "Error! Comparison variable [" "`compare'" "] has more than 2 groups."
        error 416
      }
      levelsof `compare', local(complevels)
      local lastval: word 2 of `complevels'                                     // store the second value of compare variable
      local vallab : value label `compare'
      if "`vallab'" != "" {                                                     // if the variable has value labels use them
        local i 1
        foreach level of local complevels {
          local lab`i' : lab `vallab' `level'
          local lab`i' = trim(subinstr("`lab`i''", " ", "_",.))
          local ++i
        }
      }
      else {                                                                    // if the variable has not value labels, use the actual numeric values
        local i 1
        foreach level of local complevels {
          local lab`i' `level'
          local lab`i' = trim(subinstr("`lab`i''", " ", "_",.))
          local ++i
        }
      }

      glm `2' c.`_rkmidp'##i.`compare' [pw=`_freq'], f(bin) l(logit) `vcecond' search iterate(100)
      * Pred_10 is estimate at top of wealth for level 0, Pred_11 is estimate at top wealth for level 1, etc.
      nlcom (Pred_10: 1/(1+exp(-1*(_b[_cons] + _b[`_rkmidp'])))) (Pred_00: 1/(1+exp(-1*(_b[_cons])))) ///
            (Pred_11: 1/(1+exp(-1*(_b[_cons] + _b[`_rkmidp'] + _b[`lastval'.`compare'] + _b[`lastval'.`compare'#c.`_rkmidp'])))) ///
            (Pred_01: 1/(1+exp(-1*(_b[_cons] + _b[`lastval'.`compare'])))), post noheader
      noisily display _newline "SII values for `lab1' and `lab2'"
      noisily nlcom (SII_`lab1': _b[Pred_10] - _b[Pred_00]) (SII_`lab2': _b[Pred_11] - _b[Pred_01]), post noheader
      noisily display _newline "Difference between SII in `lab2' and `lab1'"
      noisily nlcom (SII_compare: _b[SII_`lab2'] - _b[SII_`lab1']), post noheader
      calcscalars siicomp
      return scalar siicomp_p = siicomp_p
      return scalar siicomp_se = siicomp_se
      return scalar siicomp = siicomp_r

      if "`rii'"!="" {
	      glm `2' c.`_rkmidp'##i.`compare' [pw=`_freq'], f(bin) l(logit) `vcecond' search iterate(100)
	      * Pred_10 is estimate at top of wealth for level 0, Pred_11 is estimate at top wealth for level 1, etc.
	      nlcom (RII_10: 1/(1+exp(-1*(_b[_cons] + _b[`_rkmidp'])))) (RII_00: 1/(1+exp(-1*(_b[_cons])))) ///
	            (RII_11: 1/(1+exp(-1*(_b[_cons] + _b[`_rkmidp'] + _b[2.`compare'] + _b[2.`compare'#c.`_rkmidp'])))) ///
              (RII_01: 1/(1+exp(-1*(_b[_cons] + _b[2.`compare'])))), post noheader
	      noisily display _newline "RII values for `lab1' and `lab2'"
	      noisily nlcom (RII_`lab1': _b[RII_10] / _b[RII_00]) (RII_`lab2': _b[RII_11] / _b[RII_01]), post noheader
	      noisily display _newline "Ratio between RII in `lab2' and `lab1'"
	      noisily nlcom (RII_compare: _b[RII_`lab2'] / _b[RII_`lab1']), post noheader
        calcscalars riicomp
        return scalar riicomp_p = riicomp_p
        return scalar riicomp_se = riicomp_se
        return scalar riicomp = riicomp_r
  	  }
    }
  }
  else {
    noisily display _newline as error "Error calculating SII/RII" _newline
    return scalar er = 1
    error 999
  }
} // end quietly

end // program