*! version 3.1.5  12oct2024 by alexis.dinno@pdx.edu
*! perform two one-sided tests for 0th-order stochastic equivalence in unpaired
*! data

* Syntax:  tostranksum varname [if exp] [in range] , by(varname) [
*          eqvtype(type) eqvlevel(#) uppereqvlevel(#) ccontinuity alpha(#)
*          relevance]

program define tostranksum, rclass byable(recall)
  syntax varname [if] [in], BY(varname) [ EQVType(string) /*
*/       EQVLevel(real 1) UPPEReqvlevel(real 0) CContinuity Alpha(real 0.05) /*
*/       RELevance]
  if int(_caller())<8 {
    di in r "tostpr- does not support this version of Stata." _newline
    di as txt "All requests for specific version support are welcome and will be considered."
    exit
  }
  if int(_caller())<14 {
    version 8.0, missing
    * Create output constants
    local alpha_out   = "alpha"
    local delta           = "Delta"
    local delta_lower     = "Dl"
    local delta_upper     = "Du"
    local delta_inflate   = ""
    local epsilon         = "epsilon"
    local epsilon_lower   = "el"
    local epsilon_upper   = "eu"
    local epsilon_inflate = ""
    local lessthan        = "<="
    local greaterthan     = ">="
    local times           = "*"
    }
   else {
    version 14.0, missing
    * Create unicode output constants
    local alpha_out   = uchar(0945)
    local delta           = uchar(0916)
    local delta_lower     = uchar(0916)+uchar(0108)
    local delta_upper     = uchar(0916)+uchar(0117)
    local delta_inflate   = "    "
    local epsilon         = uchar(0949)
    local epsilon_lower   = uchar(0949)+uchar(0108)
    local epsilon_upper   = uchar(0949)+uchar(0117)
    local epsilon_inflate = "      "
    local lessthan        = uchar(8804)
    local greaterthan     = uchar(8805)
    local times           = uchar(0215)
    }

  * Validate eqvtype
  if lower("`eqvtype'") == "" {
    local eqvtype = "epsilon"
    }

  if !(lower("`eqvtype'") == "delta" | lower("`eqvtype'") == "epsilon") {
    noi: di as err "option eqvtype() must be either delta or epsilon"
    exit 198
    }

  * Validate eqvlevel
  if (lower("`eqvtype'") == "delta") & (`eqvlevel' == 1 & `uppereqvlevel'==0) {
    local eqvlevel = 1
    }

  if (lower("`eqvtype'") == "epsilon") & (`eqvlevel' == 1 & `uppereqvlevel'==0) {
    local eqvlevel = 2
    }

  if (lower("`eqvtype'") == "delta" || lower("`eqvtype'") == "epsilon") & (`eqvlevel' <= 0 & `uppereqvlevel' != abs(`eqvlevel')) {
    noi: di as err "option eqvlevel() incorrectly specified" _newline "the tolerance must be a positive real value"
    exit 198
    }

  * Validate uppereqvlevel
  if (`uppereqvlevel'<0) {
    noi: di as err "option uppereqvlevel() must be a positive real value"
    exit 198
    }

  if (`uppereqvlevel'==0 | `uppereqvlevel' == abs(`eqvlevel')) {
    local upper = abs(`eqvlevel')
    local lower = abs(`eqvlevel')
    }

  if (`uppereqvlevel'>0) {
    local upper = abs(`uppereqvlevel')
    local lower = abs(`eqvlevel')
    }

  * Validate alpha
  if (`alpha' < 0 | `alpha' > 1) {
    noi: di as err "option alpha() must be between 0 and 1 inclusive"
    exit 198
    }


  local PositivistConclusion = "Reject"
  local NegativistConclusion = "Reject"
  local origby "`by'"
  capture confirm numeric variable `by'
  if _rc {
    tempvar numby
    encode `by', generate(`numby')
    local by "`numby'"
  }
  marksample touse
  markout `touse' `by'

  local x "`varlist'"
  tempname g1 g2 W v unv z
  tempvar ranks
  quietly {
    summarize `by' if `touse', meanonly
    if r(N) == 0 {
      noi error 2000
      }
    if r(min) == r(max) {
      di in red "1 group found, 2 required"
      exit 499
    }
    scalar `g1' = r(min)
    scalar `g2' = r(max)

    count if `by'!=`g1' & `by'!=`g2' & `touse'
    if r(N) != 0 {
      di in red "more than 2 groups found, only 2 allowed"
      exit 499
    }

    egen double `ranks' = rank(`x') if `touse'

    summarize `ranks' if `by'==`g1' & `touse', meanonly
    local   n1   = r(N)
    scalar `W'  = r(sum)

    summarize `ranks' if `touse'
    local   n    = r(N)
    local   n2   = `n' - `n1'
    scalar `v'   = `n1'*`n2'*r(Var)/`n'
    scalar `unv' = `n1'*`n2'*(`n'+1)/12

    if `unv' < 1e7 {
      local vfmt `"%10.2f"'
      }
     else {
      local vfmt `"%10.0g"'
      }
    local se = sqrt(`v')
    local EW = `n1'*(`n'+1)/2
    if (lower("`ccontinuity'") == "") {
      local cc = 0
      local continuity_notice  di _continue
      }
    if (lower("`ccontinuity'") != "") {
      local cc = 0.5
      local continuity_notice  di as txt "Using continuity correction" _newline
      }
    local z_pos = (sign(`W'-`EW')*(abs(`W'-`EW') - `cc') )/`se'
    local p_pos = 1 - normal(abs(`z_pos'))
    local p_pos_display as txt "= " as res %6.4f `p_pos'
    if `p_pos' < 0.0001 {
      local p_pos_display as txt "< " as res "0.0001"
      }
    if `p_pos' > 0.9999 {
      local p_pos_display as txt "> " as res "0.9999"
      }
    if ( `p_pos' > `alpha'/2) {
      local PositivistConclusion = "Fail to reject"
      }
    if lower("`eqvtype'") == "delta" {
        local z1 = (`upper' - (sign(`W'-`EW')*(abs(`W'-`EW') - `cc')))/`se'
        local z2 = ((sign(`W'-`EW')*(abs(`W'-`EW') - `cc')) + `lower')/`se'
        }
    if lower("`eqvtype'") == "epsilon" {
        local z1 = `upper' - ( (sign(`W'-`EW')*(abs(`W'-`EW') - `cc'))/`se' )
        local z2 = ( (sign(`W'-`EW')*(abs(`W'-`EW') - `cc'))/`se' ) + `lower'
        }
    local neg_z_stats_out noi: di as txt _col(9) "z1 = " as res %-8.4g `z1' as txt _col(38) "z2 = " as res %-8.4g `z2' _newline
    local p1 = 1 - normal(`z1')
    local p2 = 1 - normal(`z2')
    local p1_display as txt "= " as res %6.4f `p1'
    if `p1' < 0.0001 {
      local p1_display as txt "< " as res "0.0001"
      }
    if `p1' > 0.9999 {
      local p1_display as txt "> " as res "0.9999"
      }
    local p2_display as txt "= " as res %6.4f `p2'
    if `p2' < 0.0001 {
      local p2_display as txt "< " as res "0.0001"
      }
    if `p2' > 0.9999 {
      local p2_display as txt "> " as res "0.9999"
      }
    if (`p1' > `alpha' | `p2' > `alpha') {
      local NegativistConclusion = "Fail to reject"
      }
*    }
    local neg_p_vals_out noi: di as txt _col(4) "Pr(Z > z1) " `p1_display' _col(33) as txt "Pr(Z > z2) " `p2_display'
  }

  local holdg1 = `g1'
  local g1 = `g1'
  local g2 = `g2'

  local valulab : value label `by'
  if `"`valulab'"'!=`""' {
    local g1 : label `valulab' `g1'
    local g2 : label `valulab' `g2'
  }

  local by "`origby'"
  if ("`relevance'"!="") {
    di as txt _newline "Two-sample rank-sum relevance test" _newline
    di as txt "Two-sample test for 0th-order stochastic dominance"  _newline
    di in smcl as txt %12s abbrev(`"`by'"',12) " {c |}" _col(21) "Obs" _col(28) "Rank sum" _col(40) "Expected"
    di in smcl as txt "{hline 13}{c +}{hline 33}"
    _tostranksumditablin `"`g1'"' `n1' `W' `EW'
    _tostranksumditablin `"`g2'"' `n2' `n'*(`n'+1)/2-`W' `n2'*(`n'+1)/2
    di in smcl as txt "{hline 13}{c +}{hline 33}"
    _tostranksumditablin combined `n' `n'*(`n'+1)/2 `n'*(`n'+1)/2
    di in smcl as txt _newline `"unadjusted variance"' _col(22) as res `vfmt' `unv'
    di as txt `"adjustment for ties"' _col(22) as res `vfmt' `v'-`unv'
    di as txt _col(22) "{hline 10}"
    di as txt `"adjusted variance"' _col(22) as res `vfmt' `v' _newline
    local xab = abbrev("`x'",8)
    local byab = abbrev("`by'",8)
    di as txt "Ho: P(`xab'|`byab'=`g1' > `xab'|`byab'=`g2') = 0.5"
    di as txt "Ha: P(`xab'|`byab'=`g1' > `xab'|`byab'=`g2') ≠ 0.5" _newline
    `continuity_notice'
    di as txt _col(14) `"z ="' as res %7.3f `z_pos'
    di as txt _col(2) "Pr(|Z| > |z|) " `p_pos_display' _newline
    }
  di as txt _newline "Two-sample rank-sum test for 0th-order stochastic equivalence" _newline
  di in smcl as txt %12s abbrev(`"`by'"',12) " {c |}" _col(21) "Obs" _col(28) "Rank sum" _col(40) "Expected"
  di in smcl as txt "{hline 13}{c +}{hline 33}"
  _tostranksumditablin `"`g1'"' `n1' `W' `EW'
  _tostranksumditablin `"`g2'"' `n2' `n'*(`n'+1)/2-`W' `n2'*(`n'+1)/2
  di in smcl as txt "{hline 13}{c +}{hline 33}"
  _tostranksumditablin combined `n' `n'*(`n'+1)/2 `n'*(`n'+1)/2
  local xab = abbrev("`x'",8)
  local byab = abbrev("`by'",8)
  di in smcl as txt _newline `"unadjusted variance"' _col(22) as res `vfmt' `unv'
  di as txt `"adjustment for ties"' _col(22) as res `vfmt' `v'-`unv'
  di as txt _col(22) "{hline 10}"
  di as txt `"adjusted variance"' _col(22) as res `vfmt' `v' _newline
  if (lower("`eqvtype'") == "delta") {
    if (`upper' == `lower') {
      noi: di as txt _col(5) "`delta_inflate'`delta' = " as res %-6.4f `lower' as txt _col(22) "`delta' expressed in units of summed ranks (W)"
      }
    if (`upper' != `lower') {
      noi: di as txt _col(8) "`delta_lower' = " as res %-6.4f -1*`lower' as txt _col(22) "`delta_lower' expressed in units of summed ranks (W)"
      noi: di as txt _col(8) "`delta_upper' =  " as res %-6.4f `upper' as txt _col(22) "`delta_upper' expressed in units of summed ranks (W)"
      }
    local criticalvalue = `se'*invnormal(1-`alpha')
    if (`upper' == `lower' & `lower' <= `criticalvalue') {
      noi: di _newline as res "Impossible to reject any Ho if `delta' `lessthan' z-crit`times's.e. ( " %-6.4g `criticalvalue' "). See{help tostranksum##mineqvlevel: help tostranksum}."
      }
    if (`upper' != `lower' & `lower' <= `criticalvalue') {
      noi: di _newline as res "Impossible to reject any Ho if |`delta_lower'| `lessthan' z-crit`times's.e. ( " %-6.4g `criticalvalue' "). See{help tostranksum##mineqvlevel: help tostranksum}."
      }
    if (`upper' != `lower' & `upper' <= `criticalvalue') {
      noi: di _newline as res "Impossible to reject any Ho if `delta_upper' `lessthan' z-crit`times's.e. ( " %-6.4g `criticalvalue' "). See{help tostranksum##mineqvlevel: help tostranksum}."
      }
    if (`upper' == `lower') {
      noi: di _newline as txt "Ho: |W-E(W)| `greaterthan' `delta':" _newline
      `continuity_notice'
      `neg_z_stats_out'
      noi: di as txt _col(4) "Ho1: `delta'-[W-E(W)] `lessthan' 0" _col(33) "Ho2: [W-E(W)]+`delta' `lessthan' 0"
      noi: di as txt _col(4) "Ha1: `delta'-[W-E(W)] > 0"  _col(33) "Ha2: [W-E(W)]+`delta' > 0"
      `neg_p_vals_out'
      }
    if (`upper' != `lower') {
      noi: di _newline as txt "Ho: [W-E(W)] `lessthan' `delta_lower', or [W-E(W)] `greaterthan' `delta_upper':" _newline
      `continuity_notice'
      `neg_z_stats_out'
      noi: di as txt _col(4) "Ho1: `delta_upper'-[W-E(W)] `lessthan' 0" _col(33) "Ho2: [W-E(W)]-`delta_lower' `lessthan' 0"
      noi: di as txt _col(4) "Ha1: `delta_upper'-[W-E(W)] > 0" _col(33) "Ha2: [W-E(W)]-`delta_lower' > 0"
      `neg_p_vals_out'
      }
    }
   if lower("`eqvtype'") == "epsilon" {
    if (`upper' == `lower') {
      noi: di as txt _col(5) "`epsilon_inflate'`epsilon' = " as res %-6.4f `lower' as txt _col(22) "`epsilon' expressed in units of the z distribution"
      }
    if (`upper' != `lower') {
      noi: di as txt _col(8) "`epsilon_lower' = " as res %-6.4f -1*`lower' as txt _col(22) "`epsilon_lower' expressed in units of the z distribution"
      noi: di as txt _col(8) "`epsilon_upper' =  " as res %-6.4f `upper' as txt _col(22) "`epsilon_upper'" as txt " expressed in units of the z distribution"
      }
    local criticalvalue = invnormal(1-`alpha')
    if (`upper' == `lower' & `lower' <= `criticalvalue') {
      noi: di _newline as res "Impossible to reject any Ho if `epsilon' `lessthan' t-crit (" %-5.3f `criticalvalue' "). See{help tostranksum##mineqvlevel: help tostranksum}." _newline
      }
    if (`upper' != `lower' & `lower' <= `criticalvalue') {
      noi: di _newline as res "Impossible to reject any Ho if |`epsilon_lower'| `lessthan' t-crit (" %-5.3f `criticalvalue' "). See{help tostranksum##mineqvlevel: help tostranksum}." _newline
      }
    if (`upper' != `lower' & `upper' <= `criticalvalue') {
      noi: di _newline as res "Impossible to reject any Ho if `epsilon_upper' `lessthan' t-crit (" %-5.3f `criticalvalue' "). See{help tostranksum##mineqvlevel: help tostranksum}." _newline
      }
    if (`upper' == `lower') {
      noi: di _newline as txt "Ho: |Z| `greaterthan' `epsilon':" _newline
      `continuity_notice'
      `neg_z_stats_out'
      noi: di as txt _col(4) "Ho1: `epsilon'-Z `lessthan' 0" _col(33) "Ho2: Z+`epsilon' `lessthan' 0"
      noi: di as txt _col(4) "Ha1: `epsilon'-Z > 0" _col(33) "Ha2: Z+`epsilon' < 0"
      `neg_p_vals_out'
      }
    if (`upper' != `lower') {
      noi: di _newline as txt "Ho: Z `lessthan' `epsilon_lower', or Z `greaterthan' `epsilon_upper':" _newline
      `continuity_notice'
      `neg_z_stats_out'
      noi: di as txt _col(4) "Ho1: `epsilon_upper'-Z `lessthan' 0" _col(33) "Ho2: Z-`epsilon_lower' `lessthan' 0"
      noi: di as txt _col(4) "Ha1: `epsilon_upper'-Z > 0" _col(33) "Ha2: Z-`epsilon_lower' > 0"
      `neg_p_vals_out'
      }
    }

* Output combined tests results if relevance test is requested
 if "`relevance'" != "" {
   * Format alpha to remove trailing zeros
   if (mod(`alpha'*1000, 1) == 0.0) {
     local alpha: di %6.3f `alpha'
     }
   if (mod(`alpha'*100, 1) == 0.0) {
     local alpha: di %5.2f `alpha'
     }
   if (mod(`alpha'*10, 1) == 0.0) {
     local alpha: di %4.1f `alpha'
     }
   if (mod(`alpha', 1) == 0.0) {
     local alpha: di %4.0f `alpha'
     }

   * Format Delta or epsilon to remove trailing zeros
   if (mod(`lower'*1000, 1) == 0.0) {
     local lower: di %6.3f `lower'
     }
   if (mod(`lower'*100, 1) == 0.0) {
     local lower: di %5.2f `lower'
     }
   if (mod(`lower'*10, 1) == 0.0) {
     local lower: di %4.1f `lower'
     }
   if (mod(`lower', 1) == 0.0) {
     local lower: di %4.0f `lower'
     }
   if (`upper' != `lower') {
     if (mod(`upper'*1000, 1) == 0.0) {
       local upper: di %6.3f `upper'
       }
     if (mod(`upper'*100, 1) == 0.0) {
       local upper: di %5.2f `upper'
       }
     if (mod(`upper'*10, 1) == 0.0) {
       local upper: di %4.1f `upper'
       }
     if (mod(`upper', 1) == 0.0) {
       local upper: di %4.0f `upper'
       }
     }
   if (`upper' == `lower') {

     if lower("`eqvtype'") == "delta" {
       if (`lower' < 1) {
         noi: di _newline _newline as txt "Relevance test conclusion for `alpha_out' = " as res "0" `alpha' as txt ", and `delta' = " as res "0" `lower' as txt ":"
         }
        else {
         noi: di _newline _newline as txt "Relevance test conclusion for `alpha_out' = " as res "0" `alpha' as txt ", and `delta' = " as res `lower' as txt ":"
         }
       }
     if lower("`eqvtype'") == "epsilon" {
       noi: di _newline _newline as txt "Relevance test conclusion for `alpha_out' = " as res "0" `alpha' as txt ", and `epsilon' = " as res `lower' as txt ":"
       }
     }
   if (`upper' != `lower') {

     if lower("`eqvtype'") == "delta" {
       if (`lower' < 1 & `upper' < 1) {
         noi: di _newline _newline as txt "Relevance test conclusion for `alpha_out' = " as res "0" `alpha' as txt ", `delta_lower' = " as res "-0" `lower' as txt ", and `delta_upper' = " as res "0" %-6.4f `upper' as txt ":"
         }
       if (`lower' >= 1 & `upper' < 1) {
         noi: di _newline _newline as txt "Relevance test conclusion for `alpha_out' = " as res "0" `alpha' as txt ", `delta_lower' = " as res -1*`lower' as txt  ", and `delta_upper' = " as res "0" %-6.4f `upper' as txt ":"
         }
       if (`lower' < 1 & `upper' >= 1) {
         noi: di _newline _newline as txt "Relevance test conclusion for `alpha_out' = " as res "0" `alpha' as txt ", `delta_lower' = " as res "-0" `lower' as txt  ", and `delta_upper' = " as res %-6.4f `upper' as txt ":"
         }
       if (`lower' >= 1 & `upper' >= 1) {
         noi: di _newline _newline as txt "Relevance test conclusion for `alpha_out' = " as res "0" `alpha' as txt ", `delta_lower' = " as res -1*`lower' as txt  ", and `delta_upper' = " as res %-6.4f `upper' as txt ":"
         }
       }
     if lower("`eqvtype'") == "epsilon" {
       noi: di _newline _newline as txt "Relevance test conclusion for `alpha_out' = " as res "0" `alpha' as txt ", `epsilon_lower' = " as res -1*`lower' as txt  ", and `epsilon_upper' = " as res `upper' as txt ":"
       }
     }
   noi: di as txt "  Ho test for 0th-order stochastic dominance:   " as res "`PositivistConclusion'"
   noi: di as txt "  Ho test for 0th-order stochastic equivalence: " as res "`NegativistConclusion'"
   if "`PositivistConclusion'" == "Reject" & "`NegativistConclusion'" == "Reject" {
     local RelevanceTestConclusion = "Trivial difference (overpowered test)"
     }
   if "`PositivistConclusion'" == "Reject" & "`NegativistConclusion'" == "Fail to reject" {
     local RelevanceTestConclusion = "Relevant difference"
     }
   if "`PositivistConclusion'" == "Fail to reject" & "`NegativistConclusion'" == "Reject" {
     local RelevanceTestConclusion = "Equivalence"
     }
   if "`PositivistConclusion'" == "Fail to reject" & "`NegativistConclusion'" == "Fail to reject" {
     local RelevanceTestConclusion = "Indeterminate (underpowered test)"
     }
   noi: di _newline as txt "Conclusion from combined tests: " as res "`RelevanceTestConclusion'"
   }

*******************************************************************************
* Program end. Close up shop and return things.                               *
*******************************************************************************
  if ("`relevance'" != "") {
    return local relevance = "`RelevanceTestConclusion'"
    }
  if (`upper' != `lower') {
    if "`eqvtype'" == "delta" {
      return scalar Du   = `upper'
      return scalar Dl   = `lower'
      }
    if "`eqvtype'" == "epsilon" {
      return scalar eu   = `upper'
      return scalar el   = `lower'
      }
    }

  if (`upper' == `lower') {
    if "`eqvtype'" == "delta" {
      return scalar Delta   = `eqvlevel'
      }
    if "`eqvtype'" == "epsilon" {
      return scalar epsilon   = `eqvlevel'
      }
    }

  return scalar  sum_exp  = `EW'
  return scalar  sum_obs  = `W'
  return scalar  group1   = `holdg1'
  return scalar  Var_a    = `v'
  return scalar  p2       = `p2'
  return scalar  p1       = `p1'
  return scalar  z2       = `z2'
  return scalar  z1       = `z1'
  return scalar  N_2      = `n2'
  return scalar  N_1      = `n1'

  end


program define _tostranksumditablin
  if length(`"`1'"') > 12 {
    local 1 = substr(`"`1'"',1,12)
    }
  di in smcl as txt %12s `"`1'"' " {c |}" as res _col(17) %7.0g `2' _col(26) %10.0g `3' _col(38) %10.0g `4'
  end