Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • power for comparing two ROC AUCs

    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.

    Code:
    // 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


    Code:
    // 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

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

    Comment


    • #3
      Originally posted by golnoush mahmoudi View Post
      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.
      Roman

      Comment

      Working...
      X