*! version 3.1.5  12oct2024 by alexis.dinno@pdx.edu
*! perform two one-sided tests for  the distribution of paired or matched data
*! being equivalent to one that is symmetrical & centered on zero

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

program define tostsignrank, rclass byable(recall)

  if int(_caller())<8 {
    di in r "tostranksum- 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)
    }

  /* turn "==" into "=" if needed before calling -syntax- */
  gettoken vn rest : 0, parse(" =")
  gettoken eq rest : rest, parse(" =")
  if "`eq'" == "==" {
    local 0 `vn' = `rest'
    }

  syntax varname [=/exp] [if] [in] [, EQVType(string) EQVLevel(real 1) /*
*/      UPPEReqvlevel(real 0) CContinuity Alpha(real 0.05) RELevance]



  * 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
    }


  tempname `total_positive' `total_negative' v unv z adj0
  tempvar touse diff ranks t

  quietly {
    mark `touse' `if' `in'
    gen double `diff' = `varlist'-(`exp') if `touse'
    markout `touse' `diff'
    egen double `ranks' = rank(abs(`diff')) if `touse'

    /* We do want to OMIT the ranks corresponding to `diff'==0 in the sums.  */
    gen double `t' = sum(cond(`diff'>0,`ranks',0))
    local total_positive = `t'[_N]
    replace `t' = sum(cond(`diff'<0,`ranks',0))
    local total_negative = `t'[_N]
    replace `t' = sum(cond(`diff'~=0,`ranks'*`ranks',0))
    local v = `t'[_N]/4
    local se = 2*sqrt(`v')
    if (lower("`ccontinuity'") == "") {
      local cc = 0
      local continuity_notice `""'
      }
    if (lower("`ccontinuity'") != "") {
      local cc = 0.5
      local continuity_notice `"noi: di as txt "Using continuity correction" _newline"'
      }
    local PositivistConclusion = "Reject"
    local NegativistConclusion = "Reject"
    local z_pos = ((sign(`total_positive'-`total_negative')*(abs(`total_positive'-`total_negative') - `cc')))/`se'
    local z_pos_display = `"%-6.3f `z_pos'"'
    local p_pos = 2*(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') {
      local PositivistConclusion = "Fail to reject"
      }
    if lower("`eqvtype'") == "delta" {
        local z1 = (`upper' - (sign(`total_positive' - `total_negative') * (abs(`total_positive' - `total_negative') - `cc'))) / `se'
        local z2 = ((sign(`total_positive' - `total_negative') * (abs(`total_positive' - `total_negative') - `cc')) + `lower') / `se'
        }
    if lower("`eqvtype'") == "epsilon" {
        local z1 = `upper' - ( (sign(`total_positive' - `total_negative') * (abs(`total_positive' - `total_negative') - `cc'))/`se' )
        local z2 = ( (sign(`total_positive' - `total_negative') * (abs(`total_positive' - `total_negative') - `cc'))/`se' ) + `lower'
        }
    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"
      }
    count if `touse'
    local n = r(N)
    local unv = `n'*(`n'+1)*(2*`n'+1)/24
    if `unv' < 1e7 {
      local vfmt `"%10.2f"'
      }
     else {
      local vfmt `"%10.0g"'
      }
    count if `diff' == 0 & `touse'
    local number_ties = r(N)
    local adj0 = -`number_ties'*(`number_ties'+1)*(2*`number_ties'+1)/24
    count if `diff' > 0 & `touse'
    local number_positive = r(N)
    local number_negative = `n' - `number_positive' - `number_ties'
    local table_header noi: di as txt _col(9) "Sign {c |}" _col(21) "Obs" _col(27) "Sum ranks" _col(40) "Expected"
    if "`relevance'" != "" {
      noi: di as txt _newline "Relevance Wilcoxon signed-rank test"
      noi: di _newline as txt `"Signed-rank test for the distribution of paired or matched data being"' _newline `"different from one that is symmetrical & centered on zero"' _newline
      `table_header'
      noi: di in smcl as txt "{hline 13}{c +}{hline 33}"
      noi: _tostsignrankditablein "positive" `number_positive' `total_positive' (`total_positive'+`total_negative')/2
      noi: _tostsignrankditablein "negative" `number_negative' `total_negative' (`total_positive'+`total_negative')/2
      noi: _tostsignrankditablein "zero"     `number_ties' `number_ties'*(`number_ties'+1)/2 `number_ties'*(`number_ties'+1)/2
      noi: di in smcl as txt "{hline 13}{c +}{hline 33}"
      noi: _tostsignrankditablein all `n' `n'*(`n'+1)/2 `n'*(`n'+1)/2
      noi: di in smcl as txt _newline `"Unadjusted variance"' _col(22) as res `vfmt' `unv'
      noi: di as txt `"Adjustment for ties"' _col(22) as res `vfmt' `v'-`unv'-`adj0'
      noi: di as txt `"Adjustment for zeros"' _col(22) as res `vfmt' `adj0'
      noi: di as txt _col(22) "{hline 10}"
      noi: di as txt `"Adjusted variance"' _col(22) as res `vfmt' `v' _newline
      noi: di as txt "Ho: `varlist' - `exp' has a symmetric distribution centered on zero"
      noi: di as txt "Ha: `varlist' - `exp' has an asymmetric distribution," _newline _col(5) "a distribution not centered on zero, or" _newline _col(5) "has a distribution which is asymmetric and not centered on zero" _newline
      `continuity_notice'
      noi: di as txt _col(12) "z = " as res `z_pos_display'
      noi: di as txt _col(2) "Pr(|Z|>|z|) " as res `p_pos_display' _newline
      }
    }

  di _newline as txt "Signed-rank test for the distribution of paired or matched data being"
  di as txt "equivalent to one that is symmetrical & centered on zero" _newline
  `table_header'
  di in smcl as txt "{hline 13}{c +}{hline 33}"
  _tostsignrankditablein positive `number_positive' `total_positive' (`total_positive'+`total_negative')/2
  _tostsignrankditablein negative `number_negative' `total_negative' (`total_positive'+`total_negative')/2
  _tostsignrankditablein zero     `number_ties' `number_ties'*(`number_ties'+1)/2 `number_ties'*(`number_ties'+1)/2
  di in smcl as txt "{hline 13}{c +}{hline 33}"
  _tostsignrankditablein all `n' `n'*(`n'+1)/2 `n'*(`n'+1)/2

  if `unv' < 1e7 {
    local vfmt `"%10.2f"'
    }
   else
    local vfmt `"%10.0g"'

  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'-`adj0'
  di as txt `"Adjustment for zeros"' _col(22) as res `vfmt' `adj0'
  di as txt _col(22) "{hline 10}"
  di as txt `"Adjusted variance"' _col(22) as res `vfmt' `v' _newline
  local z_stats_out `"noi: di as txt _col(8) "z1 = " as res %-8.4g `z1' as txt _col(38) "z2 = " as res %-8.4g `z2' _newline"'
  local p_vals_out `"noi: di as txt _col(4) "Pr(Z > t1) " `p1_display' _col(33) as txt "Pr(Z > t2) " `p2_display'"'
  if (lower("`eqvtype'") == "delta") {
    if (`upper' == `lower') {
      noi: di as txt _col(5) "`delta_inflate'`delta' = " as res %-8.0f `lower' as txt _col(22) "`delta' expressed in units of signed ranks (T)"
      }
    if (`upper' != `lower') {
      noi: di as txt _col(8) "`delta_lower' = " as res %-7.0g -1*`lower' as txt _col(21) "`delta_lower' expressed in units of signed ranks (T)"
      noi: di as txt _col(8) "`delta_upper' =  " as res %-6.0g `upper' as txt _col(21) "`delta_upper' expressed in units of signed ranks (T)"
      }
    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 tostsignrank##mineqvlevel: help tostsignrank}."
      }
    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 tostsignrank##mineqvlevel: help tostsignrank}."
      }
    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 tostsignrank##mineqvlevel: help tostsignrank}."
      }
    if (`upper' == `lower') {
      noi: di _newline as txt "Ho: |T-E(T)| `greaterthan' `delta':" _newline
      `continuity_notice'
      `z_stats_out'
      noi: di as txt _col(4) "Ho1: `delta'-[T-E(T)] `lessthan' 0`'" _col(33) "Ho2: [T-E(T)]+`delta' `lessthan' 0"
      noi: di as txt _col(4) "Ha1: `delta'-[T-E(T)] > 0"  _col(33) "Ha2: [T-E(T)]+`delta' > 0"
      `p_vals_out'
      }
    if (`upper' != `lower') {
      noi: di _newline as txt "Ho: [T-E(T)] `lessthan' `delta_lower', or [T-E(T)] `greaterthan' `delta_upper':" _newline
      `continuity_notice'
      `z_stats_out'
      noi: di as txt _col(4) "Ho1: `delta_upper'-[T-E(T)] `lessthan' 0" _col(33) "Ho2: [T-E(T)]-`delta_lower' `lessthan' 0"
      noi: di as txt _col(4) "Ha1: `delta_upper'-[T-E(T)] > 0" _col(33) "Ha2: [T-E(T)]-`delta_lower' > 0"
      `p_vals_out'
      }
    }
   if lower("`eqvtype'") == "epsilon" {
    if (`upper' == `lower') {
      noi: di as txt _col(5) "`epsilon_inflate'`epsilon' = " as res %-5.3f `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.3f -1*`lower' as txt _col(21) "`epsilon_lower' expressed in units of the z distribution"
      noi: di as txt _col(8) "`epsilon_upper' =  " as res %-5.3f `upper'  as txt _col(21) "`epsilon_upper' 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' z-crit (" %-5.3f `criticalvalue' "). See{help tostsignrank##mineqvlevel: help tostsignrank}."
      }
    if (`upper' != `lower' & `lower' <= `criticalvalue') {
      noi: di _newline as res "Impossible to reject any Ho if |`epsilon_lower'| `lessthan' z-crit (" %-5.3f `criticalvalue' "). See{help tostsignrank##mineqvlevel: help tostsignrank}."
      }
    if (`upper' != `lower' & `upper' <= `criticalvalue') {
      noi: di _newline as res "Impossible to reject any Ho if `epsilon_upper' `lessthan' z-crit (" %-5.3f `criticalvalue' "). See{help tostsignrank##mineqvlevel: help tostsignrank}."
      }
    if (`upper' == `lower') {
      noi: di _newline as txt "Ho: |Z| `greaterthan' `epsilon':" _newline
      `continuity_notice'
      `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"
      `p_vals_out'
      }
    if (`upper' != `lower') {
      noi: di _newline as txt "Ho: Z `lessthan' `epsilon_lower', or Z `greaterthan' `epsilon_upper':" _newline
      `continuity_notice'
      `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"
      `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" `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" `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 `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 `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 difference:  " as res "`PositivistConclusion'"
   noi: di as txt "  Ho test for 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'
      }
    }

  ret scalar p2 = `p2'
  ret scalar p1 = `p1'
  ret scalar z2 = `z2'
  ret scalar z1 = `z1'
  ret scalar Var_a = `v'
  ret scalar sum_negative = `total_negative'
  ret scalar sum_pos = `total_positive'
  ret scalar N_tie = `number_ties'
  ret scalar N_pos = `number_positive'
  ret scalar N_negative = `number_negative'
  end


program define _tostsignrankditablein
  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