*! version 2.0.1 21feb2021 daniel klein
program rioc , rclass byable(recall)
    version 11.2
    
    syntax varlist(numeric min=2 max=2) ///
    [ if ] [ in ] [ fweight ]           ///
    [ ,                                 ///
        Tab                             ///
        Detail                          ///
        STATs(passthru)                 ///
        CII(name)                       /// not documented        
        KAPPA                           ///
        SMALLsample                     ///
        SEZERO                          ///
        ASR /*synonym*/ CHi2  /// retained; not documented        
        Level(cilevel)                  ///
        CFORMAT(passthru)               ///
        PFORMAT(passthru)               ///
        SFORMAT(passthru)               ///
        PERCENT                         /// not documented
        MATCELL(name)                   /// not documented        
    ]
    
    marksample touse
    
    /*
        the notation below follows Copas and Loeber (1990)
        
        we create 0/1 indicator variables from caller's varlist
        note: we flip the coding so that zero indicates 'true'
              and nonzero and nonmissing indicates 'false' to
              create the same table as Copas and Loeber (1990) 
    */
    
    if ("`tab'" == "tab") tempname plbl olbl
    else               local quietly quietly
    
    tempvar prediction outcome
    mk_binary `prediction' `outcome' `touse' `varlist' `plbl' `olbl' , `tab'    
    
    
    tempname F
    `quietly' tabulate `prediction' `outcome' [`weight' `exp'] , matcell(`F')
    
    if ( (r(r)!=2) | (r(c)!=2) ) {
        display as err "table not 2 by 2"
        exit 459
    }
    
        /*
            Copas and Loeber (1990:303) present all formulas assuming 
            that e >= f; "if not, rows and columns are interchanged"
        */
    local a = `F'[1, 1]
    local b = max(`F'[1, 2], `F'[2, 1])
    local c = min(`F'[2, 1], `F'[1, 2])
    local d = `F'[2, 2]
    
    local e = `a' + `b'
    local f = `a' + `c'
    local n = `a' + `b' + `c' + `d'
    
    
    tempname rioc crit sr ll ul z p
    
        // Copas and Loeber (1990) (4)
    scalar `rioc' = (`n'*`a' - `e'*`f') / ( `f'*(`n'-`e') )
    
    scalar `crit' = -invnormal((1-`level'/100)/2)
    
    if ( mi("`smallsample'") ) {
            // Copas and Loeber (1990) (11)
        scalar `sr' = sqrt( `n'*`c'*( `n'*`f'*(`n'-`e')            ///
                    + `c'*(`n'*`e' + `e'*`f' -  2*`n'*`f' - `n'^2) ///
                    + 2*`n'*`c'^2 ) / ( (`n'-`e')^3*`f'^3 ) )
        scalar `z' = `rioc'/`sr'
        scalar `p' = 2*normal(-abs(`z'))
        scalar `ll' = `rioc' - `crit'*`sr'
        scalar `ul' = `rioc' + `crit'*`sr'
    }    
    else {
        tempname alpha beta delta
        
        scalar `alpha' = `e'/`n'
        scalar `beta'  = `f'/`n'
        
            // Copas and Loeber (1990) (20) 
        scalar `delta'  =  ln( ((`a'+.5)*(`d'+.5)) / ((`b'+.5)*(`c'+.5)) )
        
            // Copas and Loeber (1990) (21)
        scalar `sr' = sqrt( ((`e'+1)*(`e'+2)) / (`e'*(`a'+1)*(`b'+1)) ///
                    + ((`n'+1-`e')*(`n'+2-`e')) / ((`n'-`e')*(`c'+1)*(`d'+1)) )
        
        scalar `ll' = exp(`delta' - `crit'*`sr')
        scalar `ul' = exp(`delta' + `crit'*`sr')
        
            // Copas and Loeber (1990= (23)
        foreach phi in ll ul {
            scalar ``phi'' = ( 1+(``phi''-1)*(`alpha'+`beta'-2*`alpha'*`beta') ///
                           - sqrt( (1+(`alpha'+`beta')*(``phi''-1))^2          ///
                                  - 4*`alpha'*`beta'*``phi''*(``phi''-1) ) )   ///
                           / ( 2*(``phi''-1)*`beta'*(1-`alpha') )
        }
        scalar `sr' = .z
        scalar `z'  = .z
        scalar `p'  = .z
    }
    
    if ( ("`sezero'"=="sezero") | ("`asr'"=="asr") | ("`chi2'"=="chi2") ) {
            // Copas and Loeber (1990) (13)
        scalar `sr' = sqrt( (`e'*(`n'-`f')) / (`n'*`f'*(`n'-`e')) )
        scalar `z'  = `rioc'/`sr'
        scalar `p'  = 2*normal(-abs(`z'))
    }
    
    tempname table
    matrix `table' = `rioc'\ `sr'\ `z'\ `p'\ `ll'\ `ul'\ `crit'\ .\ 0
    matrix rownames `table' = b se z pvalue ll ul crit df eform
    matrix colnames `table' = RIOC // Copas and Loeber (1990) call it R
    
    
    if ( ("`detail'"=="detail") | (`"`stats'"'!="") ) {
        Detail `table' `F' `a' `b' `c' `d' `e' `f' `n' , level(`level') ///
            `smallsample' `sezero' `asr' `chi2' `detail' `stats' cii(`cii')
    }
    
    if ("`kappa'" == "kappa") {
        Kappa `table' `F' `e' `f' `n' , level(`level') ///
            `smallsample' `sezero' `asr' `chi2'
    }
    
    if ( c(noisily) ) Display `table' , nobs(`n') level(`level') ///
        `kappa' `sezero' `asr' `chi2' `cformat' `pformat' `sformat' `percent'
    
    
    if ("`matcell'" != "") matrix `matcell' = `F'
    
    return scalar rioc  =  `rioc'
    return scalar level = `level'
    return scalar N     =     `n'
    
    return local cmd      "rioc"
    
    return matrix table = `table'
end

program mk_binary
    syntax namelist(min=5 max=7) [ , TAB ]
    
    gettoken prediction namelist : namelist
    gettoken outcome    namelist : namelist
    gettoken touse      namelist : namelist
    
    foreach varname in `prediction' `outcome' {
        gettoken `varname' namelist : namelist
        quietly generate byte `varname' = (``varname'' == 0) if `touse'
        if ( mi("`tab'") ) continue
        local varlabel : variable label ``varname''
        if ( mi(`"`varlabel'"') ) local varlabel ``varname''
        label variable `varname' `"`varlabel'"'
    }
    
    if ( mi("`tab'") ) exit
    
    gettoken plbl olbl : namelist
    label define `plbl' 0 "True (+)" 1 "False (-)"
    label define `olbl' 0 "True"     1 "False"
    label values `prediction' `plbl'
    label values `outcome'    `olbl'
end

program Detail
    syntax anything , LEVEL(cilevel) ///
    [ SMALLSAMPLE SEZERO ASR CHI2 DETAIL STATS(passthru) CII(name) ]
    
    tokenize `anything'
    args table F a b c d e f n
    
    if ("`detail'" == "detail") {
            // Farrington and Loeber 1989 (2) (4) (5)
        local Total   `n' (`= `a' + `d'')
        local Chance  (`= `n'^2') (`= (2*`e'*`f'+`n'^2-`e'*`n'-`f'*`n')')
        local Maximum `n' (`= `f'+`n'-`e'')
        
        local dstats Correct:Total Correct:Chance Correct:Maximum
        local cnames `dstats'
    }
    
    if (`"`stats'"' != "") {
        // original table cells; orientation matters
        local b = `F'[1, 2]
        local c = `F'[2, 1]
        local e = `a' + `b'
        local f = `a' + `c'
        get_stats `a' `b' `c' `d' `e' `f' `n' , `stats' `detail'
        local i 0
        foreach sname of local snames {
            local dstats `dstats'  STAT`++i'
            local cnames `cnames' `sname'
        }
    }
    
    local colname : colnames `table'
    
    local ciicol (r(mean)\ r(se)\ .z\ .z\ r(lb)\ r(ub)\ .z\ .\ 0)
    if ( ("`sezero'"=="sezero") | ("`asr'"=="asr") | ("`chi2'"=="chi2") ) ///
        local  ciicol (r(mean)\ .z\ .z\ .z\ r(lb)\ r(ub)\ .z\ .\ 0)
    if ("`smallsample'" == "smallsample")   ///
        local  ciicol (r(mean)\ .z\ .z\ .z\  .z\ .z\ .z\ .\ 0)
    
    foreach s of local dstats {
        local s = substr("`s'", strpos("`s'", ":")+1, .)
        quietly cii ``s'' , level(`level') `cii'
        matrix `table' = `table', `ciicol'
    }
    
    matrix colnames `table' = `colname' `cnames'
end

program get_stats
    syntax anything , STATS(string asis) [ DETAIL ]
    
    tokenize `anything'
    args a  b  c  d  e  f  n
    args A  B  C  D  E  F  N
    args TP FP FN TN R1 C1
    args tp fp fn tn r1 c1
    
    local TPR a / f
    local TNR d / (n-f)
    local PPV a / e
    local NPV d / (n-e)
    local FPR b / (n-f) // not documented
    local FNR c / f     // not documented
    local ALL TPR TNR PPV NPV  // FPR FNR
    
    local shortcuts `ALL' FPR FNR
    
    local i 0
    while (`"`stats'"' != "") {
        gettoken sname stats : stats , parse(" :=")
        gettoken colon stats : stats , parse(" :=")
        if ( !inlist("`colon'", ":", "=") ) {
            local Usname = strupper("`sname'")
            local stats ``Usname'' `colon' `stats'
            if (`"`sname'"' == "all") continue
            if ( !`: list Usname in shortcuts' ) {
                display as err `"`sname' invalid"'
                error_stats
            }
            local sname `Usname'
        }
        gettoken nomin stats : stats , parse(" /") match(junk)
        gettoken slash stats : stats , parse(" /")
        gettoken denom stats : stats , parse(" ;") match(junk)
        
        capture noisily confirm name `sname'        
        if ( _rc ) error_stats
        if ("`slash'" != "/") {
            display as err `"'`slash'' found where / expected"'
            error_stats
        }
        
        foreach nd in nomin denom {
            // note: must substitute the two-letter macros first
            foreach x in tp fp fn tn r1 c1 a b c d e f n {
                local `nd' : subinstr local `nd '"`x'" "``x''" , all
                local x = strupper(`"`x'"')
                local `nd' : subinstr local `nd' "`x'" "``x''" , all
            }
            capture noisily local `nd' = ``nd''
            if ( _rc ) error_stats
            capture noisily numlist "``nd''" , integer range(>0)
            if ( _rc ) error_stats
        }
        
        gettoken semi stats : stats , parse(";")
        if (`"`semi'"'!=";") local stats `semi' `stats'
        
        c_local STAT`++i' `denom' `nomin'
          local snames `snames' `sname'
    }
    
    c_local snames : copy local snames
end

program error_stats
    display as err "above applies to option stats()"
    exit 198
end

program Kappa
    syntax anything [ , LEVEL(cilevel) SMALLSAMPLE SEZERO ASR CHI2 ]
    
    tokenize `anything'
    args table F e f n
        
    tempname po pe kappa
    scalar `po'    = (`F'[1, 1]+`F'[2, 2])/`n'
    scalar `pe'    = (2*`e'*`f'+`n'^2-`e'*`n'-`f'*`n')/`n'^2
    scalar `kappa' = (`po'-`pe') / (1-`pe')
    
    
    // below, we follow Fleiss et al. 1969
    tempname P p1_ p2_ p_1 p_2 tc
    matrix `P'   = `F'/`n'
    scalar `p1_' = `e'/`n'
    scalar `p2_' = (`n'-`e')/`n'
    scalar `p_1' = `f'/`n'
    scalar `p_2' = (`n'-`f')/`n'
    // note that wi_ = (w1_=p_1\ w2_=p_2); w_j = (w_1=p1_, w_2=p2_)
    
    tempname s z p crit ll ul
    
    if ( mi("`smallsample'") ) {
            // Fleiss et al. 1969 [8]
        scalar `s' = sqrt( (                                         ///
                     (`P'[1, 1]*((1-`pe')-(`p_1'+`p1_')*(1-`po'))^2  ///
                   +  `P'[1, 2]*((`p_1'+`p2_')*(1-`po'))^2           ///
                   +  `P'[2, 1]*((`p_2'+`p1_')*(1-`po'))^2           ///
                   +  `P'[2, 2]*((1-`pe')-(`p_2'+`p2_')*(1-`po'))^2) ///
                   - (`po'*`pe' - 2*`pe' + `po')^2  ) / (`n'*(1-`pe')^4) )
        scalar `z'    = `kappa'/`s'
        scalar `p'    = 2*normal(-abs(`z'))
        scalar `crit' = -invnormal((1-`level'/100)/2)
        scalar `ll'   = `kappa'-`crit'*`s'
        scalar `ul'   = `kappa'+`crit'*`s'
        if ( ("`asr'"=="asr") | ("`chi2'"=="chi2") ) ///
                matrix `tc' = (`kappa'\ .z\ .z\ .z\ `ll'\ `ul'\ `crit'\ .\ 0)
        else matrix `tc' = (`kappa'\ `s'\ `z'\ `p'\ `ll'\ `ul'\ `crit'\ .\ 0)
    }
    else matrix `tc' = (`kappa'\ .z\ .z\ .z\ .z\ .z\ .z\ .\ 0)
    
    if ("`sezero'" == "sezero") {
            // Fleiss et al. 1969 [9]
        scalar `s' = sqrt( ///
                     (`p1_'*`p_1'*(1-(`p_1'+`p1_'))^2  ///
                   +  `p1_'*`p_2'*(0-(`p_1'+`p2_'))^2  ///
                   +  `p2_'*`p_1'*(0-(`p_2'+`p1_'))^2  ///
                   +  `p2_'*`p_2'*(1-(`p_2'+`p2_'))^2) ///
                   -  `pe'^2) / ( (1-`pe')*sqrt(`n') )
        scalar `z' = `kappa'/`s'
        scalar `p' = 2*normal(-abs(`z'))
        matrix `tc'[2, 1] = `s'
        matrix `tc'[3, 1] = `z'
        matrix `tc'[4, 1] = `p'
    }
    
    local colnames : colnames `table'
    matrix `table' = `table', `tc'
    matrix colnames `table' = `colnames' Kappa
end

program Display
    syntax name(name=table) , nobs(integer) LEVEL(cilevel) ///
    [ KAPPA SEZERO ASR CHI2 CFORMAT(string) PFORMAT(string) SFORMAT(string) PERCENT ]
    
    tempname rtable
    matrix `rtable' = (`table')'
    local nrows = rowsof(`rtable')
    
    if ("`percent'" == "percent") {
        local cfd %6.2f
        matrix `rtable' = ///
            `rtable'[1..., 1]*100, J(`nrows', 3, .z), ///
            `rtable'[1..., 5..6]*100, `rtable'[1..., 7..9]   
    }
    else local cfd %9.0g
    
    set_fmt cf `cfd' `cformat'
    set_fmt pf %5.3f `pformat'
    set_fmt sf %8.2f `sformat'
    
    if      ( !`--nrows' ) local rspec &-
    else if ( !`--nrows' ) local rspec &--
    else {
        local and : display _dup(`--nrows') "&"
        local sep = cond("`kappa'"=="kappa", "-", "&")
        local rspec &-`and'`sep'-
    }
    
    local cspec & %12s | ///
        w10 `cf' & w9 `cf' o0& w8 `sf' & w6 `pf' & w11 `cf' & w10 `cf' &
    
    display as txt _newline "Relative improvement over chance" ///
                   _col(52) "Number of obs" _col(68) "= " as res %9.0g `nobs'
    
    display as txt "{hline 13}{c TT}{hline 64}"
    if ( (("`asr'"=="asr")|("`chi2'"=="chi2")) & mi("`sezero'") ) ///
    display as txt _col(14) "{c |}" _col(32) "ASR"
    display as txt _col(14) "{c |}"     ///
                   _col(21) "Coef."     ///
                   _col(29) "Std. Err." ///
                   _col(44) "z"         ///
                   _col(49) "P>|z|"     ///
                   _col(`= 61-strlen("`level'")') "[`level'% Conf. Interval]"
    display as txt "{hline 13}{c +}{hline 64}" _continue
    matlist `rtable'[., 1..6] , rspec(`rspec') cspec(`cspec') ///
        names(row) nodotz underscore
end

program set_fmt
    args f w ww
    c_local `f' `w'
    if ( mi("`ww'") ) exit
    if (fmtwidth(`"`ww'"') <= fmtwidth("`w'")) c_local `f' `ww'
    else display as txt "note: invalid `f'ormat(), using default"
end
exit

/* --------------------------------------
2.0.1 21feb2021 make binary indicators byte
                typo in a comment
2.0.0 02feb2021 changed default standard error for kappa
                confindence intervals for kappa
                new option -sezero-
                set se, z, and p to .z for -smallsample-
                -quietly- skips Display routine
                options -asr- and -chi2- no longer documented
                revised help file
1.1.0 29jan2021 support -by- (recall)
                new option -kappa-
                new option -matcell()-; not documented
                add FPR and FNR to stats(); not documented
                modified table header
                revised help file
1.0.0 26jan2021 new options -detail-, -stats()-, -{c s p}format()-
                new option -cii()-; not documented
                new option -percent-; not documented
                additional r(N), r(level), r(cmd)
                split main code into subroutines
                use locals to hold integers
                new immediate command -rioci-
                new help files
0.0.9 23jan2021 posted on Statalist