*! version 3.1.5  12oct2024 by alexis.dinno@pdx.edu
*! perform two one-sided tests for association in paired data

* Syntax:  tostmcci #a #b #c #d [if exp] [in range] [fweight exp], [
*          eqvlevel(#) uppereqvlevel(#) yates edwards alpha(#) relevance]

program define tostmcci, rclass byable(recall)

  if int(_caller())<8 {
    di in r "tostmcci- 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 theta           = "theta"
    local theta_inflate   = ""
    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 chisq           = "chi-sq"
    local chisq_inflate   = ""
    local sq              = "^2"
    local sq_inflate      = ""
    local lessthan        = "<="
    local greaterthan     = ">="
    local notequal        = "!="
    local rquote          = `"'"'
    local times           = "*"
    }
   else {
    version 14.0, missing
    * Create unicode output constants
    local alpha_out   = uchar(0945)
    local theta           = uchar(0952)
    local theta_inflate   = "    "
    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 chisq           = uchar(0967)+uchar(0178)
    local chisq_inflate   = "    "
    local sq              = uchar(0178)
    local sq_inflate      = " "
    local lessthan        = uchar(8804)
    local greaterthan     = uchar(8805)
    local notequal        = uchar(8800)
    local rquote          = uchar(8217)
    local times           = uchar(0215)
    }

  gettoken a 0 : 0, parse(" ,")
  gettoken b 0 : 0, parse(" ,")
  gettoken c 0 : 0, parse(" ,")
  gettoken d 0 : 0, parse(" ,")

  confirm integer number `a'
  confirm integer number `b'
  confirm integer number `c'
  confirm integer number `d'

  if `a'<0 | `b'<0 | `c'<0 | `d'<0 {
    di in red "negative numbers invalid"
    exit 498
  }

  syntax  [, EQVType(string) EQVLevel(real 1) UPPEReqvlevel(real 0) YAtes /*
*/        EDwards Alpha(real 0.05) RELevance]

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

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

  * Validate eqvlevel
  if (lower("`eqvtype'") == "delta") & (`eqvlevel' == 1 & `uppereqvlevel'==0) {
    local eqvlevel = 0.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 continuity correction option
  if ("`yates'" != "" & "`edwards'" != "") {
    noi di as err "continuity correction options must be either yates or edwards, but not both"
    exit 198
    }

  * 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 continuity = 0
  local correction_message di "" _continue
  if ("`yates'" != "") {
    local continuity = 0.5
    local correction_message di as txt _newline "Using the Yates continuity correction"
    }
  if ("`edwards'" != "") {
    local continuity = 1
    local correction_message di as txt _newline "Using the Edwards continuity correction"
    }
  local chisq_pos = ((abs(`b' - `c') - `continuity')^2)/(`b' + `c')
  local p_pos = chi2tail(1,`chisq_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 ("`relevance'" != "") {
    di as txt _newline "Relevance test for paired binary data"
    local level = 100*(1 - `alpha')
    noi: di as txt _newline "McNemar`rquote's test for difference in paired binary data"
    di as txt "{hline 17}{c TT}{hline 24}{c TT}{hline 11}"
    _crc4fld `a' `b' `c' `d' /*
      */ Controls Cases Exposed Unexposed Exposed Unexposed
    di as txt "{hline 17}{c BT}{hline 24}{c BT}{hline 11}"
    noi: di as txt "McNemar`rquote's `chisq'(1) = " as res %8.4f _col(22) `chisq_pos' as txt _col(32) "`chisq_inflate'Pr(X > `chisq') " `p_pos_display'
    noi: di as txt "Ho: Pr(Exposed|Case) = Pr(Exposed|Control)"
    noi: di as txt "Ha: Pr(Exposed|Case) `notequal' Pr(Exposed|Control)"
    noi: `correction_message'
    noi: di
    }
  noi: di as txt _newline "Liu, et al. test for equivalence in paired binary data"
  noi: di as txt "{hline 17}{c TT}{hline 24}{c TT}{hline 11}"

  _crc4fld `a' `b' `c' `d' /*
    */ Controls Cases Exposed Unexposed Exposed Unexposed
  noi: di as txt "{hline 17}{c BT}{hline 24}{c BT}{hline 11}"
  local n = `a'+`b'+`c'+`d'
  local den = `b' + `c'
  local low = min(`b',`c')
  local diff = `b' - `c'
  local se = sqrt((`b'+`c') - (`n'*(((`b'/`n')-(`c'/`n'))^2)))
  if (lower("`eqvtype'")=="delta") {
    local z1 = ( (`n'*`upper') - ((`diff') - `continuity') )/`se'
    local z2 = ( ((`diff') + `continuity') + (`n'*`lower') )/`se'
    }
  if (lower("`eqvtype'")=="epsilon") {
    local z1 = `upper' - sign(`diff')*(( abs(`diff') - `continuity' )/`se')
    local z2 = sign(`diff')*(( abs(`diff') + `continuity' )/`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"
    }
  noi: di as txt "      diff = " as res `b' as txt " - " as res `c' as txt " = " as res `diff'
  noi: di as txt "`theta_inflate's.e. `theta' = " as res %-8.4f `se'
  local 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 p_vals_out noi: di as txt _col(4) "Pr(Z > z1) " `p1_display' _col(33) as txt "Pr(Z > z2) " `p2_display'
  if lower("`eqvtype'") == "delta" {
    if (`upper' == `lower') {
      noi: di as txt "         `delta' = " as res %-8.4f `lower' as txt _col(23) "`delta' expressed in units of probability"
      }
    if (`upper' != `lower') {
      noi: di as txt "        `delta_lower' = " as res %-8.4f -1*`lower' as txt _col(23) "`delta_lower' expressed in units of probability"
      noi: di as txt "        `delta_upper' =  " as res %-8.4f `upper' as txt _col(23) "`delta_upper' expressed in units of probability"
      }
    local criticalvalue = (`se'*invnormal(1-`alpha'))/`n'
    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 tostmcc##mineqvlevel: help tostmcc}."
      }
    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 tostmcc##mineqvlevel: help tostmcc}."
      }
    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 tostmcc##mineqvlevel: help tostmcc}."
      }
    if (`upper' == `lower') {
      noi: di _newline as txt "Ho: |`theta'| `greaterthan' `delta':" _newline
      `correction_message'
      `z_stats_out'
      noi: di as txt _col(4) "Ho1: `delta'-`theta' `lessthan' 0" _col(33) "Ho2: `theta'+`delta' `lessthan' 0"
      noi: di as txt _col(4) "Ha1: `delta'-`theta' > 0"  _col(33) "Ha2: `theta'+`delta' > 0"
      `p_vals_out'
      }
    if (`upper' != `lower') {
      noi: di _newline as txt "Ho: `theta' `lessthan' `delta_lower', or `theta' `greaterthan' `delta_upper':" _newline
      `correction_message'
      `z_stats_out'
      noi: di as txt _col(4) "Ho1: `delta_upper'-`theta' `lessthan' 0" _col(33) "Ho2: `theta'-`delta_lower' `lessthan' 0"
      noi: di as txt _col(4) "Ha1: `delta_upper'-`theta' > 0"  _col(33) "Ha2: `theta'-`delta_lower' > 0"
      `p_vals_out'
      }
    }
  if lower("`eqvtype'") == "epsilon" {
    if (`upper' == `lower') {
      noi: di as txt _col(4) "`epsilon_inflate'`epsilon' = " as res %-8.4f `lower' as txt _col(23)  "`epsilon' expressed in units of the z distribution"
      }
    if (`upper' != `lower') {
      noi: di as txt _col(8) "`epsilon_lower' = " as res %-8.4f -1*`lower' as txt _col(23) "`epsilon_lower' expressed in units of the z distribution"
      noi: di as txt _col(8) "`epsilon_upper' =  " as res %-8.4f `upper' as txt _col(23) "`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 tostmcc##mineqvlevel: help tostmcc}."
      }
    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 tostmcc##mineqvlevel: help tostmcc}."
      }
    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 tostmcc##mineqvlevel: help tostmcc}."
      }
    if (`upper' == `lower') {
      noi: di _newline as txt "Ho: |Z| `greaterthan' `epsilon':" _newline
      `correction_message'
      `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
      `correction_message'
      `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'
      }
    }
  return scalar  D_f  = (`b'/`n') - (`c'/`n')
  return scalar  p2   = `p2'
  return scalar  p1   = `p1'
  return scalar  z2   = `z2'
  return scalar  z1   = `z1'
  end