*! 1.0.0 Ariel Linden 11Nov2022

capture program drop power_cmd_tworoc
program power_cmd_tworoc, 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 = diseased; n0 = non-diseased)
			ratio(real 1)			/// ratio N0 (non-diseased) / N1 (diseased)
			Power(string)  			///
			ONESIDed				///
			HANley					/// uses Hanley & McNeil formula for variance; default is Obuchowski formula
			corr(real 0)			/// correlation between auc1 and auc2
			]						///

			gettoken auc0 rest : anything
			gettoken auc1 rest : rest
			gettoken auc2 rest : rest
			
			numlist "`anything'", min(3) max(3)
			local variable_tally : word count `anything'
			
			if `auc0' < 0.50 | `auc0' > 1.0 { 
				di as err "auc0 must contain numbers between 0.50 and 1.0"
				exit 198
			}
			
			if `auc1' < 0.50 | `auc1' > 1.0 { 
				di as err "auc1 must contain numbers between 0.50 and 1.0"
				exit 198
			}
			
			if `auc2' < 0.50 | `auc2' > 1.0 { 
				di as err "auc2 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 delta A0 v0 A1 v1 A2 v2 c2 n1 n0 n q1_0 q2_0 q1_1 q2_1 q1_2 q2_2
				scalar `zalpha' = invnorm(1-`test')
				scalar `zbeta' = invnorm(`power')
				scalar `delta' = abs(`auc1' - `auc2')
				
				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')
					scalar `A2' = invnorm(`auc2') * 1.414
					scalar `v2' = (0.0099 * exp(-`A2'^2 / 2)) * ((5 * `A2'^2 + 8) + (`A2'^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 `q1_2' = `auc2' / (2 - `auc2')
					scalar `q2_2' = 2 * (`auc2'^2) / (1 + `auc2')
					scalar `v0' = `q1_0' / `ratio' + `q2_0' - `auc0'^2 * (1 /`ratio' + 1) 
					scalar `v1' = `q1_1' / `ratio' + `q2_1' - `auc1'^2 * (1 /`ratio' + 1)
					scalar `v2' = `q1_2' / `ratio' + `q2_2' - `auc2'^2 * (1 /`ratio' + 1)
				} // end Obuchowski
				
				* sample-size formula
				scalar `n1' = ceil((`zalpha' * sqrt(`v0' + `v0' - 2 * `corr' * sqrt(`v0') * sqrt(`v0')) + `zbeta' * sqrt(`v1' + `v2' - 2 * `corr' * sqrt(`v1') * sqrt(`v2')))^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 A2 v2 power q1_0 q2_0 q1_1 q2_1 q1_2 q2_2
				
				scalar `zalpha' = invnorm(1-`test')
				scalar `delta' = abs(`auc1' - `auc2')
			
				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')
					scalar `A2' = invnorm(`auc2') * 1.414
					scalar `v2' = (0.0099 * exp(-`A2'^2 / 2)) * ((5 * `A2'^2 + 8) + (`A2'^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 `q1_2' = `auc2' / (2 - `auc2')
					scalar `q2_2' = 2 * (`auc2'^2) / (1 + `auc2')
					scalar `v0' = `q1_0' / `ratio' + `q2_0' - `auc0'^2 * (1 /`ratio' + 1) 
					scalar `v1' = `q1_1' / `ratio' + `q2_1' - `auc1'^2 * (1 /`ratio' + 1)
					scalar `v2' = `q1_2' / `ratio' + `q2_2' - `auc2'^2 * (1 /`ratio' + 1)
				} // end Obuchowski
					
			* Reorganized sample size formula to get power				
			scalar `power' = normal((sqrt(`n1' * `delta'^2) - `zalpha' * sqrt(`v0' + `v0' - 2 * `corr' * sqrt(`v0') * sqrt(`v0'))) / sqrt(`v1' + `v2' - 2 * `corr' * sqrt(`v1') * sqrt(`v2')))

			} // 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 auc2 = `auc2'
			return scalar auc1 = `auc1'
			return scalar auc0 = `auc0'
			return scalar delta = `delta'
			return scalar ratio = `ratio'
			return scalar corr = `corr'
			return scalar V2 = `v2'
			return scalar V1 = `v1'
			return scalar V0 = `v0'
		
end