Confidence intervals for sensitivity/specificity/ppv/npv after estat classification

[Chris Labos]

Hello all,


I'm working on a project evaluating several biomarkers to predict pre-eclampsia (high blood pressure that occcurs in pregnancy)

The basic framework of this part of the project is to add a biomarker to a clinical model of some common clinical variables.

The code for the analysis is here.
The loop is to loop through the seven biomarkers I am investigating.
The logit model has PET as a binary outcome and uses each biomarker with age, BMI, race, and clinic mean arterial pressure as additional variables in the model
Then estat to generate the sensitivity / specificity.

[
 foreach var of varlist bin_* {
 qui logit PET `var' age bmi_pre_preg race_bin clinicMAP
 estat classification, cutoff(0.5) /*can make predicted probability greater than 50%)*/

}

The problem as you might be able to deduce is that I want to generate confidence intervals for the various test characteristics (sensitivity, specificity, positive predictive value (ppv) and negative predictive value (npv) )

I eventually read some old posts, which for some reason I cannot find, which suggested bootstrapping as a possible solution. After muddling through I was able to write a simple program to bootstrap these estimates.

            *my program for bootstrapping confidence intervals
            capture program drop cutci
            program define cutci
                version 13.0
                args varname
                logit  PET `varname' age bmi_pre_preg race_bin clinicMAP  
                estat classification, cutoff(0.5)
                end

foreach var of varlist bin_* {
bootstrap sen = r(P_p1) spe = r(P_n0) ppv =r(P_1p) npv = r(P_0n) , reps(100) seed(12345): cutci `var'
}

So, all well and good so far except that the confidence intervals generated by the bootstrapping have an upper bound that surpasses 100% in some circumstances, which of course would be theoretically impossible. Below is an example of output:


------------------------------------------------------------------------------
| Observed Bootstrap Normal-based
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
sen | 78.57143 12.72385 6.18 0.000 53.63315 103.5097
spe | 98.4375 .9245184 106.47 0.000 96.62548 100.2495
ppv | 84.61538 8.047613 10.51 0.000 68.84235 100.3884
npv | 97.67442 1.354036 72.14 0.000 95.02056 100.3283
------------------------------------------------------------------------------

Here is a simplified sample of the data. Included are three of the biomarkes under investigation. I have seven biomarkers in total but have provided three for simplicity

* Example generated by -dataex-. To install: ssc install dataex
clear
input byte(PET age) double bmi_pre_preg byte race_bin double clinicMAP float(bin_Follistatin bin_GlyFibronectingmL bin_InhibinApgmL)
1 35 31.47391933 0 110.33333333333333 0 0 1
1 45 27.76709812 1                101 0 1 0
0 40 27.00513097 1  99.33333333333333 1 0 0
1 34 28.84153181 1  97.66666666666667 1 1 1
0 43 29.99671269 1  96.33333333333333 1 0 0
0 25 21.58226077 0 103.66666666666666 0 1 0
1 30 58.51490993 0 102.16666666666666 1 0 0
1 41    35.15625 1  95.33333333333333 1 1 1
1 30 26.57103476 0 102.91666666666667 1 1 1
1 42 24.40729546 1         92.8888889 1 1 1
1 48  23.0583271 0                 99 1 1 1
0 39 54.55568942 0                 98 0 0 0
0 36 19.65983459 0                 99 0 0 0
1 29 25.15589569 0              98.75 1 1 1
1 36 22.86659805 1  91.16666666666667 1 0 1
0 42 28.28282828 0  97.33333333333333 0 0 0
1 39 33.25488986 0                 97 1 0 1
1 35 27.21730295 1  89.33333333333333 0 1 0
0 39 30.50508507 1  88.83333333333333 0 0 1
0 38 23.87511478 1  88.66666666666667 0 0 1
0 41 23.26333168 0                 96 1 0 0
0 44 29.58981476 0  95.66666666666667 1 0 0
0 39  26.7755102 0  95.55555554666665 1 0 0
0 31 22.81949008 0  95.66666666666666 0 1 0
0 31 30.83653053 0  95.16666666666667 0 0 1
0 40 31.08109659 1  87.66666666666667 1 0 0
0 28 35.11682865 0                 94 1 0 0
0 32 24.27618286 1  86.66666666666666 1 0 0
0 42 39.05273621 0  93.44444445333333 1 0 0
0 42  24.1671624 1                 86 0 0 0
0 39 24.79963244 0  93.33333333333333 1 0 0
0 38 38.70975484 0  92.66666666666667 0 1 0
0 33  40.2381133 0                 92 0 0 0
0 40 27.88518739 0  92.33333333333333 0 0 0
0 31 20.80087514 0  92.33333333333333 0 0 0
0 40 24.91349481 0  92.16666666666667 1 1 0
0 41  25.2640542 1  84.66666666666667 1 0 1
0 37 23.63281255 1  84.66666666666667 0 0 0
0 29 27.43507725 0  91.66666666666667 0 0 1
0 26 32.69537789 0  91.33333333333333 1 1 0
0 33 25.99591237 0  90.66666666666667 0 0 1
0 32 23.73866213 0  90.66666666666667 0 1 1
0 38 22.96176738 0  90.58333333333333 0 0 0
0 34 32.10720958 0                 90 1 0 0
0 32 21.11268513 0  90.33333333333333 1 0 0
0 39 26.61666922 0                 90 0 1 0
0 42  24.7480492 0                 90 1 1 1
0 35 25.23692448 0                 90 0 0 0
0 44 30.24199597 0  89.66666666666667 0 1 0
0 38 27.54469581 1  82.16666666666667 0 1 0
0 33 25.71166208 1  82.16666666666667 0 1 0
0 39 27.24614824 0               89.5 0 0 0
0 35 21.67125803 0  89.66666666666667 1 1 0
0 39 24.19159307 0  89.33333333333333 1 0 1
1 36 24.12901516 0  89.33333333333333 1 0 0
0 33 48.75725132 1                 81 1 0 0
0 33   19.605192 0  89.22222223333333 0 1 1
0 30 32.02036958 0  88.83333333333333 1 0 0
0 35 25.63115586 0  88.66666666666666 0 0 0
0 34 23.56742752 0  88.66666666666667 0 1 1
0 42 23.01573179 0  88.33333333333333 1 0 0
0 34 27.46365545 0  88.22222221333334 1 0 0
0 33 23.61830085 0  88.33333333333333 1 0 0
1 37 27.87845515 0                 88 1 0 0
0 39 30.77870114 0               87.8 0 1 0
0 35 23.32341806 0                 88 0 0 0
0 35 30.11940191 0  86.99999998666667 0 0 0
0 33 32.44936521 0  86.77777778666666 0 0 0
0 36 42.43663439 0  86.33333333333333 0 1 0
0 40 29.74972106 0  86.66666666666666 0 0 0
0 39 22.89281998 0  86.66666666666667 0 0 0
0 43  20.2020202 0  86.66666666666667 0 0 1
0 37 33.28140022 0  86.16666666666667 0 0 0
0 38 33.53067931 0                 86 0 0 0
0 42 38.32001657 1  78.33333334666666 0 0 0
0 37 23.00081144 0  86.16666666666667 0 0 0
0 39 22.75830678 0  86.16666666666667 1 0 0
0 38 21.15931349 0                 86 0 0 1
0 38 22.72451072 0         85.8888889 0 0 0
0 36 27.32260031 0  85.66666666666667 0 0 0
0 36 24.88893776 0  85.22222221333334 1 0 0
0 35 29.29456582 0  84.33333333333333 0 0 0
0 41 24.09297052 0  84.44444443333333 0 1 0
0 38 32.88888889 1  76.66666666666667 1 1 0
0 39 21.03806228 0  84.33333333333333 1 1 0
0 34  21.6591398 0         84.1111111 0 0 0
end
label values race_bin race_bin
label def race_bin 0 "White or other", modify
label def race_bin 1 "Black", modify

So I have taken the matter as far as I can on my own. My question essentially is whether there is a better way to calculate the confidence intervals for sensitivity/specificity/ppv/npv?

I suppose I could just add more reps to the bootstrapping and hope the confidence intervals shrink enough to be plausible, but I feel like there must be a more sophisticated way to do this properly.

Thanks all.

Christopher Labos

[Joseph Coveney]

You could either truncate the upper confidence bound at 100% or use an alternative to the normal approximation, such as the percentile bootstrap. I show the latter below.

.ÿ
.ÿversionÿ16.1

.ÿ
.ÿclearÿ*

.ÿ
.ÿsetÿseedÿ`=strreverse("1588489")'

.ÿ
.ÿquietlyÿinputÿbyte(PETÿage)ÿdoubleÿbmi_pre_pregÿ///
>ÿÿbyteÿrace_binÿdoubleÿclinicMAPÿ///
>ÿÿfloat(bin_Follistatinÿbin_GlyFibronectingmLÿbin_InhibinApgmL)

.ÿ
.ÿprogramÿdefineÿbootem
ÿÿ1.ÿÿÿÿÿÿÿÿÿversionÿ16.1
ÿÿ2.ÿÿÿÿÿÿÿÿÿsyntaxÿvarname(numeric)
ÿÿ3.ÿ
.ÿÿÿÿÿÿÿÿÿlogitÿPETÿi.(`varlist'ÿrace_bin)ÿc.(ageÿbmi_pre_pregÿclinicMAP)
ÿÿ4.ÿÿÿÿÿÿÿÿÿestatÿclassification,ÿcutoff(0.5)
ÿÿ5.ÿend

.ÿ
.ÿforeachÿvarÿofÿvarlistÿbin_*ÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿdisplayÿinÿsmclÿasÿtextÿ_newline(2)ÿ"`var'"
ÿÿ3.ÿÿÿÿÿÿÿÿÿquietlyÿbootstrapÿsenÿ=ÿr(P_p1)ÿspeÿ=ÿr(P_n0)ÿppvÿ=r(P_1p)ÿnpvÿ=ÿr(P_0n),ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿreps(1000)ÿnodots:ÿbootemÿ`var'
ÿÿ4.ÿÿÿÿÿÿÿÿÿestatÿbootstrap,ÿpercentile
ÿÿ5.ÿ}


bin_Follistatin

BootstrapÿresultsÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿÿÿÿÿ86
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿReplicationsÿÿÿÿÿÿ=ÿÿÿÿÿÿÿÿ939

ÿÿÿÿÿÿcommand:ÿÿbootemÿbin_Follistatin
ÿÿÿÿÿÿÿÿÿÿsen:ÿÿr(P_p1)
ÿÿÿÿÿÿÿÿÿÿspe:ÿÿr(P_n0)
ÿÿÿÿÿÿÿÿÿÿppv:ÿÿr(P_1p)
ÿÿÿÿÿÿÿÿÿÿnpv:ÿÿr(P_0n)

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿObservedÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿBootstrap
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿCoef.ÿÿÿÿÿÿÿBiasÿÿÿÿStd.ÿErr.ÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿsenÿ|ÿÿÿ57.142857ÿÿÿ8.046896ÿÿÿ17.858201ÿÿÿÿÿÿÿÿÿÿ25ÿÿÿ94.11765ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿspeÿ|ÿÿÿ97.222222ÿÿ-.2795159ÿÿÿ1.9665053ÿÿÿÿ92.64706ÿÿÿÿÿÿÿÿ100ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿppvÿ|ÿÿÿÿÿÿÿÿÿÿ80ÿÿ-.2285279ÿÿÿ14.000004ÿÿÿÿÿÿÿÿÿÿ50ÿÿÿÿÿÿÿÿ100ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿnpvÿ|ÿÿÿ92.105263ÿÿÿ1.636815ÿÿÿ2.9282979ÿÿÿÿ87.65432ÿÿÿ98.68421ÿÿÿ(P)
------------------------------------------------------------------------------
(P)ÿÿÿÿpercentileÿconfidenceÿinterval
Note:ÿOneÿorÿmoreÿparametersÿcouldÿnotÿbeÿestimatedÿinÿ61ÿbootstrapÿreplicates;
ÿÿÿÿÿÿstandard-errorÿestimatesÿincludeÿonlyÿcompleteÿreplications.


bin_GlyFibronectingmL

BootstrapÿresultsÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿÿÿÿÿ86
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿReplicationsÿÿÿÿÿÿ=ÿÿÿÿÿÿÿÿ990

ÿÿÿÿÿÿcommand:ÿÿbootemÿbin_GlyFibronectingmL
ÿÿÿÿÿÿÿÿÿÿsen:ÿÿr(P_p1)
ÿÿÿÿÿÿÿÿÿÿspe:ÿÿr(P_n0)
ÿÿÿÿÿÿÿÿÿÿppv:ÿÿr(P_1p)
ÿÿÿÿÿÿÿÿÿÿnpv:ÿÿr(P_0n)

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿObservedÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿBootstrap
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿCoef.ÿÿÿÿÿÿÿBiasÿÿÿÿStd.ÿErr.ÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿsenÿ|ÿÿÿ57.142857ÿÿÿ6.509052ÿÿÿ17.342603ÿÿÿÿ23.07692ÿÿÿ93.33334ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿspeÿ|ÿÿÿ95.833333ÿÿÿ1.690524ÿÿÿ1.8745239ÿÿÿÿ93.15069ÿÿÿÿÿÿÿÿ100ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿppvÿ|ÿÿÿ72.727273ÿÿÿ11.11828ÿÿÿ12.044387ÿÿÿÿ57.14286ÿÿÿÿÿÿÿÿ100ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿnpvÿ|ÿÿÿÿÿÿÿÿÿÿ92ÿÿÿ1.323567ÿÿÿÿ3.177988ÿÿÿÿÿ86.8421ÿÿÿÿ98.7013ÿÿÿ(P)
------------------------------------------------------------------------------
(P)ÿÿÿÿpercentileÿconfidenceÿinterval
Note:ÿOneÿorÿmoreÿparametersÿcouldÿnotÿbeÿestimatedÿinÿ10ÿbootstrapÿreplicates;
ÿÿÿÿÿÿstandard-errorÿestimatesÿincludeÿonlyÿcompleteÿreplications.


bin_InhibinApgmL

BootstrapÿresultsÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿÿÿÿÿ86
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿReplicationsÿÿÿÿÿÿ=ÿÿÿÿÿÿÿÿ991

ÿÿÿÿÿÿcommand:ÿÿbootemÿbin_InhibinApgmL
ÿÿÿÿÿÿÿÿÿÿsen:ÿÿr(P_p1)
ÿÿÿÿÿÿÿÿÿÿspe:ÿÿr(P_n0)
ÿÿÿÿÿÿÿÿÿÿppv:ÿÿr(P_1p)
ÿÿÿÿÿÿÿÿÿÿnpv:ÿÿr(P_0n)

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿObservedÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿBootstrap
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿCoef.ÿÿÿÿÿÿÿBiasÿÿÿÿStd.ÿErr.ÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿsenÿ|ÿÿÿ78.571429ÿÿ-5.650158ÿÿÿ16.305765ÿÿÿÿ36.36364ÿÿÿÿÿÿÿÿ100ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿspeÿ|ÿÿÿ98.611111ÿÿÿÿ-.12293ÿÿÿ1.6867937ÿÿÿÿ94.11765ÿÿÿÿÿÿÿÿ100ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿppvÿ|ÿÿÿ91.666667ÿÿÿÿ-.53946ÿÿÿ9.2787464ÿÿÿÿÿÿÿÿÿÿ70ÿÿÿÿÿÿÿÿ100ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿnpvÿ|ÿÿÿ95.945946ÿÿ-.8940575ÿÿÿ2.9508832ÿÿÿÿÿ88.6076ÿÿÿÿÿÿÿÿ100ÿÿÿ(P)
------------------------------------------------------------------------------
(P)ÿÿÿÿpercentileÿconfidenceÿinterval
Note:ÿOneÿorÿmoreÿparametersÿcouldÿnotÿbeÿestimatedÿinÿ9ÿbootstrapÿreplicates;
ÿÿÿÿÿÿstandard-errorÿestimatesÿincludeÿonlyÿcompleteÿreplications.

.ÿ
.ÿexit

endÿofÿdo-file


.

Couple of notes:

First, you'll notice that not all replicates were successful (see the warning messages at the bottom of the -estat bootstrap- tables). With a couple of indicator variables, you might have had quasicomplete separation or complete separation in some bootstrap draws. You'll be better off with a continuous or at least an ordered-categorical representation of the measurement values of your candidate biomarkers.

Second, unless you're probing two candidate stochastic statistical methods by exact comparison, it's good practice to set the seed only once, at the top of your do-file and not change it afterward.

Are you really using Stata Release 13? Consider upgrading; there's been a lot of enhancements.

[Chris Labos]

Thanks Joseph,

My original reply was not posted. Presumably I forgot to hit enter.

Thank you for the solution. It seems to have resolved the issue.

As per yours comments, we dichotomized the biomarkers because there is a desire to establish and prove that a specific cut-off "works" in this clinical setting. I will try the continuous biomarker measurements for curiosity but I susupect that my colleagues will want to use the binary measures for practical purposes.

I am curious as to why you set the seed with the strreverse command. I am not familiar with it and so I ask out of ignorance.

As for Stata 13, the answer is simple. I would have to pay for the update out of pocket and I am too frugal these days. I always see I will buy it as a Christmas present for myself and never do. Perhaps next year.

Thank you again.

Chris

[Joseph Coveney]

. . .we dichotomized the biomarkers because there is a desire to establish and prove that a specific cut-off "works" in this clinical setting. I will try the continuous biomarker measurements for curiosity but I susupect that my colleagues will want to use the binary measures for practical purposes.

Yeah, there's been a lot written about that in the medical field that I'm sure you're aware of.

I am curious as to why you set the seed with the strreverse command. I am not familiar with it and so I ask out of ignorance.

I used to use dates as the seed, for example,

set seed `=date("2010-01-01", "YMD")’

or whatever the date happened to be. Stata’s user manual recommends against that sort of thing as it increments in a pattern. In the interest of making the seed visible to viewers on the list, I began to use the post ID number and to reverse its sequence with the strreverse() string function in order to attenuate any pattern (incrementing). There is still the first three digits which change only slowly and, despite their becoming the least-significant digits when reversed, form a bit of a pattern, which the user’s manual recommends avoiding. The best way is probably to use the user-written command setrngseed available from SSC, but I saw the reversed post ID as a reasonable compromise for the examples that I post to illustrate a suggestion.

[Chris Labos]

Yes, dichotomania. I like that term. It is a difficult problem because on the one hand you lose a lot of information by dichotomizing a variable but on the other, at some point you have to define what is and is not an abnormal test result in order to have any usefulness clinically. I split my time between clinical work and research and so I am conflicted on this issue. I suppose much depends on whether you are trying to show correlation between two variables or whether you are trying to validate a test's clinical utility. I usually do both and put at least one in the appendix. It helps me sleep at night.

Thanks for the insight regardng setting the seed. I always just randomly picked a number out of my head and was unaware of the problems regarding this.

[Chris Labos]

Hello Joseph (and everyone),

I have one more follow-up question to thist post that has come up as I finish up this project.

My last step is I wanted to generate the confidence intervals for sensitivity/specificity/ppv/npv for the base model with only the clinical variables. In other words for the model:

 logit PET i.race_bin  c.(age bmi_pre_preg clinicMAP)

It is easy enough to get the point value for those four parameters

 logit  PET   race_bin  age  bmi_pre_preg  clinicMAP
                      estat  classification,  cutoff(0.5)

But the bootstrapping command 'bootem' will not allow me to leave the varlist blank. It provides the following error for various reasons.

'varlist required
an error occurred when bootstrap executed bootem'

I have attempted to find a work-around but I have not able to figure out how to adjust the command to run the base model with only the variables: race_bin age bmi_pre_preg clinicMAP

Any help will be greatly appreciated.

Thank you very much.

Christopher Labos

[Joseph Coveney]

. . . the bootstrapping command 'bootem' will not allow me to leave the varlist blank.

. . . I have not able to figure out how to adjust the command to run the base model with only the variables: race_bin age bmi_pre_preg clinicMAP

Any help will be greatly appreciated.

Put square brackets around the argument in order to make it optional:

program define bootem
    version 16.1
    syntax [varlist(numeric)]

    logit PET i.(`varlist' race_bin) c.(age bmi_pre_preg clinicMAP)
    estat classification, cutoff(0.5)
end

[Chris Labos]

Thanks for the prompt reply Joseph.

I'm getting a no observations error when I try to run the bootstrapping.

no observations
an error occurred when bootstrap executed bootem
r(2000);

When running the code below (the only edit I made from the original is changing the version number to 13 since that is what i am working with and removing the loop I had defined before since I am running this analysis individually)

capture program drop bootem
program define bootem
    version 13.0
    syntax [varlist(numeric)]

    logit PET i.(`varlist' race_bin) c.(age bmi_pre_preg clinicMAP)
    estat classification, cutoff(0.5)
end
 
   quietly logit PET i.( race_bin) c.(age bmi_pre_preg clinicMAP)
   quietly  bootstrap  sen  =  r(P_p1)  spe  =  r(P_n0)  ppv  =r(P_1p)  npv  =  r(P_0n),  ///
                                 reps(1000)  nodots:  bootem  
                      estat  bootstrap,  percentile

I actually did manage to get the answer I needed by just putting one of the variables from the model into the code, i.e..

quietly  bootstrap  sen  =  r(P_p1)  spe  =  r(P_n0)  ppv  =r(P_1p)  npv  =  r(P_0n),  ///
                                 reps(1000)  nodots:  bootem  i.race_bin
                                estat  bootstrap,  percentile

Which seems to have run things fine. Stata apparently just ignores the same variable appearing twice.

But I am curious to know what the proper way to do this would have been so tha I can learn.

Much obliged to you for your help.

Christopher

[Joseph Coveney]

I've been stung by this before: whenever you write a program and you have the variable list as optional, Stata defaults to _all if there are no variables specified when calling the command. That is, if you specify that you want no variables in the variable list, Stata naturally assumes that you want all variables in the dataset to be included in the variable list argument.

So, add the following to the argument in the square brackets.

program define bootem
    version 16.1
    syntax [varlist(numeric default=none)]

    logit PET i.(`varlist' race_bin) c.(age bmi_pre_preg clinicMAP)
    estat classification, cutoff(0.5)
end

.ÿ
.ÿversionÿ13.0

.ÿ
.ÿclearÿ*

.ÿ
.ÿsetÿseedÿ`=strreverse("1594299")'

.ÿ
.ÿquietlyÿsetÿobsÿ250

.ÿ
.ÿforeachÿvarÿofÿnewlistÿPETÿrace_binÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿgenerateÿbyteÿ`var'ÿ=ÿruniform()ÿ<ÿ0.5
ÿÿ3.ÿ}

.ÿ
.ÿforeachÿvarÿofÿnewlistÿageÿbmi_pre_pregÿclinicMAPÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿgenerateÿdoubleÿ`var'ÿ=ÿruniform(-0.5,ÿ0.5)
ÿÿ3.ÿ}

.ÿ
.ÿprogramÿdefineÿbootem
ÿÿ1.ÿÿÿÿÿversionÿ13.0
ÿÿ2.ÿÿÿÿÿsyntaxÿ[varlist(numericÿdefault=none)]
ÿÿ3.ÿ
.ÿÿÿÿÿlogitÿPETÿi.(`varlist'ÿrace_bin)ÿc.(ageÿbmi_pre_pregÿclinicMAP)
ÿÿ4.ÿÿÿÿÿestatÿclassification,ÿcutoff(0.5)
ÿÿ5.ÿend

.ÿ
.ÿlocalÿvar

.ÿquietlyÿbootstrapÿsenÿ=ÿr(P_p1)ÿspeÿ=ÿr(P_n0)ÿppvÿ=r(P_1p)ÿnpvÿ=ÿr(P_0n),ÿ///
>ÿÿÿÿÿÿÿÿÿreps(1000)ÿnodots:ÿbootemÿ`var'

.ÿestatÿbootstrap,ÿpercentile

BootstrapÿresultsÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿÿÿÿ250
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿReplicationsÿÿÿÿÿÿ=ÿÿÿÿÿÿÿÿ998

ÿÿÿÿÿÿcommand:ÿÿbootem
ÿÿÿÿÿÿÿÿÿÿsen:ÿÿr(P_p1)
ÿÿÿÿÿÿÿÿÿÿspe:ÿÿr(P_n0)
ÿÿÿÿÿÿÿÿÿÿppv:ÿÿr(P_1p)
ÿÿÿÿÿÿÿÿÿÿnpv:ÿÿr(P_0n)

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿObservedÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿBootstrap
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿCoef.ÿÿÿÿÿÿÿBiasÿÿÿÿStd.ÿErr.ÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿsenÿ|ÿÿÿ43.220339ÿÿÿ5.297522ÿÿÿ13.146334ÿÿÿÿ19.38775ÿÿÿ71.85185ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿspeÿ|ÿÿÿ67.424242ÿÿ-.5607649ÿÿÿÿ11.53629ÿÿÿÿ41.73913ÿÿÿ88.81119ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿppvÿ|ÿÿÿ54.255319ÿÿÿÿÿ3.1477ÿÿÿ5.1882504ÿÿÿÿ47.42268ÿÿÿ65.82278ÿÿÿ(P)
ÿÿÿÿÿÿÿÿÿnpvÿ|ÿÿÿ57.051282ÿÿÿ2.372695ÿÿÿ3.4372913ÿÿÿÿ52.67176ÿÿÿ65.94595ÿÿÿ(P)
------------------------------------------------------------------------------
(P)ÿÿÿÿpercentileÿconfidenceÿinterval
Note:ÿOneÿorÿmoreÿparametersÿcouldÿnotÿbeÿestimatedÿinÿ2ÿbootstrapÿreplicates;
ÿÿÿÿÿÿstandard-errorÿestimatesÿincludeÿonlyÿcompleteÿreplications.

.ÿ
.ÿexit

endÿofÿdo-file


.

[Chris Labos]

Once again, thank you Joseph. Works like a charm.

Christopher