*! rdcii.ado Version 1.2 2010-09-30 JRC
program define rdcii, rclass
    version 9.2
    syntax anything(id="argument numlist") [, Level(real `c(level)') Zsot Cc noBRUTEforce ///
      TOLerance(real 1e-6) LTOLerance(real 0) Verbose INITial(numlist ascending min=1 max=2)]

    tempname z z2 n0 n1 p0 p1 rd lb ub bracket chi2rd

    tokenize `anything'
    if ("`1'" == "miemar") `1' `2' `3' `4' `5' `6' `7' `8'
    else {
        local variable_tally : word count `anything'
        if (`variable_tally' > 4) exit = 103
        if (`variable_tally' < 4) exit = 102
*
        forvalues i = 1/4 {
            capture confirm integer number ``i''
            if _rc {
                display in smcl as error "`id' must all be numeric."
                exit = 499
            }
        }
        forvalues i = 1/4 {
            capture assert ``i'' >= 0
            if _rc {
                display in smcl as error "`id' must all be nonnegative."
                exit = 499
            }
        }
        if !inrange(`level', 0.1, 99.9) {
            display in smcl as error "Level must lie between 0.1 and 99.9."
            exit = 499
        }
*
        scalar define `z' = invnormal(1 - (100 - `level') / 200 )
        scalar define `z2' = invchi2(1, `level' / 100)
        scalar define `n0' = `2' + `4'
        scalar define `n1' = `1' + `3'
        scalar define `p0' = `2' / `n0'
        scalar define `p1' = `1' / `n1'
        scalar define `rd' = `p1' - `p0'
        scalar define `lb' = .
        scalar define `ub' = .
*
        capture assert ( (`n0' > 0) & (`n1' > 0) )
        if _rc {
            display in smcl as error "Each exposure group must have at least one observation."
            exit = 499
        }
*
        return scalar p0 = `p0'
        return scalar p1 = `p1'
        return scalar rd = `rd'
*
* Agresti-Caffo confidence interval
*
        agresticaffo `1' `2' `n1' `n0', lb(`lb') ub(`ub') z(`z') z2(`z2') `zsot'
        return scalar lb_ac = `lb'
        return scalar ub_ac = `ub'
*
* Newcombe's Method 10 confidence interval
*
        newcombe `1' `2' `n1' `n0' `p1' `p0' `rd', lb(`lb') ub(`ub') z(`z') z2(`z2')
        return scalar lb_ne = `lb'
        return scalar ub_ne = `ub'
*
* Wallenstein's confidence interval
*
        wallenstein `1' `2' `n1' `n0' `p1' `p0' `rd', lb(`lb') ub(`ub') z(`z') z2(`z2') `cc'
        return scalar lb_wa = `lb'
        return scalar ub_wa = `ub'
*
* Miettinen-Nurminen confidence interval
*
        scalar define `chi2rd' = .
        if (`tolerance' != 1.000e-6) local tolerance tol `tolerance'
        else local tolerance
*
        if (`ltolerance' != 0) local ltolerance ltol `ltolerance'
        else local ltolerance
*
        if ("`verbose'" != "") local verbose noisily
*
        if ("`initial'" != "") {
            gettoken lower upper : initial
            scalar define `bracket' = `lower'
        }
        else scalar define `bracket' = `rd'
*
        local one = 0.9999999999
        local minus_one = -.9999999999
*
        if (`rd' == -1) return scalar lb_mn = `rd'
        else {
            capture `verbose' ridder rdcii "miemar" `1' `2' `n1' `n0' `rd' `chi2rd' X returns exp `chi2rd' = `z2' ///
              from `minus_one' to `bracket' `tolerance' `ltolerance'
            if !_rc return scalar lb_mn = $S_1
            else {
                return scalar lb_mn = .
                if ("`bruteforce'" == "") {
                    forvalues i = -1(0.01)`=`bracket'' {
                        capture `verbose' ///
                          ridder rdcii "miemar" `1' `2' `n1' `n0' `rd' `chi2rd' X returns exp `chi2rd' = `z2' ///
                          from `minus_one' to `i' `tolerance' `ltolerance'
                        if !_rc {
                            return scalar lb_mn = $S_1
                            continue, break
                        }
                    }
                }
            }
        }
*
        if (`rd' == 1) return scalar ub_mn = `rd'
        else {
*
            if ("`upper'" != "") scalar define `bracket' = `upper'
*
            capture `verbose' ridder rdcii "miemar" `1' `2' `n1' `n0' `rd' `chi2rd' X returns exp `chi2rd' = `z2' ///
              from `bracket' to `one' `tolerance' `ltolerance'
            if !_rc return scalar ub_mn = $S_1
            else {
                return scalar ub_mn = .
                if ("`bruteforce'" == "") {
                    forvalues i = `=`bracket''(0.01)1 {
                        capture `verbose' ///
                          ridder rdcii "miemar" `1' `2' `n1' `n0' `rd' `chi2rd' X returns exp `chi2rd' = `z2' ///
                          from `i' to `one' `tolerance' `ltolerance'
                        if !_rc {
                            return scalar ub_mn = $S_1
                            continue, break
                        }
                    }
                }
            }
        }
        macro drop S_1
*
*
*
        display in smcl as text _newline(1) "Confidence intervals for risk difference" _newline(1)
        display in smcl as text "       Risk for unexposed (p0): " as result %05.3f return(p0)
        display in smcl as text "         Risk for exposed (p1): " as result %05.3f return(p1)
        display in smcl as text "     Risk difference (p1 - p0): " as result %05.3f return(rd)
        display in smcl as text _newline(1) "{hline 41}"
        display in smcl as text _column(10) "Method" _column(20) %22s "[`=string(`level', "%5.0g")'% Conf. Interval]"
        display in smcl as text "{hline 41}"
        display in smcl as text %18s "Agresti-Caffo" ///
          _column(`=cond(return(lb_ac) < 0, 22, 23)') as result %05.3f return(lb_ac) ///
          _column(`=cond(return(ub_ac) < 0, 36, 37)') as result %05.3f return(ub_ac)
        display in smcl as text %18s "Newcombe Method 10" ///
          _column(`=cond(return(lb_ne) < 0, 22, 23)') as result %05.3f return(lb_ne) ///
          _column(`=cond(return(ub_ne) < 0, 36, 37)') as result %05.3f return(ub_ne)
        display in smcl as text %18s "Wallenstein" ///
          _column(`=cond(return(lb_wa) < 0, 22, 23)') as result %05.3f return(lb_wa) ///
          _column(`=cond(return(ub_wa) < 0, 36, 37)') as result %05.3f return(ub_wa)
        display in smcl as text %18s "Miettinen-Nurminen" ///
          _column(`=cond(return(lb_mn) < 0, 22, 23)') as result %05.3f return(lb_mn) ///
          _column(`=cond(return(ub_mn) < 0, 36, 37)') as result %05.3f return(ub_mn)
        display in smcl _newline(1)
    }
end

program define agresticaffo
    version 9.2
    syntax anything, lb(name) ub(name) z(name) z2(name) [zsot]

    tempname tot N20 N21 c10 c11 rd zse tm tp

    tokenize `anything'

    if ("`zsot'" != "") scalar define `tot' = `z2' / 2 // T Over Two
    else scalar define `tot' = 2

    scalar define `N20' = `4' + `tot'
    scalar define `N21' = `3' + `tot'
    scalar define `c10' = `2' + `tot' / 2
    scalar define `c11' = `1' + `tot' / 2

    scalar define `rd' = `c11' / `N21' - `c10' / `N20'
    scalar define `zse' = `z' * sqrt( `c11' * (`N21' - `c11') / `N21' / `N21' / `N21' + ///
      `c10' * (`N20' - `c10') / `N20' / `N20' / `N20' )

    scalar define `tm' = `rd' - `zse'
    scalar define `tp' = `rd' + `zse'

    scalar define `lb' = max( -1, min(`tm', `tp') )
    scalar define `ub' = min(  1, max(`tm', `tp') )
end

program define newcombe
    version 9.2
    syntax anything, lb(name) ub(name) z(name) z2(name)

    tempname a b c lb1 lb2 ub1 ub2 cl1 cl2
    tokenize `anything'

* Wilson intervals of the individual proportions (adapted from -ciwi- by N. J. Cox)
    forvalues i = 1/2 {
        local n`i' = ``=`i' + 2''
        local p = `i' + 4
        scalar define `a' = 2 * ``i'' + `z2'
        scalar define `b' = `z' * sqrt( `z2' + 4 * `n`i'' * ``p'' * (1 - ``p'') )
        scalar define `c' = 2 * (`n`i'' + `z2')
        scalar define `lb`i'' = (`a' - `b') / `c'
        scalar define `ub`i'' = (`a' + `b') / `c'
    }
*
    scalar define `cl1' = `7' - `z' * sqrt( `lb1' * (1 - `lb1') / `n1' + `ub2' * (1 - `ub2') / `n2' )
    scalar define `cl2' = `7' + `z' * sqrt( `ub1' * (1 - `ub1') / `n1' + `lb2' * (1 - `lb2') / `n2' )

    if missing(`cl1') {
        scalar define `cl1' = sign(`7')
    }

    if missing(`cl2') {
        scalar define `cl2' = sign(`7')
    }

    scalar define `lb' = min(`cl1', `cl2')
    scalar define `ub' = max(`cl1', `cl2')
end

program define wallenstein
    version 9.2
    syntax anything, lb(name) ub(name) z(name) z2(name) [recursion(integer 0) cc]

    tempname epsilon delta p_bar d1 d2 p1dl p1du p2du minus_seven swap
    tokenize `anything'

    local N = `3' + `4'
    scalar define `p_bar' = (`1' + `2') / `N'
    scalar define `d1' = .
    scalar define `d2' = .
    if ("`cc'" == "") {
        wallenquad , d1(`d1') d2(`d2') three(`3') four(`4') delta(`7') z2(`z2') p_bar(`p_bar') en(`N')
        scalar define `lb' = min(`d1', `d2')
        scalar define `ub' = max(`d1', `d2')
        scalar define `epsilon' = 0
    }
    else {
        scalar define `epsilon' = 1 / 2 / `3' + 1 / 2 / `4'
        scalar define `delta' = sign(`7') * min(1, abs(`7' + `epsilon'))
        wallenquad , d1(`d1') d2(`d2') three(`3') four(`4') delta(`delta') z2(`z2') p_bar(`p_bar') en(`N')
        scalar define `ub' = max(`d1', `d2')
        scalar define `delta' = sign(`7') * min(1, abs(`7' - `epsilon'))
        wallenquad , d1(`d1') d2(`d2') three(`3') four(`4') delta(`delta') z2(`z2') p_bar(`p_bar') en(`N')
        scalar define `lb' = min(`d1', `d2')
    }
*
    if (`5' >= `6') {
        scalar define `p1du' = `p_bar' + `ub' * `4' / `N'
        scalar define `p2du' = `p_bar' - `ub' * `3' / `N'
        scalar define `p1dl' = `p_bar' + `lb' * `4' / `N'
        if (`p2du' < 0) {
            scalar define `ub' = `7' + `epsilon' + `z2' / 2 / `3' + (`z' / sqrt(`3')) * ///
                sqrt((`7' + `epsilon') * (1 - (`7' + `epsilon')) + `z2' / 4 / `3')
            scalar define `ub' = `ub' / (1 + `z2' / `3')
        }
        if (`p1du' > 1) {
            scalar define `ub' = `7' + `epsilon' + `z2' / 2 / `4' + (`z' / sqrt(`4')) * ///
                sqrt((`7' + `epsilon') * (1 - (`7' + `epsilon')) + `z2' / 4 / `4')
            scalar define `ub' = `ub' / (1 + `z2' / `4')
        }
        if (`p1dl' < 0) {
            scalar define `lb' = `7' - `epsilon' + `z2' / 2 / `3' + (`z' / sqrt(`3')) * ///
                sqrt((`7' - `epsilon') * (1 - (`7' - `epsilon')) + `z2' / 4 / `3')
            scalar define `lb' = `lb' / (1 + `z2' / `3')
            scalar define `lb' = -`lb'
        }
        if (`recursion' == 1) {
            scalar define `swap' = `ub'
            scalar define `ub' = -`lb'
            scalar define `lb' = -`swap'
        }
    }
    else { // "If p_hat2 > p_hat1, the above logic is somewhat reversed, and it may be easier simply to reverse indices."
        scalar define `minus_seven' = -`7'
        if (`recursion' == 0) wallenstein `2' `1' `4' `3' `6' `5' `minus_seven', ///
          lb(`lb') ub(`ub') z(`z') z2(`z2') recursion(1) `cc'
    }
end

program define wallenquad
    version 9.2
    syntax , d1(name) d2(name) three(name) four(name) delta(name) z2(name) p_bar(name) en(integer)

    tempname a b c d

    scalar define `a' = 1 + `z2' / `en' * ( 1 + (`three' - `four')^2 / `three' / `four' )
    scalar define `b' = -2 * `delta' + `z2' * (1 - 2 * `p_bar') * (`three' - `four') / `three' / `four'
    scalar define `c' = `delta' * `delta' - `z2' * `en' * `p_bar' * (1 - `p_bar') / `three' / `four'
    scalar define `d' = -( `b' + ( 2 * (`b' >= 0) - 1 )* sqrt(`b' * `b' - 4 * `a' * `c') ) / 2
    scalar define `d1' = `d' / `a'
    scalar define `d2' = `c' / `d'
end

program define miemar
    version 9.2
    syntax anything

    tempname L3 L2 L1 L0 L23 q p a R0 R1 V

    tokenize `anything'

    local c = `1' + `2'
    local S = `3' + `4'

    scalar define `L3' = `S'
    scalar define `L2' = (`4' + 2 * `3') * (-`7') - `S' - `c'
    scalar define `L1' = (`3' * (-`7') - `S' - 2 * `1') * (-`7') + `c'
    scalar define `L0' = `1' * (-`7') * (1 - (-`7'))

    scalar define `L23' = `L2' / 3 / `L3'

    scalar define `q' = `L23'^3 - `L1' * `L2' / 6 / `L3'^2 + `L0' / 2 / `L3'
    scalar define `p' = sign(`q') * sqrt(`L23'^2 - `L1' / 3 / `L3')
    scalar define `a' = ( _pi + acos(`q' / `p'^3) ) / 3
    scalar define `R0' = 2 * `p' * cos(`a') - `L23'
    scalar define `R1' = `R0' + (-`7')
    scalar define `V' = ( `R1' * (1 - `R1') / `4' + `R0' * (1 - `R0') / `3' ) * `S' / (`S' - 1)
    scalar define `6' = (`5' - `7')^2 / `V'
end