meaning of pg.compute_bootci() with custom function #274

FlorinAndrei · 2022-06-16T03:26:29Z

FlorinAndrei
Jun 16, 2022

This is the dataset:

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?

Answered by raphaelvallat

Jun 17, 2022

@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

View full answer

FlorinAndrei · 2022-06-16T03:36:14Z

FlorinAndrei
Jun 16, 2022
Author

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.

0 replies

raphaelvallat · 2022-06-16T19:24:46Z

raphaelvallat
Jun 16, 2022
Maintainer

@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:

Unlike your code, Pingouin replaces the pairs (Xi, Yi) and not Xi and Yi separately
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

0 replies

FlorinAndrei · 2022-06-17T01:00:57Z

FlorinAndrei
Jun 17, 2022
Author

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?).

0 replies

raphaelvallat · 2022-06-17T15:24:06Z

raphaelvallat
Jun 17, 2022
Maintainer

@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

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

meaning of pg.compute_bootci() with custom function #274

{{title}}

Replies: 4 comments

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

Select a reply

meaning of pg.compute_bootci() with custom function #274

FlorinAndrei Jun 16, 2022

Replies: 4 comments

FlorinAndrei Jun 16, 2022 Author

raphaelvallat Jun 16, 2022 Maintainer

FlorinAndrei Jun 17, 2022 Author

raphaelvallat Jun 17, 2022 Maintainer

FlorinAndrei
Jun 16, 2022

FlorinAndrei
Jun 16, 2022
Author

raphaelvallat
Jun 16, 2022
Maintainer

FlorinAndrei
Jun 17, 2022
Author

raphaelvallat
Jun 17, 2022
Maintainer