*! version 3.1.5  12oct2024 by alexis.dinno@pdx.edu
*! perform two one-sided tests for proportion equivalence

* Syntax:  tostpr varname [=exp] [if exp] [in range] [, eqvtype(type)
*          eqvlevel(#) uppereqvlevel(#) by(groupname) yates ha alpha(#)
*          relevance]

program define tostpr, rclass byable(recall)

  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 theta           = "theta"
    local theta_inflate   = "{space 4}"
    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 notequal        = "!="
    local caret           = "^"
    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 lessthan        = uchar(8804)
    local greaterthan     = uchar(8805)
    local notequal        = uchar(8800)
    local caret           = uchar(0710)
    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) BY(varname) YAtes ha Alpha(real 0.05) /*
*/       RELevance]

  quietly {

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

  if lower("`eqvtype'") == "delta" & (`eqvlevel' >= 1 | `uppereqvlevel' >= 1) {
    noi: di as err "option eqvlevel() incorrectly specified" _newline "you are likely to find all proportions equivalent within an interval of plus or minus 1 or more"
    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')
    }

  * Invalidate options specified with one-sample test
  capture confirm number `exp'
  if _rc == 0 {
    if "`by'" != "" {
      noi di as err "may not combine = and option by()"
      exit 198
      }

    if ("`ha'" != "") {
      noi di _newline as res "Hauck-Anderson continuity correction not relevant in one-sample tests"
      local ha = ""
      }
    }

  * Validate continuity correction option
  if ("`yates'" != "" & "`ha'" != "") {
    noi di as err "continuity correction options must be either yates or ha, but not both"
    exit 198
    }

  local continuity = 0

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

*******************************************************************************
* The business starts here                                                    *
*******************************************************************************

**********
* conduct the positivist z test for proportion difference
  local PositivistConclusion = "Reject"
  local NegativistConclusion = "Reject"
  if "`exp'" != "" & "`by'" == "" {
    local eexp = "= `exp'"
    }
  local level: di %3.2f (100*(1-`alpha'))
  * No continuity correction
  if ("`yates'"=="") {
    local correction = 0
    local correction_message di "" _continue
    }
  quietly: prtest `varlist' `eexp' `if' `in', by("`by'") level(`level')
  capture confirm number `exp'

**********
* one-sample test for proportion equivalence
  if _rc == 0 {
    local n1 = r(N_1)
    local m1 = r(P_1)
    local se = sqrt(`m1'*(1-`m1')/(`n1'))
    local se_p = sqrt(`exp'*(1-`exp')/`n1')
    if (1-normal(abs( (`m1' - `exp')/sqrt(`m1'*(1-`m1')/(`n1')) ))) > (`alpha'/2) {
      local PositivistConclusion = "Fail to reject"
      }
    * No continuity correction for one-sample test
    if ("`yates'"=="") {
      local correction = 0
      local correction_message di "" _continue
      }
    * Yates continuity correction for one-sample test
    if ("`yates'"!="") {
      local correction = 0.5*(1/`n1')
      local correction_message di as txt "Using the Yates continuity correction" _newline
      }
    local z_pos = sign(`m1'-`exp') * (abs(`m1' - `exp') - `correction') / `se_p'
    if "`relevance'" != "" {
      local name1 = trim(substr(trim(`"`varlist'"'),1,12))
      local c1 = 49 - length("`name1'")
      noi: di as txt _newline "Relevance test of sample proportions" _newline
      noi: di as txt `"One-sample z test for proportion difference"' _newline /*
        */ _col(`c1') "`name1'" _col(49) /*
        */ `": Number of obs = "' as res %14.0g `n1'
      noi: _tostprtest1 `varlist' `n1' `m1' `se' `alpha'
      noi: di in smcl as txt "{hline 13}{c BT}{hline 66}"
      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"
        }
      noi: di as txt _col(7) "`theta_inflate'`theta'`caret' = prop(`name1'=1)" _col(70) "z = " as res %7.0g `z_pos'
      noi: di as txt _col(4) "`theta_inflate'Ho: `theta' - " as res "`exp'" as txt " = 0" _newline
      noi: `correction_message'
      noi: di as txt _col(37) "Ha: `theta' - " as res "`exp' `notequal' 0"
      noi: di as txt _col(31) "Pr(|Z| > |z|) " `p_pos_display' _newline
      }
    if lower("`eqvtype'") == "delta" {
      local z1 = (`upper' - (sign(`m1' - `exp') * ((`m1' - `exp') - `continuity'))) / `se_p'
      local z2 = ((sign(`m1' - `exp') * ((`m1' - `exp') - `continuity')) + `lower') / `se_p'
      }
    if lower("`eqvtype'") == "epsilon" {
      local z1 = `upper' - `z_pos'
      local z2 = `z_pos' + `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"
      }
    local name1 = trim(substr(trim(`"`varlist'"'),1,12))
    local c1 = 49 - length("`name1'")
    noi: di
    noi: di as txt `"One-sample z test for proportion equivalence"' _newline /*
        */ _col(`c1') "`name1'" _col(49) /*
        */ `": Number of obs = "' as res %14.0g `n1'
     noi: _tostprtest1 `varlist' `n1' `m1' `se' `alpha'
    if lower("`eqvtype'") == "delta" {
      noi: di in smcl as txt "{hline 13}{c +}{hline 66}"
      if (`upper' == `lower') {
        noi: _tostprtest2 "`delta'-`theta'`caret'" `n1' `upper'-(`m1'-`exp') `se_p' `alpha'
        noi: _tostprtest2 "`theta'`caret'+`delta'" `n1' (`m1'-`exp')+`lower' `se_p' `alpha'
        }
      if (`upper' != `lower') {
        noi: _tostprtest2 "`delta_upper'-`theta'`caret'" `n1' `upper'-(`m1'-`exp') `se_p' `alpha'
        noi: _tostprtest2 "`theta'`caret'-`delta_lower'" `n1' (`m1'-`exp')+`lower' `se_p' `alpha'
        }
      }
    noi: di in smcl as txt "{hline 13}{c BT}{hline 66}"
    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 neg_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" {
      noi: di as txt _col(7) "`theta_inflate'`theta'`caret' = prop(" abbrev(`"`varlist'"', 12) "=1) - " as res `exp' as txt " = " as res `m1' - `exp'
      if (`upper' == `lower') {
        noi: di as txt _col(8) "`delta_inflate'`delta' = " as res %-6.3g `lower' as txt _col(24) "`delta' expressed in same units as prop(" abbrev(`"`varlist'"', 12) as txt "=1)"
        }
      if (`upper' != `lower') {
        noi: di as txt _col(11) "`delta_lower' = " as res %-7.3g -1*`lower' as txt _col(23) "`delta_lower' expressed in same units as prop(" abbrev(`"`varlist'"', 12) as txt "=1)"
        noi: di as txt _col(11) "`delta_upper' =  " as res %-6.3g `upper' as txt _col(23) "`delta_upper' expressed in same units as prop(" abbrev(`"`varlist'"', 12) as txt "=1)"
        }
      local criticalvalue = `se_p'*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 tostpr##mineqvlevel: help tostpr}."
        }
      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 tostpr##mineqvlevel: help tostpr}."
        }
      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 tostpr##mineqvlevel: help tostpr}."
        }
      if (`upper' == `lower') {
        noi: di _newline as txt "Ho: |`theta'| `greaterthan' `delta':" _newline
        noi: `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"
        `neg_p_vals_out'
        }
     if (`upper' != `lower') {
        noi: di _newline as txt "Ho: `theta' `lessthan' `delta_lower', or `theta' `greaterthan' `delta_upper':" _newline
        noi: `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"
        `neg_p_vals_out'
        }
      }
    if lower("`eqvtype'") == "epsilon" {
      noi: di as txt _col(7) "`theta_inflate'`theta'`caret' = prop(" abbrev(`"`varlist'"', 12) "=1) - " as res `exp' as txt " = " as res `m1' - `exp'
      if (`upper' == `lower') {
        noi: di as txt _col(6) "`epsilon_inflate'`epsilon' = " as res %-5.3f `lower' as txt _col(24) "`epsilon' expressed in units of the z distribution"
        }
      if (`upper' != `lower') {
        noi: di as txt _col(11) "`epsilon_lower' = " as res %-6.3f -1*`lower' as txt _col(23) "`epsilon_lower' expressed in units of the z distribution"
        noi: di as txt _col(11) "`epsilon_upper' =  " as res %-5.3f `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 tostpr##mineqvlevel: help tostpr}."
        }
      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 tostpr##mineqvlevel: help tostpr}."
        }
      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 tostpr##mineqvlevel: help tostpr}."
        }
      if (`upper' == `lower') {
        noi: di _newline as txt "Ho: |Z| `greaterthan' `epsilon':" _newline
        noi: `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"
        `neg_p_vals_out'
        }
      if (`upper' != `lower') {
        noi: di _newline as txt "Ho: Z `lessthan' `epsilon_lower', or Z `greaterthan' `epsilon_upper':" _newline
        noi: `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"
        `neg_p_vals_out'
        }
      }
    }

**********
* Two-sample z test for proportion difference
  if _rc != 0 {
    local n1 = r(N_1)
    local m1 = r(P_1)
    local se1= sqrt(`m1'*(1-`m1')/(`n1'))
    local n2 = r(N_2)
    local m2 = r(P_2)
    local se2= sqrt(`m2'*(1-`m2')/(`n2'))
    local N  = `n1'+`n2'
    local p  = (`m1'*`n1' +  `m2'*`n2')/(`N')
    local se_p = sqrt(`p' * (1 - `p') * ((1 / `n1') + (1 / `n2')) )
    local correction_message di "" _continue
    * No continuity correction for one-sample test
    if ("`yates'"=="") {
      local correction = 0
      local correction_message di "" _continue
      }
    * Yates' continuity correction
    if ("`yates'" != "") {
      local continuity = 0.5 * ((1 / `n1') + (1 / `n2'))
      local correction_message di as txt "Using the Yates continuity correction" _newline
      }
    * Hauck-Anderson continuity correction
    if ("`ha'" != "") {
      local continuity = 1 / (2 * min(`n1', `n2'))
      local se_p = sqrt( (((`m1') * (1 - `m1'))/(`n1' - 1)) + ((`m2' * (1 - `m2'))/(`n2' - 1)) )
      local correction_message di as txt "Using the Hauck-Anderson continuity correction" _newline
      }
    local z_pos = sign(`m1' - `m2') * (abs(`m1' - `m2') - `continuity') / `se_p'
    if "`by'" != "" {
      sum `by'
      local min = r(min)
      local max = r(max)
      }
    if "`relevance'" != "" {
      noi: di as txt _newline "Relevance test of sample proportions" _newline
      if "`by'" != "" {
        local name1 = "`: label (`by')`min''"
        local name2 = "`: label (`by')`max''"
        local c1 = 48 - length("`by'") - length(`"`name1'"')
        local c2 = 48 - length("`by'") - length(`"`name2'"')
        noi: di "Two-sample z test for proportion difference" _newline /*
          */ _col(`c1') "`by'=`name1'" _col(48) /*
          */ `": Number of obs = "' as res %14.0g `n1'
        noi di as txt _col(`c2') "`by'=`name2'"  _col(48) /*
	  */ `": Number of obs = "' as res %14.0g `n2'
        noi: _tostprtest0 "`by'=`name1'" `n1' `m1' `se1' `alpha'
        noi: _tostprtest2 "`by'=`name2'" `n2' `m2' `se2' `alpha'
        }
       else {
        local name1 = ustrtrim(ustrtrim(substr(ustrtrim(`"`varlist'"'),1,12)))
        local name2 = ustrtrim(ustrtrim(substr(ustrtrim(`"`exp'"'),1,12)))
        local c1 = 49 - length(`"`name1'"')
        local c2 = 49 - length(`"`name2'"')
	noi: di as txt "Two-sample z test for proportion difference" _newline /*
          */ _col(`c1') "`name1'" _col(49) /*
          */ `": Number of obs = "' as res %14.0g `n1'
        noi: di as txt _col(`c2') "`name2'" _col(49) /*
          */ `": Number of obs = "' as res %14.0g `n2'
        noi: _tostprtest0 `varlist' `n1' `m1' `se1' `alpha'
        noi: _tostprtest2 `exp' `n2' `m2' `se2' `alpha'
        }
      noi: di in smcl as txt "{hline 13}{c +}{hline 66}"
      local se = sqrt(`se1'^2 + `se2'^2)
      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"
        }
      noi: _tostprtest2 "`theta'`caret'" `N' `m1'-`m2' `se' `alpha'
      noi: di as txt _col(14) "{c |}  under Ho:" as res %9.0g _col(28) `se_p' %5.0g _col(40) `z_pos' %06.4f _col(48) `p_pos'
      noi: di in smcl as txt "{hline 13}{c BT}{hline 66}"
      if "`by'" == "" {
        noi: di as txt _col(7) "`theta_inflate'`theta'`caret' = prop(" abbrev(`"`varlist'"', 12) "=1) - prop(`exp'=1) = "  as res %9.0g `m1' - `m2' as txt _col(70) "z = " as res %-8.4f _col(73) `z_pos'
        }
      if "`by'" != "" {
        noi: di as txt _col(7) "`theta_inflate'`theta'`caret' = prop(" abbrev(`"`varlist'"', 12) "|`by'=`name1'" as txt ") - prop(" abbrev(`"`varlist'"', 12) "|`by'=`name2'" as txt ")" as txt _col(70) "z = " as res %-8.4f _col(73) `z_pos'
        noi: di as txt _col(17) "= " as res %9.0g `m1' - `m2'
        }
      noi: di as txt _col(4) "`theta_inflate'Ho: `theta' = 0" _newline
      noi: `correction_message'
      noi: di as txt _col(37) "Ha: `theta' `notequal' 0"
      noi: di as txt _col(31) "Pr(|Z| > |z|) " `p_pos_display' _newline
      if (1-normal(abs(((`m1' - `m2') - `continuity') / `se_p'))) > (`alpha'/2) {
        local PositivistConclusion = "Fail to reject"
        }
      }
    if lower("`eqvtype'") == "delta" {
      local z1 = (`upper' - (sign(`m1' - `m2') * ((`m1' - `m2') - `continuity'))) / `se_p'
      local z2 = ((sign(`m1' - `m2') * ((`m1' - `m2') - `continuity')) + `lower') / `se_p'
      }
    if lower("`eqvtype'") == "epsilon" {
      local z1 = `upper' - `z_pos'
      local z2 = `z_pos' + `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"
      }
    if "`by'" != "" {
      local name1 = "`: label (`by')`min''"
      local name2 = "`: label (`by')`max''"
      local c1 = 48 - length("`by'") - length(`"`name1'"')
      local c2 = 48 - length("`by'") - length(`"`name2'"')
      noi: di as txt _newline "Two-sample z test for proportion equivalence" _newline /*
        */ _col(`c1') "`by'=`name1'" _col(48) /*
        */ `": Number of obs = "' as res %14.0g `n1'
      noi di as txt _col(`c2') "`by'=`name2'" _col(48) /*
        */ `": Number of obs = "' as res %14.0g `n2'
      noi: _tostprtest1 "`by'=`name1'" `n1' `m1' `se1' `alpha'
      noi: _tostprtest2 "`by'=`name2'" `n2' `m2' `se2' `alpha'
      }
    if "`by'" == "" {
      local name1 = trim(substr(trim(`"`varlist'"'),1,12))
      local name2 = trim(substr(trim(`"`exp'"'),1,12))
      local c1 = 49 - length(`"`name1'"')
      local c2 = 49 - length(`"`name2'"')
      noi: di as txt _newline "Two-sample z test for proportion equivalence" _newline /*
        */ _col(`c1') "`name1'" _col(49) /*
        */ `": Number of obs = "' as res %14.0g `n1'
      noi: di as txt _col(`c2') "`name2'" _col(49) /*
        */ `": Number of obs = "' as res %14.0g `n2'
      noi: _tostprtest1 `varlist' `n1' `m1' `se1' `alpha'
      noi: _tostprtest2 `exp' `n2' `m2' `se2' `alpha'
      }
    noi: di in smcl as txt "{hline 13}{c +}{hline 66}"
    if lower("`eqvtype'") == "delta" {
      if (`upper' == `lower') {
        noi: _tostprtest2 "`delta'-`theta'`caret'" `N' `upper'-(`m1'-`m2') `se_p' `alpha'
        noi: _tostprtest2 "`theta'`caret'+`delta'" `N' (`m1'-`m2')+`lower' `se_p' `alpha'
        }
      if (`upper' != `lower') {
        noi: _tostprtest2 "`delta_upper'-`theta'`caret'" `N' `upper'-(`m1'-`m2') `se_p' `alpha'
        noi: _tostprtest2 "`theta'`caret'-`delta_lower'" `N' (`m1'-`m2')+`lower' `se_p' `alpha'
        }
      }
    if lower("`eqvtype'") == "epsilon" {
      noi: _tostprtest2 "`theta'`caret'" `N' `m1'-`m2' `se_p' `alpha'
      }
    if (`p1' > `alpha' | `p2' > `alpha') {
      local NegativistConclusion = "Fail to reject"
      }
    noi: di in smcl as txt "{hline 13}{c BT}{hline 66}"
    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 neg_p_vals_out noi: di as txt _col(4) "Pr(Z > z1) " `p1_display' _col(33) as txt "Pr(Z > z2) " as res %6.4f `p2_display'
    if lower("`eqvtype'") == "delta" {
      if "`by'" == "" {
        noi: di as txt _col(7) "`theta_inflate'`theta'`caret' = prop(" abbrev(`"`varlist'"', 12) "=1) - prop(`exp'=1) = "  as res %9.0g `m1' - `m2'
        }
      if "`by'" != "" {
        noi: di as txt _col(7) "`theta_inflate'`theta'`caret' = prop(" abbrev(`"`varlist'"', 12) "|`by'=`name1'" as txt ") - prop(" abbrev(`"`varlist'"', 12) "|`by'=`name2'" as txt ")"
        noi: di as txt _col(17) "= " as res %9.0g `m1' - `m2'
        }
      if (`upper' == `lower') {
        noi: di as txt _col(8) "`delta_inflate'`delta' = " as res %-6.3f `lower' as txt _col(24) "`delta' expressed in same units as prop(" abbrev(`"`varlist'"', 12) as txt "=1)"
        }
      if (`upper' != `lower') {
        noi: di as txt _col(11) "`delta_lower' = " as res %-7.3g -1*`lower' as txt _col(23) "`delta_lower' expressed in same units as prop(" abbrev(`"`varlist'"', 12) as txt "=1)"
        noi: di as txt _col(11) "`delta_upper' =  " as res %-6.3g `upper' as txt _col(23) "`delta_upper' expressed in same units as prop(" abbrev(`"`varlist'"', 12) as txt "=1)"
        }
      local criticalvalue = `se_p'*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 tostpr##mineqvlevel: help tostpr}."
        }
      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 tostpr##mineqvlevel: help tostpr}."
        }
      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 tostpr##mineqvlevel: help tostpr}."
        }
      if (`upper' == `lower') {
        noi: di _newline as txt "Ho: |`theta'| `greaterthan' `delta':" _newline
        noi: `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"
        `neg_p_vals_out'
        }
      if (`upper' != `lower') {
        noi: di _newline as txt "Ho: `theta' `lessthan' `delta_lower', or `theta' `greaterthan' `delta_upper':" _newline
        noi: `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"
        `neg_p_vals_out'
        }
      }
    if lower("`eqvtype'") == "epsilon" {
      if "`by'" == "" {
        noi: di as txt _col(7) "`theta_inflate'`theta'`caret' = prop(" abbrev(`"`varlist'"', 12) "=1) - prop(`exp'=1) = "  as res %9.0g `m1' - `m2'
        }
      if "`by'" != "" {
        noi: di as txt _col(7) "`theta_inflate'`theta'`caret' = prop(" abbrev(`"`varlist'"', 12) "|`by' = `name1') - prop(" abbrev(`"`varlist'"', 12) "|`by' = `name2')"
        noi: di as txt _col(17) "= " as res %9.0g `m1' - `m2'
        }
      if (`upper' == `lower') {
        noi: di as txt _col(6) "`epsilon_inflate'`epsilon' = " as res %-5.3f `lower' as txt _col(24) "`epsilon' expressed in units of the z distribution"
        }
      if (`upper' != `lower') {
        noi: di as txt _col(11) "`epsilon_lower' = " as res %-6.3f -1*`lower' as txt _col(23) "`epsilon_lower' expressed in units of the z distribution"
        noi: di as txt _col(11) "`epsilon_upper' =  " as res %-5.3f `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 tostpr##mineqvlevel: help tostpr}."
        }
      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 tostpr##mineqvlevel: help tostpr}."
        }
      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 tostpr##mineqvlevel: help tostpr}."
        }
      if (`upper' == `lower') {
        noi: di _newline as txt "Ho: |Z| `greaterthan' `epsilon':" _newline
        noi: `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"
        `neg_p_vals_out'
        }
      if (`upper' != `lower') {
        noi: di _newline as txt "Ho: Z `lessthan' `epsilon_lower', or Z `greaterthan' `epsilon_upper':" _newline
        noi: `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"
        `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" `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 = " 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'"
    }
  capture confirm number `exp'
  if (`upper' == `lower') {
    if "`eqvtype'" == "delta" {
      return scalar Delta   = `eqvlevel'
      }
    if "`eqvtype'" == "epsilon" {
      return scalar epsilon   = `eqvlevel'
      }
    }
  if (`upper' != `lower') {
    if "`eqvtype'" == "delta" {
      return scalar Dl   = `lower'
      return scalar Du   = `upper'
      }
    if "`eqvtype'" == "epsilon" {
      return scalar el   = `lower'
      return scalar eu   = `upper'
      }
    }
  return scalar N_1     = `n1'
  return scalar P_1     = `m1'
  if _rc != 0 {
    return scalar N_2   = `n2'
    return scalar P_2   = `m2'
    }
  return scalar z1      = `z1'
  return scalar z2      = `z2'
  end


program define _tostprtest0
 local name = abbrev(`"`1'"', 12)
 local n "`2'"
 local mean "`3'"
 local se "`4'"
 local alpha "`5'"
 local show = "`6'"
 if `n' == 1 | `n' >= . {
   local se = .
   }
 local beg = 13 - length(`"`name'"')
 if "`show'" != "" {
   local z z
 	 local zp P>|Z|
   }
 local level : di %3.2f (100*(1-`alpha'))
 if (mod(`level'*10, 1) == 0.0) {
   local level: di %3.1f `level'
   }
 if (mod(`level', 1) == 0.0) {
   local level: di %2.0f `level'
   }
 local cil `=string(`level')'
 local cil `=length("`cil'")'
 noi: di in smcl as txt "{hline 13}{c TT}{hline 66}"
 noi: di in smcl as txt _col(5) "Variable {c |}" /*
 */ _col(22) "Mean" _col(29) /*
 */ "Std. Err." _col(42) "z" _col(46) /*
 */ "Pr(|Z|>|z|)" _col(`=63-`cil'') `"[`=strsubdp("`level'")'% Conf. Interval]"'
 noi: di in smcl as txt "{hline 13}{c +}{hline 66}"
 local vval = (100-(100-`level')/2)/100
 noi: di in smcl as txt _col(`beg') `"`name'"' /*
 */ as txt _col(14) "{c |}" as res /*
 */ _col(17) %9.0g  `mean'   /*
 */ _col(28) %9.0g  `se'     /*
 */ _col(60) %9.0g  `mean'-invnorm(`vval')*`se'   /*
 */ _col(72) %9.0g  `mean'+invnorm(`vval')*`se'
 end

program define _tostprtest1
 local name = abbrev(`"`1'"', 12)
 local n "`2'"
 local mean "`3'"
 local se "`4'"
 local alpha "`5'"
 local show = "`6'"
 if `n' == 1 | `n' >= . {
   local se = .
   }
 local beg = 13 - length(`"`name'"')
 if "`show'" != "" {
   local z z
 	 local zp P>|Z|
   }
 local level : di %3.2f (100*(1-`alpha'))
 if (mod(`level'*10, 1) == 0.0) {
   local level: di %3.1f `level'
   }
 if (mod(`level', 1) == 0.0) {
   local level: di %2.0f `level'
   }
 local cil `=string(`level')'
 local cil `=length("`cil'")'
 noi: di in smcl as txt "{hline 13}{c TT}{hline 66}"
 noi: di in smcl as txt _col(5) "Variable {c |}" /*
 */ _col(22) "Mean" _col(29) /*
 */ "Std. Err." _col(44) "`z'" _col(49) /*
 */ "`zp'" _col(`=63-`cil'') `"[`=strsubdp("`level'")'% Conf. Interval]"'
 noi: di in smcl as txt "{hline 13}{c +}{hline 66}"
 local vval = (100-(100-`level')/2)/100
 noi: di in smcl as txt _col(`beg') `"`name'"' /*
 */ as txt _col(14) "{c |}" as res /*
 */ _col(17) %9.0g  `mean'   /*
 */ _col(28) %9.0g  `se'     /*
 */ _col(60) %9.0g  `mean'-invnorm(`vval')*`se'   /*
 */ _col(72) %9.0g  `mean'+invnorm(`vval')*`se'
 end

program define _tostprtest2
 local name = abbrev(`"`1'"', 12)
 local n "`2'"
 local mean "`3'"
 local se "`4'"
 if `n' == 1 | `n' == . {
   local se = .
   }
 local alpha "`5'"
 local level : di %3.2f (100*(1-`alpha'))
 if (mod(`level'*10, 1) == 0.0) {
   local level: di %3.1f `level'
   }
 if (mod(`level', 1) == 0.0) {
   local level: di %2.0f `level'
   }
 local vval = (100 - (100-`level')/2)/100
 noi: di in smcl as txt %12s `"`name'"' " {c |}" as res /*
 		*/ _col(17) %9.0g  `mean'   /*
		*/ _col(28) %9.0g  `se'     /*
		*/ _col(60) %9.0g  `mean'-invnorm(`vval')/*
		*/ *`se'   /*
		*/ _col(72) %9.0g  `mean'+invnorm(`vval')*`se'
 end