Several methods for meta-analysis of diagnostic tests have been proposed. However, when several diagnostic tests are evaluated, the situation is problematic since the within-studies correlation needs to be taken into account.
We present an extension of the bivariate random effects meta-analysis for the log-transformed sensitivity and specificity that can be applied for two or more diagnostic tests. The advantage of this method is that a closed-form expression is derived for the calculation of the within-studies covariances. The method allows the direct calculation of sensitivity and specificity, as well as, the diagnostic odds ratio (DOR), the area under curve (AUC) and the parameters of the Summary Receiver Operator's Characteristic (SROC) curve, along with the means for a formal comparison of these quantities for different tests. There is no need for individual patient data or the simultaneous evaluation of both diagnostic tests in all studies. The method is simple and fast, it can be extended for several diagnostic tests and can be fitted in nearly all statistical packages.
Niki L. Dimou, Maria Adam and Pantelis G. Bagos A Multivariate Method for Meta-Analysis and Comparison of Diagnostic Tests. 2015, Statistics in Medicine (in press).
clear
set more off
input study tp1 fp1 fn1 tn1 tp2 fp2 fn2 tn2 n111 n110 n101 n100 n001 n000 n010 n011
1 115 53 67 63 147 10 35 106 107 8 40 27 5 58 48 5
2 30 2 23 73 31 0 22 75 23 7 8 15 0 73 2 0
3 77 16 52 52 71 2 58 66 60 17 11 41 . . . .
4 75 42 12 191 70 5 17 228 67 8 3 9 2 189 39 3
5 161 89 70 380 171 26 60 443 151 10 20 50 . . . .
6 20 32 26 25 26 1 20 56 16 4 10 16 . . . .
7 . . . . 115 17 16 73 . . . . . . . .
8 261 54 63 197 236 20 88 231 221 40 15 48 8 189 42 12
9 . . . . 69 8 18 38 . . . . . . . .
10 56 11 46 87 42 2 60 96 . . . . . . . .
11 22 2 20 80 24 3 18 79 19 3 5 15 . . . .
12 . . . . 82 13 71 301 . . . . . . . .
13 . . . . 74 11 8 130 . . . . . . . .
14 93 28 25 118 89 4 29 142 . . . . . . . .
15 32 29 3 13 25 2 10 40 25 7 0 3 0 13 27 2
16 64 18 29 73 167 8 98 88 . . . . . . . .
17 32 10 9 13 26 8 15 15 . . . . . . . .
18 . . . . 64 14 15 293 . . . . . . . .
19 36 3 5 70 38 3 0 73 . . . . . . . .
20 27 6 3 33 23 0 7 39 . . . . . . . .
21 . . . . 149 7 109 114 . . . . . . . .
22 84 41 56 90 90 2 50 129 69 15 12 44 . . . .
23 5 4 11 49 7 2 9 51 . . . . . . . .
24 196 75 99 345 190 12 105 408 152 44 38 61 6 339 69 6
25 . . . . 865 79 252 2218 . . . . . . . .
26 62 11 114 127 69 5 107 133 52 10 17 97 5 122 11 0
27 143 43 53 196 . . . . . . . . . . . .
28 73 22 29 90 68 14 35 132 58 15 10 19 1 89 12 10
29 49 23 18 359 47 7 20 375 39 10 8 10 3 356 19 4
30 35 8 51 149 . . . . . . . . . . . .
31 383 38 166 170 481 23 68 185 366 17 115 51 . . . .
32 130 8 128 113 . . . . . . . . . . . .
33 61 36 37 196 . . . . . . . . . . . .
34 70 39 36 93 . . . . . . . . . . . .
35 482 2 82 153 . . . . . . . . . . . .
36 57 25 6 111 . . . . . . . . . . . .
37 42 46 14 127 . . . . . . . . . . . .
38 57 9 33 93 . . . . . . . . . . . .
39 75 21 21 106 . . . . . . . . . . . .
40 26 1 32 29 . . . . . . . . . . . .
41 80 28 69 284 . . . . . . . . . . . .
42 64 16 27 153 . . . . . . . . . . . .
43 50 14 20 191 . . . . . . . . . . . .
44 89 3 9 39 . . . . . . . . . . . .
45 25 1 14 20 . . . . . . . . . . . .
46 8 8 0 31 . . . . . . . . . . . .
47 20 2 29 18 . . . . . . . . . . . .
48 157 287 78 1466 . . . . . . . . . . . .
49 80 50 39 45 . . . . . . . . . . . .
50 8 8 31 91 . . . . . . . . . . . .
51 143 39 63 130 . . . . . . . . . . . .
52 48 1 40 99 . . . . . . . . . . . .
53 36 6 41 313 . . . . . . . . . . . .
54 60 8 20 119 . . . . . . . . . . . .
55 18 3 31 25 . . . . . . . . . . . .
56 113 19 29 481 . . . . . . . . . . . .
57 163 10 28 140 . . . . . . . . . . . .
end
***Replace zero values***
gen zero=0
replace zero=1 if tp1==0|fp1==0|fn1==0|tn1==0|tp2==0|fp2==0|fn2==0|tn2==0
replace tp1=tp1+.5 if zero==1
replace fp1=fp1+.5 if zero==1
replace fn1=fn1+.5 if zero==1
replace tn1=tn1+.5 if zero==1
replace tp2=tp2+.5 if zero==1
replace fp2=fp2+.5 if zero==1
replace fn2=fn2+.5 if zero==1
replace tn2=tn2+.5 if zero==1
***Calculation of logit Sensitivity-(1-Specificity)***
gen b1=log(tp1/fn1)
gen b2=log(fp1/tn1)
gen b3=log(tp2/fn2)
gen b4=log(fp2/tn2)
***Calculation of TPR-FPR***
gen tpr1=exp(b1)/(1+exp(b1))
gen fpr1=exp(b2)/(1+exp(b2))
gen tpr2=exp(b3)/(1+exp(b3))
gen fpr2=exp(b4)/(1+exp(b4))
***Calculaton of the Variance***
gen V11=1/tp1+1/1/fn1
gen V22=1/tn1+1/fp1
gen V33=1/tp2+1/1/fn2
gen V44=1/tn2+1/fp2
***Extrapolation of the data***
gen total1=tp1+fp1+fn1+tn1
gen total2=tp2+fp2+fn2+tn2
gen total_cases1=tp1+fn1
gen total_cases2=tp2+fn2
gen total_controls1=fp1+tn1
gen total_controls2=fp2+tn2
gen both=1 if total1!=. & total2!=.
gen not_equal1=1 if total1<total2 & both==1
gen not_equal2=1 if total1>total2 & both==1
gen extrapolated_tp1=total_cases2*tp1/total_cases1 if not_equal1==1
gen extrapolated_fn1=total_cases2*fn1/total_cases1 if not_equal1==1
gen extrapolated_fp1=total_controls2*fp1/total_controls1 if not_equal1==1
gen extrapolated_tn1=total_controls2*tn1/total_controls1 if not_equal1==1
gen extrapolated_tp2=total_cases1*tp2/total_cases2 if not_equal2==1
gen extrapolated_fn2=total_cases1*fn2/total_cases2 if not_equal2==1
gen extrapolated_fp2=total_controls1*fp2/total_controls2 if not_equal2==1
gen extrapolated_tn2=total_controls1*tn2/total_controls2 if not_equal2==1
replace tp1=extrapolated_tp1 if not_equal1==1
replace fn1=extrapolated_fn1 if not_equal1==1
replace fp1=extrapolated_fp1 if not_equal1==1
replace tn1=extrapolated_tn1 if not_equal1==1
replace tp2=extrapolated_tp2 if not_equal2==1
replace fn2=extrapolated_fn2 if not_equal2==1
replace fp2=extrapolated_fp2 if not_equal2==1
replace tn2=extrapolated_tn2 if not_equal2==1
replace total1=total2 if not_equal1==1
replace total_cases1=total_cases2 if not_equal1==1
replace total_controls1=total_controls2 if not_equal1==1
replace total2=total1 if not_equal2==1
replace total_cases2=total_cases1 if not_equal2==1
replace total_controls2=total_controls1 if not_equal2==1
***Fitting missing cells for Yi=1***
***Calculation of common OR given Yi=1***
gen logOR_Y_1=log((n111*n100)/(n110*n101))
replace logOR_Y_1=log(((0.5+n111)*(0.5+n100))/((0.5+n110)*(0.5+n101))) if n111==0|n110==0|n101==0|n100==0
gen std_logOR_Y_1=sqrt(1/n111+1/n110+1/n101+1/n100)
replace std_logOR_Y_1=sqrt(1/(n111+0.5)+1/(n110+0.5)+1/(n101+0.5)+1/(n100+0.5)) if n111==0|n110==0|n101==0|n100==0
metan logOR_Y_1 std_logOR_Y_1, eform random nograph
rename tp1 a
rename fn1 b
rename tp2 c
rename fn2 d
gen e=.
replace e=r(ES)
***Generating values from equations (B11-14)***
cap drop x1-x4
gen x1=(-sqrt((-a*e-b-c*e+c)^2-4*a*c*(e-1)*e)+a*e+b+c*e-c)/(2*(e-1))
gen x2=(sqrt((-a*e-b-c*e+c)^2-4*a*c*(e-1)*e)+a*e-2*a-b-c*e+c)/(2*(e-1))
gen x3=(sqrt((-a*e-b-c*e+c)^2-4*a*c*(e-1)*e)-a*e-b+c*e-c)/(2*(e-1))
gen x4=(-sqrt((-a*e-b-c*e+c)^2-4*a*c*(e-1)*e)+a*e+2*b*e-b-c*e+c)/(2*(e-1))
replace n111=x1 if n111==.
replace n110=x2 if n110==.
replace n101=x3 if n101==.
replace n100=x4 if n100==.
***Fitting missing cells for Yi=0***
***Calculation of common OR given Yi=0***
gen logOR_Y_0=log((n011*n000)/(n010*n001))
replace logOR_Y_0=log(((0.5+n011)*(0.5+n000))/((0.5+n010)*(0.5+n001))) if n011==0|n010==0|n001==0|n000==0
gen std_logOR_Y_0=sqrt(1/n011+1/n010+1/n001+1/n000)
replace std_logOR_Y_0=sqrt(1/(n011+0.5)+1/(n010+0.5)+1/(n001+0.5)+1/(n000+0.5)) if n011==0|n010==0|n001==0|n000==0
metan logOR_Y_0 std_logOR_Y_0, eform random nograph
replace e=r(ES)
rename a tp1
rename b fn1
rename c tp2
rename d fn2
rename fp1 a
rename tn1 b
rename fp2 c
rename tn2 d
***Generating values from equations (B15-18)***
cap drop x1-x4
gen x1=(-sqrt((-a*e-b-c*e+c)^2-4*a*c*(e-1)*e)+a*e+b+c*e-c)/(2*(e-1))
gen x2=(sqrt((-a*e-b-c*e+c)^2-4*a*c*(e-1)*e)+a*e-2*a-b-c*e+c)/(2*(e-1))
gen x3=(sqrt((-a*e-b-c*e+c)^2-4*a*c*(e-1)*e)-a*e-b+c*e-c)/(2*(e-1))
gen x4=(-sqrt((-a*e-b-c*e+c)^2-4*a*c*(e-1)*e)+a*e+2*b*e-b-c*e+c)/(2*(e-1))
replace n011=x1 if n011==.
replace n010=x2 if n010==.
replace n001=x3 if n001==.
replace n000=x4 if n000==.
rename a fp1
rename b tn1
rename c fp2
rename d tn2
***Calculation of the Covariance***
drop zero
gen zero=1 if n111==0|n110==0|n101==0|n100==0|n011==0|n010==0|n001==0|n000==0
replace n111=n111+0.5 if zero==1
replace n110=n110+0.5 if zero==1
replace n101=n101+0.5 if zero==1
replace n100=n100+0.5 if zero==1
replace n011=n011+0.5 if zero==1
replace n010=n010+0.5 if zero==1
replace n001=n001+0.5 if zero==1
replace n000=n000+0.5 if zero==1
gen V13=(n111/((n111+n110)*(n111+n101)))-(n110/((n111+n110)*(n110+n100)))-(n101/((n101+n100)*(n111+n101)))+(n100/((n101+n100)*(n110+n100)))
gen V24=(n000/((n001+n000)*(n010+n000)))-(n001/((n001+n000)*(n011+n001)))-(n010/((n011+n010)*(n010+n000)))+(n011/((n011+n010)*(n011+n001)))
gen V14=0
gen V23=0
gen V12=0
gen V34=0
***Bivariate model fitting***
mvmeta b V, vars(b1 b2)
mvmeta b V, vars(b3 b4)
***Model fitting of Eq. 10***
mvmeta b V, vars(b1 b2 b3 b4)
test b1=b3
test b2=b4
***Calculation of Di, Si***
gen D1=b1-b2
gen D2=b3-b4
gen S1=b1+b2
gen S2=b3+b4
sum D1
local D1_mean=r(mean)
sum D2
local D2_mean=r(mean)
sum S1
local S1_mean=r(mean)
sum S2
local S2_mean=r(mean)
***Calculation of the slope of the regression of Di on Si***
gen b1_reg=((`e(Sigma11)')-(`e(Sigma22)'))/((`e(Sigma11)')+(`e(Sigma22)')+2*(`e(Sigma12)'))
gen b2_reg=((`e(Sigma33)')-(`e(Sigma44)'))/((`e(Sigma33)')+(`e(Sigma44)')+2*(`e(Sigma34)'))
***Calculation of the intercept of the regression of Di on Si***
gen a1_reg=(_b[b1]-_b[b2])-(_b[b1]+_b[b2])*((`e(Sigma11)')-(`e(Sigma22)'))/((`e(Sigma11)')+(`e(Sigma22)')+2*(`e(Sigma12)'))
gen a2_reg=(_b[b3]-_b[b4])-(_b[b3]+_b[b4])*((`e(Sigma33)')-(`e(Sigma44)'))/((`e(Sigma33)')+(`e(Sigma44)')+2*(`e(Sigma34)'))
***Construction of the plot of Di on Si***
gen D1_pred=a1_reg+b1_reg*S1
gen D2_pred=a2_reg+b2_reg*S2
twoway (sc D1 S1 [fw=total1],msize(small) msymbol(Oh)) (sc D1_pred S1,c(line) ms(i) lwidth(thin) lpattern(dot) ) (sc D2 S2 [fw=total2],msymbol(Th) msize(small) ) (sc D2_pred S2,c(line) ms(i) lwidth(thin))
***Construction of the SROC curve***
gen tse1=exp(a1_reg/(1-b1_reg))*(fpr1/(1-fpr1))^((1+b1_reg)/(1-b1_reg))/(1+(exp(a1_reg/(1-b1_reg))*(fpr1/(1-fpr1))^((1+b1_reg)/(1-b1_reg))))
gen tse2=exp(a2_reg/(1-b2_reg))*(fpr2/(1-fpr2))^((1+b2_reg)/(1-b2_reg))/(1+(exp(a2_reg/(1-b2_reg))*(fpr2/(1-fpr2))^((1+b2_reg)/(1-b2_reg))))
***Adding values (0,0) and (1,1) in our plot
local new = _N + 1
set obs `new'
replace tse1=0 in `new'
replace fpr1=0 in `new'
replace tse2=0 in `new'
replace fpr2=0 in `new'
replace total1=1 in `new'
replace total2=1 in `new'
local new = _N + 1
set obs `new'
replace tse1=1 in `new'
replace fpr1=1 in `new'
replace tse2=1 in `new'
replace fpr2=1 in `new'
replace total1=1 in `new'
replace total2=1 in `new'
gr7 tpr1 tse1 fpr1 [fw=total1], ysize(6) xsize(6) noaxis xline(0(.1)1) yline(0(.1)1) ///
tlab(0(.1)1) xlab(0(.1)1) ylab(0(.1)1) s(Oi) c(.s) ///
l1(Sensitivity) b2(1-Specificity) ti(Summary ROC Curve) key1(" ") key2(" ")
gr7 tpr2 tse2 fpr2 [fw=total2], ysize(6) xsize(6) noaxis xline(0(.1)1) yline(0(.1)1) ///
tlab(0(.1)1) xlab(0(.1)1) ylab(0(.1)1) s(Ti) c(.s[-]) ///
l1(Sensitivity) b2(1-Specificity) ti(Summary ROC Curve) key1(" ") key2(" ")
***Calculation of the AUC1,AUC2***
sum a1_reg
sum a2_reg
sum b1_reg
sum b2_reg
local a1_reg=a1_reg
local a2_reg=a2_reg
local b1_reg=b1_reg
local b2_reg=b2_reg
preserve
clear
set obs 10000
local t1=0
local t2=1
local t=`t2' - `t1'
gen x= `t1'+`t'*_n/10000
gen y=exp(`a1_reg'/(1-`b1_reg'))*(x/(1-x))^((1+`b1_reg')/(1-`b1_reg'))/(1+(exp(`a1_reg'/(1-`b1_reg'))*(x/(1-x))^((1+`b1_reg')/(1-`b1_reg'))))
integ y x
local auc1=r(integral)
dis `auc1'
restore
preserve
clear
set obs 10000
local t1=0
local t2=1
local t=`t2' - `t1'
gen x= `t1'+`t'*_n/10000
gen y=exp(`a2_reg'/(1-`b2_reg'))*(x/(1-x))^((1+`b2_reg')/(1-`b2_reg'))/(1+(exp(`a2_reg'/(1-`b2_reg'))*(x/(1-x))^((1+`b2_reg')/(1-`b2_reg'))))
integ y x
local auc2=r(integral)
dis `auc2'
restore
***Calculation of the Variance of AUC1-AUC2***
***Calculation of Var(b1)***
nlcom ((`e(Sigma11)')-(`e(Sigma22)'))/((`e(Sigma11)')+(`e(Sigma22)')+2*(`e(Sigma12)'))
matrix V=r(V)
local var_b1=V[1,1]
dis `var_b1'
***Calculation of Var(b2)***
nlcom ((`e(Sigma33)')-(`e(Sigma44)'))/((`e(Sigma33)')+(`e(Sigma44)')+2*(`e(Sigma34)'))
matrix V=r(V)
local var_b2=V[1,1]
dis `var_b2'
***Calculation of Var(a1)***
nlcom (_b[b1]-_b[b2])-(_b[b1]+_b[b2])*((`e(Sigma11)')-(`e(Sigma22)'))/((`e(Sigma11)')+(`e(Sigma22)')+2*(`e(Sigma12)'))
matrix V=r(V)
local var_a1=V[1,1]
dis `var_a1'
***Calculation of Var(a2)***
nlcom (_b[b3]-_b[b4])-(_b[b3]+_b[b4])*((`e(Sigma33)')-(`e(Sigma44)'))/((`e(Sigma33)')+(`e(Sigma44)')+2*(`e(Sigma34)'))
matrix V=r(V)
local var_a2=V[1,1]
dis `var_a2'
***Calculation of Cov(b1,b2)***
nlcom (((`e(Sigma11)')-(`e(Sigma22)'))/((`e(Sigma11)')+(`e(Sigma22)')+2*(`e(Sigma12)')))+ ///
(((`e(Sigma33)')-(`e(Sigma44)'))/((`e(Sigma33)')+(`e(Sigma44)')+2*(`e(Sigma34)')))
matrix V=r(V)
local var_b1_b2=V[1,1]
dis `var_b1_b2'
local cov_b1_b2=((`var_b1_b2')-(`var_b1')-(`var_b2'))/2
dis `cov_b1_b2'
***Calculation of Cov(a1,a2)***
nlcom ((_b[b1]-_b[b2])-(_b[b1]+_b[b2])*((`e(Sigma11)')-(`e(Sigma22)'))/((`e(Sigma11)')+(`e(Sigma22)')+2*(`e(Sigma12)')))+ ///
((_b[b3]-_b[b4])-(_b[b3]+_b[b4])*((`e(Sigma33)')-(`e(Sigma44)'))/((`e(Sigma33)')+(`e(Sigma44)')+2*(`e(Sigma34)')))
matrix V=r(V)
local var_a1_a2=V[1,1]
dis `var_a1_a2'
local cov_a1_a2=((`var_a1_a2')-(`var_a1')-(`var_a2'))/2
dis `cov_a1_a2'
***Calculation of Cov(b1,a1)***
nlcom (((`e(Sigma11)')-(`e(Sigma22)'))/((`e(Sigma11)')+(`e(Sigma22)')+2*(`e(Sigma12)')))+ ///
((_b[b1]-_b[b2])-(_b[b1]+_b[b2])*((`e(Sigma11)')-(`e(Sigma22)'))/((`e(Sigma11)')+(`e(Sigma22)')+2*(`e(Sigma12)')))
matrix V=r(V)
local var_b1_a1=V[1,1]
dis `var_b1_a1'
local cov_b1_a1=((`var_b1_a1')-(`var_b1')-(`var_a1'))/2
dis `cov_b1_a1'
***Calculation of Cov(b1,a2)***
nlcom (((`e(Sigma11)')-(`e(Sigma22)'))/((`e(Sigma11)')+(`e(Sigma22)')+2*(`e(Sigma12)')))+ ///
((_b[b3]-_b[b4])-(_b[b3]+_b[b4])*((`e(Sigma33)')-(`e(Sigma44)'))/((`e(Sigma33)')+(`e(Sigma44)')+2*(`e(Sigma34)')))
matrix V=r(V)
local var_b1_a2=V[1,1]
dis `var_b1_a2'
local cov_b1_a2=((`var_b1_a2')-(`var_b1')-(`var_a2'))/2
dis `cov_b1_a2'
***Calculation of Cov(b2,a1)***
nlcom (((`e(Sigma33)')-(`e(Sigma44)'))/((`e(Sigma33)')+(`e(Sigma44)')+2*(`e(Sigma34)')))+ ///
((_b[b1]-_b[b2])-(_b[b1]+_b[b2])*((`e(Sigma11)')-(`e(Sigma22)'))/((`e(Sigma11)')+(`e(Sigma22)')+2*(`e(Sigma12)')))
matrix V=r(V)
local var_b2_a1=V[1,1]
dis `var_b2_a1'
local cov_b2_a1=((`var_b2_a1')-(`var_b2')-(`var_a1'))/2
dis `cov_b2_a1'
***Calculation of Cov(b2,a2)***
nlcom (((`e(Sigma33)')-(`e(Sigma44)'))/((`e(Sigma33)')+(`e(Sigma44)')+2*(`e(Sigma34)')))+ ///
((_b[b3]-_b[b4])-(_b[b3]+_b[b4])*((`e(Sigma33)')-(`e(Sigma44)'))/((`e(Sigma33)')+(`e(Sigma44)')+2*(`e(Sigma34)')))
matrix V=r(V)
local var_b2_a2=V[1,1]
dis `var_b2_a2'
local cov_b2_a2=((`var_b2_a2')-(`var_b2')-(`var_a2'))/2
dis `cov_b2_a2'
***Calculation of partial derivative a1***
preserve
clear
set obs 10000
local t1=0
local t2=1
local t=`t2' - `t1'
gen x= `t1'+`t'*_n/10000
gen y=((x/(1-x))^((1+`b1_reg')/(1-`b1_reg')))/((1+((x/(1-x))^((1+`b1_reg')/(1-`b1_reg')))*exp(`a1_reg'/(1-`b1_reg')))^2)
integ y x
local integ_a1=r(integral)
restore
local deriv_a1=(1/(1-(`b1_reg')))*exp((`a1_reg')/(1-(`b1_reg')))*(`integ_a1')
dis `deriv_a1'
***Calculation of partial derivative a2****
preserve
clear
set obs 10000
local t1=0
local t2=1
local t=`t2' - `t1'
gen x= `t1'+`t'*_n/10000
gen y=((x/(1-x))^((1+`b2_reg')/(1-`b2_reg')))/((1+((x/(1-x))^((1+`b2_reg')/(1-`b2_reg')))*exp(`a2_reg'/(1-`b2_reg')))^2)
integ y x
local integ_a2=r(integral)
restore
local deriv_a2=-(1/(1-(`b2_reg')))*exp((`a2_reg')/(1-(`b2_reg')))*(`integ_a2')
dis `deriv_a2'
***Calculation of partial derivative b1****
preserve
clear
set obs 10000
local t1=0
local t2=1
local t=`t2' - `t1'
gen x= `t1'+`t'*_n/10000
gen y=(((x/(1-x))^((1+`b1_reg')/(1-`b1_reg')))*(`a1_reg'+2*log(x/(1-x))))/((1+((x/(1-x))^((1+`b1_reg')/(1-`b1_reg')))*exp(`a1_reg'/(1-`b1_reg')))^2)
integ y x
local integ_b1=r(integral)
restore
local deriv_b1=((1/(1-(`b1_reg')))^2)*exp((`a1_reg')/(1-(`b1_reg')))*(`integ_b1')
dis `deriv_b1'
***Calculation of partial derivative b2****
preserve
clear
set obs 10000
local t1=0
local t2=1
local t=`t2' - `t1'
gen x= `t1'+`t'*_n/10000
gen y=(((x/(1-x))^((1+`b2_reg')/(1-`b2_reg')))*(`a2_reg'+2*log(x/(1-x))))/((1+((x/(1-x))^((1+`b2_reg')/(1-`b2_reg')))*exp(`a2_reg'/(1-`b2_reg')))^2)
integ y x
local integ_b2=r(integral)
restore
local deriv_b2=-((1/(1-(`b2_reg')))^2)*exp((`a2_reg')/(1-(`b2_reg')))*(`integ_b2')
dis `deriv_b2'
***Retrieve Var(AUC1-AUC2)***
matrix G = (`deriv_a1',`deriv_b1',`deriv_a2',`deriv_b2')
matrix G1=G'
matrix V=((`var_a1'),(`cov_b1_a1'),(`cov_a1_a2'), (`cov_b2_a1')\(`cov_b1_a1'),(`var_b1'),(`cov_b1_a2'),(`cov_b1_b2')\(`cov_a1_a2'),(`cov_b1_a2'),(`var_a2') ,(`cov_b2_a2')\(`cov_b2_a1'),(`cov_b1_b2'),(`cov_b2_a2'),(`var_b2'))
mat v1=G*V*G1
svmat v1, names(var_auc1_auc2)
gen se_auc1_auc2=sqrt(var_auc1_auc2)
***Calculation of the test for the equality of the AUC for the two tests***
gen z=(`auc1'-`auc2')/se_auc1_auc2
gen p_value=2*(1-normal(abs(z)))