U Qv f=x @spddlZddlZddlmZddlmZmZdddddgZdddZ dddZ dddZ dddZ dddZ dS)Npearsonr) _check_eftype remove_na compute_escicompute_bootciconvert_effsizecompute_effsizecompute_effsize_from_tFcohenffffff? two-sidedcCsddlm}m} |dks t|dks0td|dk r@|dk sDtt|tsRtd|krfdkspntd|d kr8t|} dt |d } |d krt | d|d } t | | | | | | g} nR|d kr| |} t | | | tj g} n$| |} t tj | | | g} t| }n|dksH|rpt d||d d |} |d}n6t |||||d d ||} ||d }t | d|d |} t || | || | g}t||S)a{Parametric confidence intervals around a Cohen d or a correlation coefficient. Parameters ---------- stat : float Original effect size. Must be either a correlation coefficient or a Cohen-type effect size (Cohen d or Hedges g). nx, ny : int Length of vector x and y. paired : bool Indicates if the effect size was estimated from a paired sample. This is only relevant for cohen or hedges effect size. eftype : string Effect size type. Must be "r" (correlation) or "cohen" (Cohen d or Hedges g). confidence : float Confidence level (0.95 = 95%) decimals : int Number of rounded decimals. alternative : string Defines the alternative hypothesis, or tail for the correlation coefficient. Must be one of "two-sided" (default), "greater" or "less". This parameter only has an effect if ``eftype`` is "r". Returns ------- ci : array Desired converted effect size Notes ----- To compute the parametric confidence interval around a **Pearson r correlation** coefficient, one must first apply a Fisher's r-to-z transformation: .. math:: z = 0.5 \cdot \ln \frac{1 + r}{1 - r} = \text{arctanh}(r) and compute the standard error: .. math:: \text{SE} = \frac{1}{\sqrt{n - 3}} where :math:`n` is the sample size. The lower and upper confidence intervals - *in z-space* - are then given by: .. math:: \text{ci}_z = z \pm \text{crit} \cdot \text{SE} where :math:`\text{crit}` is the critical value of the normal distribution corresponding to the desired confidence level (e.g. 1.96 in case of a 95% confidence interval). These confidence intervals can then be easily converted back to *r-space*: .. math:: \text{ci}_r = \frac{\exp(2 \cdot \text{ci}_z) - 1} {\exp(2 \cdot \text{ci}_z) + 1} = \text{tanh}(\text{ci}_z) A formula for calculating the confidence interval for a **Cohen d effect size** is given by Hedges and Olkin (1985, p86). If the effect size estimate from the sample is :math:`d`, then it follows a T distribution with standard error: .. math:: \text{SE} = \sqrt{\frac{n_x + n_y}{n_x \cdot n_y} + \frac{d^2}{2 (n_x + n_y)}} where :math:`n_x` and :math:`n_y` are the sample sizes of the two groups. In one-sample test or paired test, this becomes: .. math:: \text{SE} = \sqrt{\frac{1}{n_x} + \frac{d^2}{2 n_x}} The lower and upper confidence intervals are then given by: .. math:: \text{ci}_d = d \pm \text{crit} \cdot \text{SE} where :math:`\text{crit}` is the critical value of the T distribution corresponding to the desired confidence level. References ---------- * https://en.wikipedia.org/wiki/Fisher_transformation * Hedges, L., and Ingram Olkin. "Statistical models for meta-analysis." (1985). * http://www.leeds.ac.uk/educol/documents/00002182.htm * https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5133225/ Examples -------- 1. Confidence interval of a Pearson correlation coefficient >>> import pingouin as pg >>> x = [3, 4, 6, 7, 5, 6, 7, 3, 5, 4, 2] >>> y = [4, 6, 6, 7, 6, 5, 5, 2, 3, 4, 1] >>> nx, ny = len(x), len(y) >>> stat = pg.compute_effsize(x, y, eftype='r') >>> ci = pg.compute_esci(stat=stat, nx=nx, ny=ny, eftype='r') >>> print(round(stat, 4), ci) 0.7468 [0.27 0.93] 2. Confidence interval of a Cohen d >>> stat = pg.compute_effsize(x, y, eftype='cohen') >>> ci = pg.compute_esci(stat, nx=nx, ny=ny, eftype='cohen', decimals=3) >>> print(round(stat, 4), ci) 0.1538 [-0.737 1.045] r)normt)rpearsonspearmanr dghedges)rgreaterZlesszFAlternative must be one of 'two-sided' (default), 'greater' or 'less'.N#confidence must be between 0 and 1.)rrrrr r) scipy.statsrrlowerAssertionError isinstancefloatnpZarctanhsqrtabsppfarrayinftanhround)statnxnypairedeftype confidencedecimals alternativerrzsecritZci_zcidofr54/tmp/pip-unpacked-wheel-2te3nxqf/pingouin/effsize.pyrs8w        * cperc sddlm} m} ddlm} t|ts*td|kr>dksHntd|dksTt| |svt|tsv| |svtdd} t |}|j }|j dkstd |dkstd |d k rt |}|j }|j dkstd |dkstd r||kstdt|trd||dkr0s&tddd}n|dkrbddlm sTtdfdd}nx|dkrddlmfdd}nP|dkrd} dd}n8|dkrd} dd}n |dkrd} d d}ntd!t |}t j|}|jt |||fdd"}|d k r|||}rdt|D]4}|||d d |f||d d |f||<q,nX|jt |||fdd"}t|D]4}|||d d |f||d d |f||<qnH||}| r|||}n,t|D]"}|||d d |f||<qd|d#}|d$krj| |}t j|dd%}t ||}t ||||||||g}nL|d&krd't |d|g}t ||}ndd(lm}|||d|d)}t ||}| r||fS|Sd S)*uBootstrapped confidence intervals of univariate and bivariate functions. Parameters ---------- x : 1D-array or list First sample. Required for both bivariate and univariate functions. y : 1D-array, list, or None Second sample. Required only for bivariate functions. func : str or custom function Function to compute the bootstrapped statistic. Accepted string values are: * ``'pearson'``: Pearson correlation (bivariate, paired x and y) * ``'spearman'``: Spearman correlation (bivariate, paired x and y) * ``'cohen'``: Cohen d effect size (bivariate, paired or unpaired x and y) * ``'hedges'``: Hedges g effect size (bivariate, paired or unpaired x and y) * ``'mean'``: Mean (univariate = only x) * ``'std'``: Standard deviation (univariate) * ``'var'``: Variance (univariate) method : str Method to compute the confidence intervals (see Notes): * ``'cper'``: Bias-corrected percentile method (default) * ``'norm'``: Normal approximation with bootstrapped bias and standard error * ``'per'``: Simple percentile paired : boolean Indicates whether x and y are paired or not. For example, for correlation functions or paired T-test, x and y are assumed to be paired. Pingouin will resample the pairs (x_i, y_i) when paired=True, and resample x and y separately when paired=False. If paired=True, x and y must have the same number of elements. confidence : float Confidence level (0.95 = 95%) n_boot : int Number of bootstrap iterations. The higher, the better, the slower. decimals : int Number of rounded decimals. seed : int or None Random seed for generating bootstrap samples. return_dist : boolean If True, return the confidence intervals and the bootstrapped distribution (e.g. for plotting purposes). Returns ------- ci : array Bootstrapped confidence intervals. Notes ----- Results have been tested against the `bootci `_ Matlab function. Since version 1.7, SciPy also includes a built-in bootstrap function :py:func:`scipy.stats.bootstrap`. The SciPy implementation has two advantages over Pingouin: it is faster when using ``vectorized=True``, and it supports the bias-corrected and accelerated (BCa) confidence intervals for univariate functions. However, unlike Pingouin, it does not return the bootstrap distribution. The percentile bootstrap method (``per``) is defined as the :math:`100 \times \frac{\alpha}{2}` and :math:`100 \times \frac{1 - \alpha}{2}` percentiles of the distribution of :math:`\theta` estimates obtained from resampling, where :math:`\alpha` is the level of significance (1 - confidence, default = 0.05 for 95% CIs). The bias-corrected percentile method (``cper``) corrects for bias of the bootstrap distribution. This method is different from the BCa method — the default in Matlab and SciPy — which corrects for both bias and skewness of the bootstrap distribution using jackknife resampling. The normal approximation method (``norm``) calculates the confidence intervals with the standard normal distribution using bootstrapped bias and standard error. References ---------- * DiCiccio, T. J., & Efron, B. (1996). Bootstrap confidence intervals. Statistical science, 189-212. * Davison, A. C., & Hinkley, D. V. (1997). Bootstrap methods and their application (Vol. 1). Cambridge university press. * Jung, Lee, Gupta, & Cho (2019). Comparison of bootstrap confidence interval methods for GSCA using a Monte Carlo simulation. Frontiers in psychology, 10, 2215. Examples -------- 1. Bootstrapped 95% confidence interval of a Pearson correlation >>> import pingouin as pg >>> import numpy as np >>> rng = np.random.default_rng(42) >>> x = rng.normal(loc=4, scale=2, size=100) >>> y = rng.normal(loc=3, scale=1, size=100) >>> stat = np.corrcoef(x, y)[0][1] >>> ci = pg.compute_bootci(x, y, func='pearson', paired=True, seed=42, decimals=4) >>> print(round(stat, 4), ci) 0.0945 [-0.098 0.2738] Let's compare to SciPy's built-in bootstrap function >>> from scipy.stats import bootstrap >>> bt_scipy = bootstrap( ... data=(x, y), statistic=lambda x, y: np.corrcoef(x, y)[0][1], ... method="basic", vectorized=False, n_resamples=2000, paired=True, random_state=42) >>> np.round(bt_scipy.confidence_interval, 4) array([-0.0952, 0.2883]) 2. Bootstrapped 95% confidence interval of a Cohen d >>> stat = pg.compute_effsize(x, y, eftype='cohen') >>> ci = pg.compute_bootci(x, y, func='cohen', seed=42, decimals=3) >>> print(round(stat, 4), ci) 0.7009 [0.403 1.009] 3. Bootstrapped confidence interval of a standard deviation (univariate) >>> import numpy as np >>> stat = np.std(x, ddof=1) >>> ci = pg.compute_bootci(x, func='std', seed=123) >>> print(round(stat, 4), ci) 1.5534 [1.38 1.8 ] Compare to SciPy's built-in bootstrap function, which returns the bias-corrected and accelerated CIs (see Notes). >>> def std(x, axis): ... return np.std(x, ddof=1, axis=axis) >>> bt_scipy = bootstrap(data=(x, ), statistic=std, n_resamples=2000, random_state=123) >>> np.round(bt_scipy.confidence_interval, 2) array([1.39, 1.81]) Changing the confidence intervals type in Pingouin >>> pg.compute_bootci(x, func='std', seed=123, method="norm") array([1.37, 1.76]) >>> pg.compute_bootci(x, func='std', seed=123, method="percentile") array([1.35, 1.75]) 4. Bootstrapped confidence interval using a custom univariate function >>> from scipy.stats import skew >>> round(skew(x), 4), pg.compute_bootci(x, func=skew, n_boot=10000, seed=123) (-0.137, array([-0.55, 0.32])) 5. Bootstrapped confidence interval using a custom bivariate function. Here, x and y are not paired and can therefore have different sizes. >>> def mean_diff(x, y): ... return np.mean(x) - np.mean(y) >>> y2 = rng.normal(loc=3, scale=1, size=200) # y2 has 200 samples, x has 100 >>> ci = pg.compute_bootci(x, y2, func=mean_diff, n_boot=10000, seed=123) >>> print(round(mean_diff(x, y2), 2), ci) 0.88 [0.54 1.21] We can also get the bootstrapped distribution >>> ci, bt = pg.compute_bootci(x, y2, func=mean_diff, n_boot=10000, return_dist=True, seed=9) >>> print(f"The bootstrap distribution has {bt.size} samples. The mean and standard " ... f"{bt.mean():.4f} ± {bt.std():.4f}") The bootstrap distribution has 10000 samples. The mean and standard 0.8807 ± 0.1704 r) isfunction isroutinerrr)rnormal percentileperZ cpercentiler7zyfunc must be a function (e.g. np.mean, custom function) or a string (e.g. 'pearson'). See documentation for more details.Fzx must be one-dimensional.z"x must have more than one element.Nzy must be one-dimensional.z"y must have more than one element.z?x and y must have the same number of elements when paired=True.z%srz5Paired should be True if using correlation functions.cSst||dSNrrxyr5r5r6funcszcompute_bootci..funcr spearmanrcs||dSr?r5r@rDr5r6rCs)r r)r cs||dS)N)r+r,r5r@)r func_strr+r5r6rCsmeanTcSstj|ddS)Nr)axis)r rGrAr5r5r6rCsstdcSstj|dddSNrr)ddofrH)r rJrIr5r5r6rCsvarcSstj|dddSrK)r rMrIr5r5r6rCszFunction string not recognized.)sizereplacer )rr<rL)r=r>d)_bias_corrected_ci)alpha)inspectr9r:rrrrrstrr asarrayrNndimrEZpingouin.effsizer ValueErroremptyrandomZ default_rngchoiceZarangeranger#rJrGr$r=Zpingouin.regressionrRr')rArBrCmethodr+r-Zn_bootr.seedZ return_distr9r:rZ vectorizabler)r*ZbootstatrngZboot_x referenceiZboot_yrSZzar1Zbiasr3intervalrRr5)r rFr+rEr6rs,                      4 4     (    c Cs|}|}||fD] }t|sd|d}t|q|dkrJtd||ksZ|dkr^|S|dkrd|td|dn|} |d kr| S|d krtd d ||fDr| dd d||dStd| Sn|dkr8tdd ||fDr||dd||||} nd} | t| d| S|dkr^| ddd| ddS|dkrt| tj td S|dkrtdn ddl m } | | tdSdS)a Conversion between effect sizes. Parameters ---------- ef : float Original effect size. input_type : string Effect size type of ef. Must be ``'cohen'`` or ``'pointbiserialr'``. output_type : string Desired effect size type. Available methods are: * ``'cohen'``: Unbiased Cohen d * ``'hedges'``: Hedges g * ``'pointbiserialr'``: Point-biserial correlation * ``'eta-square'``: Eta-square * ``'odds-ratio'``: Odds ratio * ``'AUC'``: Area Under the Curve * ``'none'``: pass-through (return ``ef``) nx, ny : int, optional Length of vector x and y. Required to convert to Hedges g. Returns ------- ef : float Desired converted effect size See Also -------- compute_effsize : Calculate effect size between two set of observations. compute_effsize_from_t : Convert a T-statistic to an effect size. Notes ----- The formula to convert from a`point-biserial correlation `_ **r** to **d** is given in [1]_: .. math:: d = \frac{2r_{pb}}{\sqrt{1 - r_{pb}^2}} The formula to convert **d** to a point-biserial correlation **r** is given in [2]_: .. math:: r_{pb} = \frac{d}{\sqrt{d^2 + \frac{(n_x + n_y)^2 - 2(n_x + n_y)} {n_xn_y}}} The formula to convert **d** to :math:`\eta^2` is given in [3]_: .. math:: \eta^2 = \frac{(0.5 d)^2}{1 + (0.5 d)^2} The formula to convert **d** to an odds-ratio is given in [4]_: .. math:: \text{OR} = \exp (\frac{d \pi}{\sqrt{3}}) The formula to convert **d** to area under the curve is given in [5]_: .. math:: \text{AUC} = \mathcal{N}_{cdf}(\frac{d}{\sqrt{2}}) References ---------- .. [1] Rosenthal, Robert. "Parametric measures of effect size." The handbook of research synthesis 621 (1994): 231-244. .. [2] McGrath, Robert E., and Gregory J. Meyer. "When effect sizes disagree: the case of r and d." Psychological methods 11.4 (2006): 386. .. [3] Cohen, Jacob. "Statistical power analysis for the behavioral sciences. 2nd." (1988). .. [4] Borenstein, Michael, et al. "Effect sizes for continuous data." The handbook of research synthesis and meta-analysis 2 (2009): 221-235. .. [5] Ruscio, John. "A probability-based measure of effect size: Robustness to base rates and other factors." Psychological methods 1 3.1 (2008): 19. Examples -------- 1. Convert from Cohen d to eta-square >>> import pingouin as pg >>> d = .45 >>> eta = pg.convert_effsize(d, 'cohen', 'eta-square') >>> print(eta) 0.048185603807257595 2. Convert from Cohen d to Hegdes g (requires the sample sizes of each group) >>> pg.convert_effsize(.45, 'cohen', 'hedges', nx=10, ny=10) 0.4309859154929578 3. Convert a point-biserial correlation to Cohen d >>> rpb = 0.40 >>> d = pg.convert_effsize(rpb, 'pointbiserialr', 'cohen') >>> print(d) 0.8728715609439696 4. Reverse operation: convert Cohen d to a point-biserial correlation >>> pg.convert_effsize(d, 'cohen', 'pointbiserialr') 0.4000000000000001 Could not interpret input '')pointbiserialrr z.Input type must be 'cohen' or 'pointbiserialr'nonerer rr rcss|]}|dk VqdSNr5.0vr5r5r6 isz"convert_effsize..r zUYou need to pass nx and ny arguments to compute Hedges g. Returning Cohen's d insteadcss|]}|dk VqdSrgr5rhr5r5r6rktsz eta-squarez odds-ratiorzmUsing effect size 'r' in `pingouin.convert_effsize` has been deprecated. Please use 'pointbiserialr' instead.rr;N) rrrXr r!allwarningswarnexppirrZcdf) ZefZ input_type output_typer)r*itZotinperrrarr5r5r6rsFj   & "    c Cst|sd|d}t|t|}t|}|j|jkrN|rNtdd}t|||d\}}|j|j}}|dkr|||j dd}|S| dkrt ||\}} |S| d kr|d d d f|} t | d kd | d kS|sD||d } t |d|jdd|d|jdd| } ||| }n2||t |jdd|jddd }t|d|||dSd S)aCalculate effect size between two set of observations. Parameters ---------- x : np.array or list First set of observations. y : np.array or list Second set of observations. paired : boolean If True, uses Cohen d-avg formula to correct for repeated measurements (see Notes). eftype : string Desired output effect size. Available methods are: * ``'none'``: no effect size * ``'cohen'``: Unbiased Cohen d * ``'hedges'``: Hedges g * ``'r'``: Pearson correlation coefficient * ``'pointbiserialr'``: Point-biserial correlation * ``'eta-square'``: Eta-square * ``'odds-ratio'``: Odds ratio * ``'AUC'``: Area Under the Curve * ``'CLES'``: Common Language Effect Size Returns ------- ef : float Effect size See Also -------- convert_effsize : Conversion between effect sizes. compute_effsize_from_t : Convert a T-statistic to an effect size. Notes ----- Missing values are automatically removed from the data. If ``x`` and ``y`` are paired, the entire row is removed. If ``x`` and ``y`` are independent, the Cohen :math:`d` is: .. math:: d = \frac{\overline{X} - \overline{Y}} {\sqrt{\frac{(n_{1} - 1)\sigma_{1}^{2} + (n_{2} - 1) \sigma_{2}^{2}}{n1 + n2 - 2}}} If ``x`` and ``y`` are paired, the Cohen :math:`d_{avg}` is computed: .. math:: d_{avg} = \frac{\overline{X} - \overline{Y}} {\sqrt{\frac{(\sigma_1^2 + \sigma_2^2)}{2}}} The Cohen's d is a biased estimate of the population effect size, especially for small samples (n < 20). It is often preferable to use the corrected Hedges :math:`g` instead: .. math:: g = d \times (1 - \frac{3}{4(n_1 + n_2) - 9}) The common language effect size is the proportion of pairs where ``x`` is higher than ``y`` (calculated with a brute-force approach where each observation of ``x`` is paired to each observation of ``y``, see :py:func:`pingouin.wilcoxon` for more details): .. math:: \text{CL} = P(X > Y) + .5 \times P(X = Y) For other effect sizes, Pingouin will first calculate a Cohen :math:`d` and then use the :py:func:`pingouin.convert_effsize` to convert to the desired effect size. References ---------- * Lakens, D., 2013. Calculating and reporting effect sizes to facilitate cumulative science: a practical primer for t-tests and ANOVAs. Front. Psychol. 4, 863. https://doi.org/10.3389/fpsyg.2013.00863 * Cumming, Geoff. Understanding the new statistics: Effect sizes, confidence intervals, and meta-analysis. Routledge, 2013. * https://osf.io/vbdah/ Examples -------- 1. Cohen d from two independent samples. >>> import numpy as np >>> import pingouin as pg >>> x = [1, 2, 3, 4] >>> y = [3, 4, 5, 6, 7] >>> pg.compute_effsize(x, y, paired=False, eftype='cohen') -1.707825127659933 The sign of the Cohen d will be opposite if we reverse the order of ``x`` and ``y``: >>> pg.compute_effsize(y, x, paired=False, eftype='cohen') 1.707825127659933 2. Hedges g from two paired samples. >>> x = [1, 2, 3, 4, 5, 6, 7] >>> y = [1, 3, 5, 7, 9, 11, 13] >>> pg.compute_effsize(x, y, paired=True, eftype='hedges') -0.8222477210374874 3. Common Language Effect Size. >>> pg.compute_effsize(x, y, eftype='cles') 0.2857142857142857 In other words, there are ~29% of pairs where ``x`` is higher than ``y``, which means that there are ~71% of pairs where ``x`` is *lower* than ``y``. This can be easily verified by changing the order of ``x`` and ``y``: >>> pg.compute_effsize(y, x, eftype='cles') 0.7142857142857143 rcrdz9x and y have unequal sizes. Switching to paired == False.F)r+rrPrZclesNrg?r r r)r*)rrXr rVrNrorprrGrJrrwherer!rMr) rArBr+r,rvr)r*rr_Zdiffr4Zpoolsdr5r5r6r s2v       22cCst|sd|d}t|t|ts2d}t||dk r^|dk r^|td|d|}n$|dk rzd|t|}ntdt|d|||d S) aCompute effect size from a T-value. Parameters ---------- tval : float T-value nx, ny : int, optional Group sample sizes. N : int, optional Total sample size (will not be used if nx and ny are specified) eftype : string, optional Desired output effect size. Returns ------- ef : float Effect size See Also -------- compute_effsize : Calculate effect size between two set of observations. convert_effsize : Conversion between effect sizes. Notes ----- If both nx and ny are specified, the formula to convert from *t* to *d* is: .. math:: d = t * \sqrt{\frac{1}{n_x} + \frac{1}{n_y}} If only N (total sample size) is specified, the formula is: .. math:: d = \frac{2t}{\sqrt{N}} Examples -------- 1. Compute effect size from a T-value when both sample sizes are known. >>> from pingouin import compute_effsize_from_t >>> tval, nx, ny = 2.90, 35, 25 >>> d = compute_effsize_from_t(tval, nx=nx, ny=ny, eftype='cohen') >>> print(d) 0.7593982580212534 2. Compute effect size when only total sample size is known (nx+ny) >>> tval, N = 2.90, 60 >>> d = compute_effsize_from_t(tval, N=N, eftype='cohen') >>> print(d) 0.7487767802667672 rcrdzT-value must be floatNrr z*You must specify either nx + ny, or just Nr rx)rrXrrr r!r)Ztvalr)r*Nr,rvrr5r5r6r -s3  )NNNFr r r r) NNr7Fr r8r NF)NN)Fr )NNNr )roZnumpyr rrZpingouin.utilsrr__all__rrrr r r5r5r5r6sH   + / $ "