power for comparing two ROC AUCs

[Dave Airey]

I had the chance to make a usermethod yesterday for Stata's power command.

This is to compare two diagnostic's AUCs. The program creates binormal data for each diagnostic measure for the given AUCs, and then uses Stata's roccomp command to test the difference in AUCs. Power is computed over many data sets. For example,

power roc, iter(500) alpha(0.05) n(80(20)140) kappa(0.5 1.0) ///
auc1(.9) auc2(.75) rhoc(.6) rhod(.6)

Computes power over 500 datasets comparing a diagnostic with an AUC of .9 to and AUC or .75, assuming a correlation between diagnostics of .6 between both cases and controls. kappa is the ratio of controls to cases. rhoc and rhod are the correlations between diagnostics within control or disease.

The first code is saved to power_cmd_roc.ado and the second code is saved to power_cmc_roc_init.ado. Both are then put in your Stata personal directory or where Stata can find them.

// evaluator
program power_cmd_roc, rclass
    version 15.1
    
    syntax, n(integer)  /// controls
            kappa(real) /// control/disease
            rhoc(real)  ///
            rhod(real)  ///
            auc1(real)  ///
            auc2(real)  ///
            alpha(real) ///
            iter(integer)

quietly {
            
    local nd = floor(`n'/(`kappa'+1))
    local nc = floor(`n' - `nd')
    matrix C1 = (1, `rhoc' \ `rhoc', 1)
    matrix C2 = (1, `rhod' \ `rhod', 1)
    local diag1_mu = invnorm(`auc1')*sqrt(2)
    local diag2_mu = invnorm(`auc2')*sqrt(2)

    tempname memhold
    tempfile results control
    postfile `memhold' pvalue using "`results'"

    foreach i of numlist 1/`iter' {

        // make control data for diag1 and diag2
        clear
        drawnorm diag1 diag2, n(`nc') corr(C1) means(0 0) sds(1 1)
        generate refclass = 0
        save "`control'", replace

        // make disease data for diag1 and diag2
        clear
        drawnorm diag1 diag2, n(`nd') corr(C2) ///
            means(`diag1_mu' `diag2_mu') sds(1 1)
        generate refclass = 1

        // append control and disease data
        append using "`control'"

        // compare ROC AUCs for diag1 and diag2
        roccomp refclass diag1 diag2
        local p = r(p)
        post `memhold' (`p')

    }

    postclose `memhold'
    use "`results'", clear
    generate significant = pvalue < `alpha'
    summarize significant, meanonly
    tempname power
    scalar `power' = r(mean)

}

return scalar alpha = `alpha'
return scalar power = `power'
return scalar N = `n'
return scalar kappa = `kappa'
return scalar nc = `nc'
return scalar nd = `nd'
return scalar auc1 = `auc1'
return scalar auc2 = `auc2'
return scalar rhoc = `rhoc'
return scalar rhod = `rhod'

end

// initializer
program power_cmd_roc_init, sclass
    version 15.1
    sreturn clear
    sreturn local pss_numopts "kappa"
    sreturn local pss_colnames "kappa nc nd"
end

[golnoush mahmoudi]

Hi,
Thanks for the code. Could you define please what are n(80(20)140) in the code?

[Roman Mostazir]

Hi,
Thanks for the code. Could you define please what are n(80(20)140) in the code?

`n' represents the number of sample size for which power of the study will be calculated i.e., in this code, power will be calculated based on sample size 80, 100, 120 and 140.