*! 2.0.1 Ariel Linden 27Sep2022 // renamed "kappa" to "ratio" for consistency with official power *! 2.0.0 Ariel Linden 22Sep2022 // revised formula for computing v0; added "ordinal" option *! 1.0.1 Ariel Linden 01Jul2022 // fixed default null hypothesis, various error messages, and added return scalar for variance function *! 1.0.0 Ariel Linden 21Jun2022 capture program drop power_cmd_oneroc program power_cmd_oneroc, rclass version 11.0 /* obtain settings */ syntax anything(id="numlist"), /// [ Alpha(real 0.05) /// significance level n(string) /// total sample size n1(string) n0(string) /// group sample sizes (n1 = cases; n0 = controls) ratio(real 1) /// ratio N0 (controls) / N1 (cases) Power(string) /// ONESIDed /// HANley /// uses Hanley & McNeil formula for variance; default is Obuchowski formula ] /// gettoken auc0 rest : anything gettoken auc1 rest : rest numlist "`anything'", min(2) max(2) local variable_tally : word count `anything' if `auc1' < 0.50 | `auc1' > 1.0 { di as err "auc1 must contain numbers between 0.50 and 1.0" exit 198 } if `auc0' < 0.50 | `auc0' > 1.0 { di as err "auc0 must contain numbers between 0.50 and 1.0" exit 198 } // set default test type to "twosided" if "`onesided'" == "" { local test = `alpha'/ 2 } else local test = `alpha' ******************************* **** compute sample size ****** ******************************* if (`"`n'`n1'`n0'"' == "") { tempname zalpha zbeta v0 delta A0 A1 v1 n1 n0 n q1_0 q2_0 q1_1 q2_1 scalar `zalpha' = invnorm(1-`test') scalar `zbeta' = invnorm(`power') scalar `delta' = `auc1' - `auc0' if "`hanley'" =="" { scalar `A0' = invnorm(`auc0') * 1.414 scalar `v0' = (0.0099 * exp(-`A0'^2 / 2)) * ((5 * `A0'^2 + 8) + (`A0'^2 + 8) / `ratio') scalar `A1' = invnorm(`auc1') * 1.414 scalar `v1' = (0.0099 * exp(-`A1'^2 / 2)) * ((5 * `A1'^2 + 8) + (`A1'^2 + 8) / `ratio') } // end Hanley else { scalar `q1_0' = `auc0' / (2 - `auc0') scalar `q2_0' = 2 * (`auc0'^2) / (1 + `auc0') scalar `q1_1' = `auc1' / (2 - `auc1') scalar `q2_1' = 2 * (`auc1'^2) / (1 + `auc1') scalar `v0' = `q1_0' / `ratio' + `q2_0' - `auc0'^2 * (1 /`ratio' + 1) scalar `v1' = `q1_1' / `ratio' + `q2_1' - `auc1'^2 * (1 /`ratio' + 1) } // end Obuchowski scalar `n1' = ceil(((`zalpha' * sqrt(`v0') + `zbeta' * sqrt(`v1')) ^ 2 / `delta' ^ 2)) scalar `n0' = ceil(`n1' * `ratio') scalar `n' = `n1' + `n0' } // end sample size ************************* **** compute power ****** ************************* else if (`"`n'`n1'`n0'"' != "") { //complete sample-size specifications if (`"`n'"'=="") { tempname n if (`"`n0'"'=="") { tempname n0 scalar `n0' = ceil(`ratio'*`n1') } else if (`"`n1'"'=="") { tempname n1 scalar `n1' = ceil(`n0'/`ratio') } scalar `n' = `n1'+`n0' local ratio = `n0'/`n1' } // end n == "" else { if (`"`n1'"'!="") { tempname n0 scalar `n0' = `n' - `n1' } else if (`"`n0'"'!="") { tempname n1 scalar `n1' = `n' - `n0' } else { tempname n1 n0 scalar `n1' = ceil(`n'/(1 + `ratio')) scalar `n0' = `n'-`n1' } } // end sample-size specifications tempname zalpha A0 v0 delta A1 v1 power q1_0 q2_0 q1_1 q2_1 scalar `zalpha' = invnorm(1-`test') scalar `delta' = `auc1' - `auc0' if "`hanley'" =="" { scalar `A0' = invnorm(`auc0') * 1.414 scalar `v0' = (0.0099 * exp(-`A0'^2 / 2)) * ((5 * `A0'^2 + 8) + (`A0'^2 + 8) / `ratio') scalar `A1' = invnorm(`auc1') * 1.414 scalar `v1' = (0.0099 * exp(-`A1'^2 / 2)) * ((5 * `A1'^2 + 8) + (`A1'^2 + 8) / `ratio') } // end Hanley else { scalar `q1_0' = `auc0' / (2 - `auc0') scalar `q2_0' = 2 * (`auc0'^2) / (1 + `auc0') scalar `q1_1' = `auc1' / (2 - `auc1') scalar `q2_1' = 2 * (`auc1'^2) / (1 + `auc1') scalar `v0' = `q1_0' / `ratio' + `q2_0' - `auc0'^2 * (1 /`ratio' + 1) scalar `v1' = `q1_1' / `ratio' + `q2_1' - `auc1'^2 * (1 /`ratio' + 1) } // end Obuchowski scalar `power' = normal(((sqrt(`n1' * `delta' ^ 2) - `zalpha' * sqrt(`v0')) / sqrt(`v1'))) } // end power // saved results return scalar N = `n' return scalar N1 = `n1' return scalar N0 = `n0' return scalar alpha = `alpha' return scalar power = `power' return scalar auc1 = `auc1' return scalar auc0 = `auc0' return scalar delta = `delta' return scalar ratio = `ratio' return scalar V1 = `v1' return scalar V0 = `v0' end