Bias corrected and accelerated bootstrap of summation of two proportions

[Abdullah Algarni]

Dear Statalisters,

I want to get the bias corrected and accelerated bootstrap (BCa) 95% CI of the summation of two proportions. I know how to run BCa for one proportion but need to learn how to do it for the summation of two proportions. Here is my code to get BCa of one proportion:

cls
program drop _all

{ // BCa 95% CI of proportion of patient who died within 30-days
    { // define the "pro_bca" program
        program define pro_bca, rclass
            version 17.0
            ci proportion died30in
            return scalar died30in_pro = r(proportion)
        end
    }
    
    { // run the "pro_bca" program in bootstrap command
        bootstrap pro1 = r(died30in_pro), reps (1000) bca seed(123): pro_bca
            estat bootstrap, all
    }
}

Any help would be so appreciated .

Thank you

[Joseph Coveney]

What do you mean by "the summation of two proportions"? As in 0.95 + 0.99 = 1.94?

Oo you mean the mean (average) of two independent proportions?

Your program above calls the official Stata command ci (?!) with a variable name. This implies that the two proportions are dependent (matched by observation).

So, do you want the mean (sum?) of two dependent proportions?

[Abdullah Algarni]

What do you mean by "the summation of two proportions"? As in 0.95 + 0.99 = 1.94?

Thank you Joseph Coveney

Yes, I meant as you mentioned "As in 0.95 + 0.99 = 1.94?" with a maximum of ±2 or (±200 if multiplied by 100).
The two variables are not dependent

[Joseph Coveney]

Yes, I meant as you mentioned . . . The two variables are not dependent

Maybe something like the following? (Begin at the "Begin here" comment.)

version 17.0

clear *

// seedem
set seed 458973425

quietly set obs 200

generate byte out = rbinomial(1, 0.5)
label variable out "Outcome, binomial"

generate byte pid = mod(_n, 2)
label define Group 0 "First independent proportion" 1 "Second independent proportion"
label values pid Group
label variable pid "Proportion ID"

*
* Begin here
*
program define propEm, rclass
    version 17.0
    syntax [if]

    tempname Matcell
    tabulate out pid `if', matcell(`Matcell')
    
    return scalar sup = ///
        `Matcell'[2, 1] / (`Matcell'[1, 1] + `Matcell'[2, 1]) + ///
        `Matcell'[2, 2] / (`Matcell'[1, 2] + `Matcell'[2, 2])
end

quietly bootstrap sup = r(sup), bca /* strata(pid) */ reps(400): propEm
estat bootstrap, bca

exit

[Abdullah Algarni]

Thank you so much Joseph Coveney

I wrote a program, however error keeps appearing, see below please, I don't know what is the problem in my code.

cls
scalar drop _all
program drop _all

program define Dnri, rclass
    version 17.0
    syntax [if]
    
    quietly: count `if'
        scalar total_ct = r(N)
    
    quietly: count if died30in==1
        scalar died_ct = r(N)
        scalar died_pc = (died_ct / total_ct)
        
    quietly: count if died30in==0
        scalar alive_ct = r(N)
        scalar alive_pc = (alive_ct / total_ct)
    // Patient with event
        // Define x⁺ (No. of patient w/ event correctly reclassified)
            quietly: count if qsofa2==1 & sirs2==0 & died30in==1
                scalar x⁺ = r(N)
        

        // Define y⁺ (No. of patient w/ event incorrectly reclassified)
            quietly: count if qsofa2==0 & sirs2==1 & died30in==1
                scalar y⁺ = r(N)
        

        // Define z⁺ (Net reclassification in patient w/ event)
            scalar z⁺ = (x⁺ - y⁺)
            
        

        // Event NRI
            scalar event_nri = (z⁺ / died_ct)*100
        

        
    // Patient without event
        // Define x⁻ (No. of patient w/o event correctly reclassified)
            quietly: count if qsofa2==0 & sirs2==1 & died30in==0
                scalar x⁻ = r(N)
    

        // Define y⁻ (No. of patient w/o event incorrectly reclassified)
            quietly: count if qsofa2==1 & sirs2==0 & died30in==0
                scalar y⁻ = r(N)
    
        
        // Define z⁻ (Net reclassification in patient w/o event)
            scalar z⁻ = (x⁻ - y⁻)
    

        // Non-event NRI
            scalar noevent_nri = (z⁻ / alive_ct)*100
    

    
    // Additive NRI
        scalar aDnri = (event_nri + noevent_nri)
end

    scalar list
.         scalar list
     aDnri =  10.061369
noevent_nri =  40.374958
      z⁻ =       1206
      y⁻ =         67
      x⁻ =       1273
 event_nri = -30.313589
      z⁺ =        -87
      y⁺ =        105
      x⁺ =         18
  alive_pc =  .91233965
  alive_ct =       2987
   died_pc =  .08766035
   died_ct =        287
  total_ct =       3274


quietly bootstrap aDnri = r(aDnri), bca reps(1000): Dnri
estat bootstrap, bca
'r(aDnri)' evaluated to missing in full sample
r(322);

end of do-file

r(322);

[Joseph Coveney]

It appears that you never returned your scalar.

Maybe you intended

   // Additive NRI
        return scalar aDnri = (event_nri + noevent_nri)

[Abdullah Algarni]

It appears that you never returned your scalar.

Maybe you intended

 // Additive NRI
return scalar aDnri = (event_nri + noevent_nri)

You are right, I modify my code. However another error appeared

estat bootstrap, bca
option bca requires BCa confidence intervals to be saved by the bootstrap prefix command
r(198);

Thank you so much again

[Joseph Coveney]

I modify my code. However another error appeared

You don't show what you typed, but you might've neglected to specify the bca option to the bootstrap command immediately preceding the command that you did show.

[Abdullah Algarni]

You don't show what you typed, but you might've neglected to specify the bca option to the bootstrap command immediately preceding the command that you did show.

Sorry Joseph Coveney , here my code AND the results; "option bca requires BCa confidence intervals to be saved by the bootstrap prefix command" appeared because Jackknife can not be estimated; hence no BCa. I found this in internet: https://www.stata.com/statalist/arch.../msg01078.html. Is that true?

cls
scalar drop _all
program drop _all

program define Dnri, rclass
    version 17.0
    syntax [if]
    
    quietly: count `if'
        scalar total_ct = r(N)
    
    quietly: count if died30in==1
        scalar died_ct = r(N)
        scalar died_pc = (died_ct / total_ct)
        
    quietly: count if died30in==0
        scalar alive_ct = r(N)
        scalar alive_pc = (alive_ct / total_ct)
    // Patient with event
        // Define x⁺ (No. of patient w/ event correctly reclassified)
            quietly: count if qsofa2==1 & sirs2==0 & died30in==1
                scalar x⁺ = r(N)

        // Define y⁺ (No. of patient w/ event incorrectly reclassified)
            quietly: count if qsofa2==0 & sirs2==1 & died30in==1
                scalar y⁺ = r(N)

        // Define z⁺ (Net reclassification in patient w/ event)
            scalar z⁺ = (x⁺ - y⁺)

        // Event NRI
            scalar event_nri = (z⁺ / died_ct)*100

    // Patient without event
        // Define x⁻ (No. of patient w/o event correctly reclassified)
            quietly: count if qsofa2==0 & sirs2==1 & died30in==0
                scalar x⁻ = r(N)
    
        // Define y⁻ (No. of patient w/o event incorrectly reclassified)
            quietly: count if qsofa2==1 & sirs2==0 & died30in==0
                scalar y⁻ = r(N)
    
        // Define z⁻ (Net reclassification in patient w/o event)
            scalar z⁻ = (x⁻ - y⁻)
    
        // Non-event NRI
            scalar noevent_nri = (z⁻ / alive_ct)*100
    
    // Additive NRI
        return scalar aDnri = (event_nri + noevent_nri)
end

bootstrap aDnri = r(aDnri), bca reps(1000) : Dnri
estat bootstrap, all

. bootstrap aDnri = r(aDnri), bca reps(1000) : Dnri
(running Dnri on estimation sample)

Jackknife replications (3,274)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..................................................    50
..................................................   100
..................................................   150
..................................................   200
..................................................   250
..................................................   300
..................................................   350
..................................................   400
..................................................   450
..................................................   500
..................................................   550
..................................................   600
..................................................   650
..................................................   700
..................................................   750
..................................................   800
..................................................   850
..................................................   900
..................................................   950
.................................................. 1,000
.................................................. 1,050
.................................................. 1,100
.................................................. 1,150
.................................................. 1,200
.................................................. 1,250
.................................................. 1,300
.................................................. 1,350
.................................................. 1,400
.................................................. 1,450
.................................................. 1,500
.................................................. 1,550
.................................................. 1,600
.................................................. 1,650
.................................................. 1,700
.................................................. 1,750
.................................................. 1,800
.................................................. 1,850
.................................................. 1,900
.................................................. 1,950
.................................................. 2,000
.................................................. 2,050
.................................................. 2,100
.................................................. 2,150
.................................................. 2,200
.................................................. 2,250
.................................................. 2,300
.................................................. 2,350
.................................................. 2,400
.................................................. 2,450
.................................................. 2,500
.................................................. 2,550
.................................................. 2,600
.................................................. 2,650
.................................................. 2,700
.................................................. 2,750
.................................................. 2,800
.................................................. 2,850
.................................................. 2,900
.................................................. 2,950
.................................................. 3,000
.................................................. 3,050
.................................................. 3,100
.................................................. 3,150
.................................................. 3,200
.................................................. 3,250
........................

warning: jackknife returned missing acceleration estimates. BCa confidence intervals cannot be computed.

Bootstrap replications (1,000)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..................................................    50
..................................................   100
..................................................   150
..................................................   200
..................................................   250
..................................................   300
..................................................   350
..................................................   400
..................................................   450
..................................................   500
..................................................   550
..................................................   600
..................................................   650
..................................................   700
..................................................   750
..................................................   800
..................................................   850
..................................................   900
..................................................   950
.................................................. 1,000

Bootstrap results                                        Number of obs = 3,274
                                                         Replications  = 1,000

      Command: Dnri
        aDnri: r(aDnri)

------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       aDnri |   10.06137   3.643285     2.76   0.0058     2.920663    17.20208
------------------------------------------------------------------------------

. estat bootstrap, all

Bootstrap results                               Number of obs     =      3,274
                                                Replications      =       1000

      Command: Dnri
        aDnri: r(aDnri)

------------------------------------------------------------------------------
             |    Observed               Bootstrap
             | coefficient       Bias    std. err.  [95% conf. interval]
-------------+----------------------------------------------------------------
       aDnri |   10.061369    .042408   3.6432846    2.920663   17.20208   (N)
             |                                       3.045116   17.35165   (P)
             |                                       3.177128   17.35165  (BC)
------------------------------------------------------------------------------
Key:  N: Normal
      P: Percentile
     BC: Bias-corrected

. estat bootstrap, bca
option bca requires BCa confidence intervals to be saved by the bootstrap prefix command
r(198);

end of do-file

r(198);

[Joseph Coveney]

I found this in internet: https://www.stata.com/statalist/arch.../msg01078.html. Is that true?

Seems.

I don't have time to pore through your code in detail, but with respect to the jackknife estimates, you neglect to carry the `if' all the way through to each time you count something or other. Imposing the `if' throughout is essential for jackknife estimates, which are in turn essential for the bias-corrected, accelerated bootstrap estimate.

You might want to take advantage of the handy official Stata marksample programming aid in order to enforce the necessary leave-one-out that jackknife needs to impose on your estimation steps.

help mark

for more information.

[Abdullah Algarni]

Thank you so much Joseph Coveney , you help me solve multiple puzzles; would you allow me acknowledging you upon publication of my research?

[Joseph Coveney]

I think that it would be better acknowledge Statalist in general rather than me in particular. So, i would prefer instead that you to mention having received suggestions on Statalist for overcoming some problems during programming of the bootstrap estimator.

[Abdullah Algarni]

... but with respect to the jackknife estimates, you neglect to carry the `if' all the way through to each time you count something or other. Imposing the `if' throughout is essential for jackknife estimates, which are in turn essential for the bias-corrected, accelerated bootstrap estimate.

You might want to take advantage of the handy official Stata marksample programming aid in order to enforce the necessary leave-one-out that jackknife needs to impose on your estimation steps.

help mark

for more information.

You are right Joseph Coveney, the problem solved after adding marksample touse.

cls
scalar drop _all
program drop _all

program define Dnri, rclass
    version 17.0
    syntax [if]
    marksample touse
    quietly: count if `touse'
        scalar total_ct = r(N)
    
    quietly: count if died30in==1 & `touse'
        scalar died_ct = r(N)
        scalar died_pc = (died_ct / total_ct)
        
    quietly: count if died30in==0 & `touse'
        scalar alive_ct = r(N)
        scalar alive_pc = (alive_ct / total_ct)
    // Patient with event
        // Define x⁺ (No. of patient w/ event correctly reclassified)
            quietly: count if qsofa2==1 & sirs2==0 & died30in==1 & `touse'
                scalar x⁺ = r(N)

        // Define y⁺ (No. of patient w/ event incorrectly reclassified)
            quietly: count if qsofa2==0 & sirs2==1 & died30in==1 & `touse'
                scalar y⁺ = r(N)

        // Define z⁺ (Net reclassification in patient w/ event)
            scalar z⁺ = (x⁺ - y⁺)

        // Event NRI
            scalar event_nri = (z⁺ / died_ct)*100

    // Patient without event
        // Define x⁻ (No. of patient w/o event correctly reclassified)
            quietly: count if qsofa2==0 & sirs2==1 & died30in==0 & `touse'
                scalar x⁻ = r(N)
    
        // Define y⁻ (No. of patient w/o event incorrectly reclassified)
            quietly: count if qsofa2==1 & sirs2==0 & died30in==0 & `touse'
                scalar y⁻ = r(N)
    
        // Define z⁻ (Net reclassification in patient w/o event)
            scalar z⁻ = (x⁻ - y⁻)
    
        // Non-event NRI
            scalar noevent_nri = (z⁻ / alive_ct)*100
    
    // Additive NRI
        return scalar aDnri = (event_nri + noevent_nri)
end

bootstrap aDnri = r(aDnri), bca reps(1000) seed(123) nodots cluster(hospital): Dnri
estat bootstrap, all

returns

warning: Dnri does not set e(sample), so no observations will be excluded
         from the resampling because of missing values or other reasons. To
         exclude observations, press Break, save the data, drop any
         observations that are to be excluded, and rerun bootstrap.

Bootstrap results                                        Number of obs = 3,274
                                                         Replications  = 1,000

      Command: Dnri
        aDnri: r(aDnri)

                               (Replications based on 21 clusters in hospital)
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       aDnri |   10.06137   4.470131     2.25   0.0244     1.300073    18.82267
------------------------------------------------------------------------------

. estat bootstrap, all

Bootstrap results                               Number of obs     =      3,274
                                                Replications      =       1000

      Command: Dnri
        aDnri: r(aDnri)

                               (Replications based on 21 clusters in hospital)
------------------------------------------------------------------------------
             |    Observed               Bootstrap
             | coefficient       Bias    std. err.  [95% conf. interval]
-------------+----------------------------------------------------------------
       aDnri |   10.061369   .4567508   4.4701313    1.300073   18.82267   (N)
             |                                       1.676063   19.49649   (P)
             |                                       .0806428   17.72687  (BC)
             |                                      -.7476254   17.17115 (BCa)
------------------------------------------------------------------------------
Key:   N: Normal
       P: Percentile
      BC: Bias-corrected
     BCa: Bias-corrected and accelerated