*! version 3.1.5  12oct2024 by alexis.dinno@pdx.edu
*! perform two one-sided tests for equivalence of RR and unity in paired binary data

* Syntax:  tostrrpi #a #b #c #n #delta0, [ deltaupper(#) alpha(#) exactchisq
* treatment1(string) treatment2(string) outcome(string) nooutcome(string) relevance]

program define tostrrpi, rclass

  if int(_caller())<8 {
    di in r "tostrrpi- 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           = "delta0"
    local delta_lower     = "dl"
    local delta_upper     = "du"
    local delta_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 delta           = uchar(0948)+uchar(8320)
    local delta_lower     = uchar(0948)+uchar(0108)
    local delta_upper     = uchar(0948)+uchar(0117)
    local delta_inflate   = "    "
    local chisq           = uchar(0967)+uchar(8306)
    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 n 0 : 0, parse(" ,")
  gettoken delta0 0 : 0, parse(" ,")

  confirm integer number `a'
  confirm integer number `b'
  confirm integer number `c'
  confirm integer number `n'
  confirm number `delta0'
  local d = `n' - `a' - `b' - `c'

  syntax [, DELTAUpper(real 0) Alpha(real 0.05) EXactchisq  /*
  */      treatment1(string) treatment2(string) OUTcome(string) /*
  */      NOOutcome(string) RELevance]

* Validate delta0
  if (`delta0' < 0 | `delta0' > 1 ) {
    noisily: di as err "delta0 must be a positive real value that is less than or equal to 1"
    exit 198
    }

* Validate deltaupper
  if (`deltaupper' < 1 & `deltaupper' !=0) {
    noisily: di as err "deltaupper must be a positive real value that is greater than or equal to 1"
    exit 198
    }
   else {
    local upper = 1/`deltaupper'
	}
  if `deltaupper' == 0 {
    local upper = `delta0'
    }
  local lower = `delta0'

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

* Set up treatment and outcome status names, as user requests
  if "`treatment1'"=="" {
    local treatname1 = "Treatment 1"
   	}
   else {
    local treatname1 = "`treatment1'"
    }
  local treatname1 = substr("`treatname1'",1,9)
  local treatname1long = substr("`treatment1'",1,16)

  if "`treatment2'"=="" {
    local treatname2 = "Treatment 2"
    }
   else {
    local treatname2 = "`treatment2'"
    }
  local treatname2 = substr("`treatname2'",1,9)
  local treatname2long = substr("`treatment2'",1,16)

  if "`outcome'" =="" {
    local positive = "Positive"
    local rowpositive = "Positive"
    local colpositive = "Positive"
    }
   else {
    local positive = substr("`outcome'",1,9)
    local rowpositive = substr("`outcome'",1,16)
    local colpositive = substr("`outcome'",1,9)
    }
  if "`nooutcome'" =="" {
    local negative = "Negative"
    local rownegative = "Negative"
    local colnegative = "Negative"
    }
   else {
    local rownegative = substr("`nooutcome'",1,16)
    local colnegative = substr("`nooutcome'",1,9)
    }


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

**********
* conduct McNemar's test of proportion difference (because if OR != 1, RR != 1)
  local PositivistConclusion = "Reject"
  local NegativistConclusion = "Reject"
  local chisq_pos = ((`b' - `c')^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"
    }
  local chisq_p_val_notice di as txt "" _continue
  * Use exact chi-sq p value if requested
  if ("`exactchisq'"!="") {
    quietly {
      local level = 100*(1 - `alpha')
      mcci `a' `b' `c' `d', level(`level')
      local p_pos = r(p_exact)
      local chisq_p_val_notice di as txt _newline "Using exact p-value for McNemar`rquote's `chisq'"
      }
    }
  if (`p_pos' > `alpha'/2) {
    local PositivistConclusion = "Fail to reject"
    }
  if ("`relevance'" != "") {
    di as txt _newline "Relevance test for relative risk and unity in paired designs" _newline
    noi: di as txt "McNemar`rquote's test for difference between relative risk and unity in paired designs"
    di as txt "{hline 17}{c TT}{hline 24}{c TT}{hline 11}"

    _crc4fld `a' `b' `c' `d' "`treatname1long'" "`treatname2long'" "`colpositive'" "`colnegative'" "`rowpositive'" "`rownegative'"
    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') " as res %6.4f `p_pos_display'
    `chisq_p_val_notice'
    noi: di as txt "Reject Ho if p `lessthan' `alpha'/2 (two-sided test)"
    noi: di as txt _newline "Ho: Pr(`positive'|`treatname1') = Pr(`positive'|`treatname2')"
    noi: di as txt "Ha: Pr(`positive'|`treatname1') `notequal' Pr(`positive'|`treatname2')"
   di ""
   }
  if `a' == 0 & `b' == 0 & `c' == 0 {
    local RR = .
   	local sdRR = .
    local z1 = .
	   local z2 = .
	   noi: di as res _newline "Note: neither treatment has any observed responses; relative risk undefined"
   	}
   else if (`a'> 0 & `b' == 0 & `c' == 0) {
    local RR = 1
   	local riskprefix = ""
	   local sdRR = sqrt((`a'+`b')*(`b'+`c')/(`a'+`c')^3)
    local z1 = sqrt(`a'*(1-`lower')/`lower')
    local z2 = sqrt(`a'*(1-`upper')/`upper')
    noi: di as res _newline "Note: treatments have complete concordance"
    }
   else {
    local RR = (`a'+`b')/(`a'+`c')
   	local riskprefix = "0"
	   local sdRR = sqrt((`a'+`b')*(`b'+`c')/(`a'+`c')^3)
    _TangTangChanZ `a' `b' `c' `n' `lower'
    local se = r(se)
	   local z1 = r(TTC)
    _TangTangChanZ `a' `c' `b' `n' `upper'
	   local z2 = r(TTC)
    }

  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"
    }
  di as txt _newline "Tang-Tang-Chan test for equivalence of relative risk and unity in paired designs"
  di as txt "{hline 17}{c TT}{hline 24}{c TT}{hline 12}"
  _crc4fld `a' `b' `c' `d' "`treatname1long'" "`treatname2long'" "`colpositive'" "`colnegative'" "`rowpositive'" "`rownegative'"
  di as txt "{hline 17}{c BT}{hline 24}{c BT}{hline 12}"

  noisily: di as text "Relative risk of `positive' in `treatname2' vs. `treatname1':"
  noisily: di as text _col(9) "RR = " as res `a'+`b' as text " / " as res `a'+`c' as text " = "  as res "`riskprefix'" %-8.4g `RR'
  noisily: di as text _col(4) "s.e. RR = " as res "0" %-8.4g `sdRR'
  noisily: di as text _col(5) "`delta_inflate'`delta_lower' = " as res  %-8.4f `lower'
  noisily: di as text _col(5) "`delta_inflate'`delta_upper' = " as res  %-8.4f 1/`upper'

  noisily: di _newline as txt "Ho: true RR `lessthan' `delta_lower', or true RR `greaterthan' `delta_upper':" _newline
  noisily: di as txt _col(9) "z1 = " as res %-8.4g `z1' as txt _col(38) "z2 = " as res %-8.4g `z2' _newline
  noisily: di as txt _col(4) "Ho1: true RR `lessthan' `delta_lower'" _col(33) "Ho2: true RR `greaterthan' `delta_upper'"
  noisily: di as txt _col(4) "Ha1: true RR > `delta_lower'"  _col(33) "Ha2: true RR < `delta_upper'"
  noisily: di as txt _col(4) "Pr(Z > z1) " as res %6.4f `p1_display' _col(33) as txt "Pr(Z > z2) " `p2_display'


* 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 delta0 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') {
     noi: di _newline _newline as txt "Relevance test conclusion for `alpha_out' = " as res "0" `alpha' as txt ", `delta_lower' =" as res %5.2f `lower' as txt ", and `delta_upper' =" as res %5.2f 1/`lower' as txt ":"
     }
   if (`upper' != `lower') {
     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 1/`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') {
    return scalar delta0   = `lower'
    }
  if (`upper' != `lower') {
    return scalar deltaupper   = 1/`upper'
    return scalar deltalower   = `lower'
    }
  return scalar   p2   = `p2'
  return scalar   p1   = `p1'
  return scalar   z2   = `z2'
  return scalar   z1   = `z1'
  return scalar sdRR   = `sdRR'
  return scalar   RR   = `RR'

  end

  *
program define _q1, rclass sort
/*
    _q1 returns the primary root of the "q1" term of f(x) at the top of page
	1221 in Tang, Tang, and Chan (2003)

    Syntax:
            _q1 #a #b #c #n #delta0
*/

  gettoken a 0 : 0, parse(" ,")
  gettoken b 0 : 0, parse(" ,")
  gettoken c 0 : 0, parse(" ,")
  gettoken n 0 : 0, parse(" ,")
  gettoken delta0 0 : 0, parse(" ,")

  confirm integer number `a'
  confirm integer number `b'
  confirm integer number `c'
  confirm integer number `n'
  confirm number `delta0'

  local A = `n'*(1+`delta0')
  local B = (`a'+`c')*(`delta0'^2)-(`a'+`b'+2*`c')
  local C = `c'*(1-`delta0')*(`a'+`b'+`c')/`n'
  local q1 = (sqrt(((`B')^2) - (4*`A'*`C')) - `B')/(2*`A')

  ret scalar q1 = `q1'

  end


program define _TangTangChanZ, rclass sort


  gettoken a 0 : 0, parse(" ,")
  gettoken b 0 : 0, parse(" ,")
  gettoken c 0 : 0, parse(" ,")
  gettoken n 0 : 0, parse(" ,")
  gettoken delta0 0 : 0, parse(" ,")

  confirm integer number `a'
  confirm integer number `b'
  confirm integer number `c'
  confirm integer number `n'
  confirm number `delta0'

  local numerator = `a' + `b' - (`a' + `c')*`delta0'
  _q1 `a' `b' `c' `n' `delta0'
  local q1 = r(q1)
  local denominator = sqrt(`n'*((1+`delta0')*`q1' + ((`a'+`b'+`c')*(`delta0'-1)/`n')))
  local TTC = `numerator'/`denominator'

  return scalar TTC = `TTC'
  return scalar se = `denominator'

  end