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()
