# Coverage of a Bootstrap Confidence Interval for a Change in a Binomial Proportion

How can I estimate the coverage of a bootstrap confidence interval for a change in a binomial proportion please?

For example, if the results from two tests (A and B) are:

              A
---------
B         OK    Not       Tot
-----------------------------
OK        11      6        17
Not        0     11        11
-----------------------------
T0t       11     17        28


Then I can estimate improvement (I):

Improvement (I) = pB/pA-1


using a bootstrap technique defined by the Python code below.

1000 runs
bootstrap size = 28


This gave an estimate and a 95% CI of I:

I = 55% [14%, 150%]


How can I get the coverage of that CI please?

In Python:

print('Generate the sample data')
data = pd.DataFrame({'A':[1]*11+[0]*6+[0]*11,
'B':[1]*11+[1]*6+[0]*11})
print('sample size: ',len(data))
print('')
print('A B X')
print('1 1',len(data[((data.A==1)&(data.B==1))]))
print('1 0',len(data[((data.A==1)&(data.B==0))]))
print('0 1',len(data[((data.A==0)&(data.B==1))]))
print('0 0',len(data[((data.A==0)&(data.B==0))]))
print('')

# Results
Lower = {}
Media = {}
Upper = {}

# Control Parameters
Runs_Max = 1000
Runs = range(Runs_Max)

BS = len(data)
print('bootstrap size: ',BS)

# Results
I_R = []

for R in Runs:

# Bootstrap
BooP = data.sample(BS, replace=True)

# Data
X_11 = len(BooP[((BooP.A==1)&(BooP.B==1))])
X_10 = len(BooP[((BooP.A==1)&(BooP.B==0))])
X_01 = len(BooP[((BooP.A==0)&(BooP.B==1))])
X_00 = len(BooP[((BooP.A==0)&(BooP.B==0))])

# Improvement (I) = pB/pA-1
if X_11+X_10 == 0:
I_x = 10101 # approx infinity!
else:
I_x = (X_11+X_01)/(X_11+X_10)-1

# Results
I_R.append(I_x)

# CI
Lower[R] = np.percentile(I_R,  2.5)
Media[R] = np.percentile(I_R, 50  )
Upper[R] = np.percentile(I_R, 97.5)

Low = Lower[max(list(Lower.keys()))]
Med = Media[max(list(Lower.keys()))]
Hig = Upper[max(list(Lower.keys()))]

print('I = ',Med,Low,Hig)
$$$$


I think the coverage is around 93%, slightly under the target 95%.

Figure 1, Coverage Probability (CP) against P(10)

I assessed it with a simulation.

Coverage changes with probability, so it would be good to try it with various values of the probabilities of the four possible outcomes (P11, P10, P01, P00). Unfortunately, running many values of each P would take too long. Instead, I used MultinomCI to get the Wilson interval for each P. This gave:

      Est   Low     High
P(11) 0.393 0.236   0.576
P(10) 0.000 0.000   0.121
P(01) 0.214 0.102   0.395
P(00) 0.393 0.236   0.576


I gave the name Prob to the vector [P11, P10, P01, P00]

sum (Prob) = 1


I took the range of likely values for P10 and assumed that as P10 increases, P01 decreases equally.

I assessed 4 values of Prob:

Prob = [0.3, 0.0,  0.3,  0.4]
Prob = [0.3, 0.04, 0.26, 0.4]
Prob = [0.3, 0.08, 0.22, 0.4]
Prob = [0.3, 0.12, 0.18, 0.4]


I ran it with:

n       =  28
runs    =  1000
reps    =  10000


which took 41.5 hours.

Last, I estimated the 95% CI for the coverage of the 95% CI, again using the Wilson interval.

In Python:

import numpy  as np
import pandas as pd
import time
import rpy2.robjects as ro
import statsmodels.api
import matplotlib.pyplot as plt

start = time.time()

from rpy2.robjects.packages import importr

package_name = "DescTools"

try:
pkg = importr(package_name)
except:
ro.r(f'install.packages("{package_name}")')
pkg = importr(package_name)
pkg

r_string = """CI = MultinomCI(c(11,0,6,11), conf.level=0.95, method="wilson")
"""
ro.r(r_string)
A_C = np.array(ro.r['CI'])

print(' ')
print('Estimate, CI')
print(A_C)
print(' ')

# Control Parameters
n = 28
print('n       = ',n)

nrep = 10 #10000
print('reps    = ',nrep)

runs = 10 #1000
print('runs    = ',runs)

P_10s = [0.00, 0.04, 0.08, 0.12]

d_CP = {}
d_Re = {}

for P_10 in P_10s:
pvals = [.3, P_10, (.3-P_10), .4]
print('Prob ',pvals)
print('total P = ',sum(pvals))

P_11 = pvals[0]
P_10 = pvals[1]
P_01 = pvals[2]
P_00 = pvals[3]

I_T = (P_11+P_01)/(P_11+P_10)-1
print('True I  = ',I_T)

print('Estimate the Coverage Probability using simulation')
CP = []

for it in range(nrep):

# Generate the sample data
li = np.random.multinomial(n, pvals)
data = pd.DataFrame({'A':[1]*li[0] +[1]*li[1] +[0]*li[2] +[0]*li[3],
'B':[1]*li[0] +[0]*li[1] +[1]*li[2] +[0]*li[3]})

# Results
Lower = {}
Media = {}
Upper = {}

# Results
I_R = []

for R in range(runs):

# Bootstrap size = n
BooP = data.sample(n, replace=True)

# Data
X_11 = len(BooP[((BooP.A==1)&(BooP.B==1))])
X_10 = len(BooP[((BooP.A==1)&(BooP.B==0))])
X_01 = len(BooP[((BooP.A==0)&(BooP.B==1))])
X_00 = len(BooP[((BooP.A==0)&(BooP.B==0))])

# Improvement (I) = pB/pA-1
if X_11+X_10 == 0:
I_x = 10101 # approx infinity!
else:
I_x = (X_11+X_01)/(X_11+X_10)-1

# Results
I_R.append(I_x)

# CI
Lower[R] = np.percentile(I_R,  2.5)
Media[R] = np.percentile(I_R, 50  )
Upper[R] = np.percentile(I_R, 97.5)

Low = Lower[max(list(Lower.keys()))]
Med = Media[max(list(Lower.keys()))]
Hig = Upper[max(list(Lower.keys()))]

# Check whether the interval contains the true value
if (I_T < Hig) and (I_T > Low):
CP.append(1)
else:
CP.append(0)

end = time.time()
print('time    = ',end - start)
print('CP      = ',sum(CP)/len(CP))

# results
d_Re[P_10] = CP
d_CP[P_10] = sum(CP)/len(CP)

# CI
CI_Low  = []
CI_High = []
for P_10 in d_CP.keys():
low, high = statsmodels.stats.proportion.proportion_confint(d_CP[float(P_10)]*nrep,
nrep,
alpha=1-0.95,
method='wilson')
CI_Low.append(low)
CI_High.append(high)

print('Plot')
df_G1 = pd.DataFrame({'P_10' : list(d_CP.keys()),
'CP'   : list(d_CP.values()),
'Lo'   : CI_Low,
'Hi'   : CI_High})

fig, ax1 = plt.subplots(1,1)

df_G1.plot(x='P_10', y='Hi', legend=False, ax=ax1, label='95% CI',  linewidth=5, color='k', linestyle='--')
df_G1.plot(x='P_10', y='CP', legend=False, ax=ax1, label='CP',      linewidth=5, color='k', linestyle='-')
df_G1.plot(x='P_10', y='Lo', legend=False, ax=ax1, label='95% CI',  linewidth=5, color='k', linestyle='--')

for item in ([ax1.title, ax1.xaxis.label, ax1.yaxis.label] +
ax1.get_xticklabels() + ax1.get_yticklabels()):
item.set_fontsize(22)

plt.xlabel('$$P(10)$$')
plt.ylabel('$$CP$$')

plt.xlim([0,0.12])
plt.ylim([0.9,1])

ax1.set_yticks([.9,.95,1])
plt.xticks([0.00, 0.04, 0.08, 0.12])

plt.grid(which='both', color='b')

fig = plt.gcf()
fig.set_size_inches(4,4)
plt.show()
plt.clf()
`
• Anyone know why it under covers please? Commented Feb 15, 2021 at 12:19