meaning of pg.compute_bootci() with custom function

[FlorinAndrei]

This is the dataset:

BirdPecks.csv

Prep work:

BirdPecks = pd.read_csv("BirdPecks.csv")

Nboot = 10000

pg1 = BirdPecks[BirdPecks["group"] == 1]["pecks"].to_numpy()
pg2 = BirdPecks[BirdPecks["group"] == 2]["pecks"].to_numpy()
pg1_shape = pg1.shape[0]
pg2_shape = pg2.shape[0]

sns.boxplot(data=BirdPecks, x="group", y="pecks");

I am trying to bootstrap a confidence interval for the difference in means between pg1 and pg2. Here's the "manual" version:

bd = []
for i in tqdm(range(Nboot)):
    res_pg1 = np.random.choice(pg1, size=pg1_shape, replace=True, p=None)
    res_pg2 = np.random.choice(pg2, size=pg2_shape, replace=True, p=None)
    bd.append(np.mean(res_pg1) - np.mean(res_pg2))

sns.histplot(bd)

Here's what I think is the equivalent via Pingouin:

ci, dist = pg.compute_bootci(
    x=pg1,
    y=pg2,
    func=lambda x, y: np.mean(x) - np.mean(y),
    confidence=0.9,
    n_boot=Nboot,
    paired=False,
    return_dist=True
)
print(ci)

sns.histplot(dist);

Why are the histograms so different?

The distribution returned by pg.compute_bootci() - is it not the distribution of the output of the func parameter?

[FlorinAndrei]

Another question is that the confidence interval returned by pg.compute_bootci() varies a lot if I use different options for the method parameter. I'm not sure why there is so much change. And they all differ quite a lot from:

pg.ttest(pg1, pg2, confidence=0.9)

Which gives:

               T        dof alternative     p-val  \
T-test  1.570018  29.902445   two-sided  0.126933   

                                          CI90%   cohen-d   BF10     power  
T-test  [-2.618045665165672, 67.13569272398922]  0.549579  0.851  0.323706 

In the past, I've noticed method="cper" gives the closest results to boot.ci() from R (used with type="bca"). But here the results from PG do not match at all what I get from R.

[raphaelvallat]

@FlorinAndrei please see the source code of Pingouin here:

pingouin/pingouin/effsize.py

Lines 343 to 363 in 8825a21

	# Bootstrap
	rng = np.random.RandomState(seed) # Random seed
	bootsam = rng.choice(np.arange(n), size=(n, n_boot), replace=True)
	bootstat = np.empty(n_boot)
	if y is not None:
	reference = func(x, y)
	for i in range(n_boot):
	# Note that here we use a bootstrapping procedure with replacement
	# of all the pairs (Xi, Yi). This is NOT suited for
	# hypothesis testing such as p-value estimation). Instead, for the
	# latter, one must only shuffle the Y values while keeping the X
	# values constant, i.e.:
	# >>> bootsam = rng.random_sample((n_boot, n)).argsort(axis=1)
	# >>> for i in range(n_boot):
	# >>> bootstat[i] = func(x, y[bootsam[i, :]])
	bootstat[i] = func(x[bootsam[:, i]], y[bootsam[:, i]])
	else:
	reference = func(x)
	for i in range(n_boot):
	bootstat[i] = func(x[bootsam[:, i]])

A few differences:

1. Unlike your code, Pingouin replaces the pairs (Xi, Yi) and not Xi and Yi separately
2. Because of this, n will be the size of the smallest array:

   pingouin/pingouin/effsize.py

   Line 298 in 8825a21

   n = min(n, ny)

As for the different confidence intervals, please refer to the source code, which should hopefully be self-explanatory:

pingouin/pingouin/effsize.py

Lines 365 to 393 in 8825a21

	# CONFIDENCE INTERVALS
	alpha = 1 - confidence
	dist_sorted = np.sort(bootstat)
	if method in ['norm', 'normal']:
	# Normal approximation
	za = norm.ppf(alpha / 2)
	se = np.std(bootstat, ddof=1)
	bias = np.mean(bootstat - reference)
	ll = reference - bias + se * za
	ul = reference - bias - se * za
	ci = [ll, ul]
	elif method in ['percentile', 'per']:
	# Uncorrected percentile
	pct_ll = int(n_boot * (alpha / 2))
	pct_ul = int(n_boot * (1 - alpha / 2))
	ci = [dist_sorted[pct_ll], dist_sorted[pct_ul]]
	else:
	# Corrected percentile bootstrap
	# Compute bias-correction constant z0
	z_0 = norm.ppf(np.mean(bootstat < reference) +
	np.mean(bootstat == reference) / 2)
	z_alpha = norm.ppf(alpha / 2)
	pct_ul = 100 * norm.cdf(2 * z_0 - z_alpha)
	pct_ll = 100 * norm.cdf(2 * z_0 + z_alpha)
	ll = np.percentile(bootstat, pct_ll)
	ul = np.percentile(bootstat, pct_ul)
	ci = [ll, ul]

For the unit tests of the confidence intervals against Matlab, see:

pingouin/pingouin/tests/test_effsize.py

Lines 80 to 94 in 8825a21

	# 1. bootci around a pearson correlation coefficient
	# Matlab: bootci(n_boot, {@corr, x_m, y_m}, 'type', 'norm');
	ci = compute_bootci(x_m, y_m, method='norm', seed=123)
	assert ci[0] == 0.52 and ci[1] == 1.05
	ci = compute_bootci(x_m, y_m, method='per', seed=123)
	assert ci[0] == 0.45 and ci[1] == 0.96
	ci = compute_bootci(x_m, y_m, method='cper', seed=123)
	assert ci[0] == 0.39 and ci[1] == 0.95
	# 2. Univariate function: mean
	ci_n = compute_bootci(x_m, func='mean', method='norm', seed=42)
	ci_p = compute_bootci(x_m, func='mean', method='per', seed=42)
	ci_c = compute_bootci(x_m, func='mean', method='cper', seed=42)
	assert ci_n[0] == 2.98 and ci_n[1] == 3.21
	assert ci_p[0] == 2.98 and ci_p[1] == 3.21
	assert ci_c[0] == 2.98 and round(ci_c[1], 1) == 3.2

[FlorinAndrei]

Quick question: if paired=False then why replacements are still done paired? My understanding for the generic bootstrap is that you resample each group separately - at least if no pairing is involved (with pairing I am not sure how bootstrap is supposed to work, perhaps it's supposed to be how PG is doing it?).

[raphaelvallat]

@FlorinAndrei you are right — I have just looked at the implementation of the recently-added scipy.stats.bootstrap and they indeed do resample each group separately when paired=False, which makes more sense because otherwise you're discarding data if x and y do not have the same length. My implementation was based on Matlab's bootci which only support the case where x and y have the same length:

When you use multiple data input arguments d1,...,dN, you can specify some arguments as scalar values, but all nonscalar arguments must have the same number of rows.

I'll open an issue for this. Thanks