rwskit.stats.bayesian.ttest

An implementation of Bayesian t-tests.

Attributes

log
InferenceData	Data generated during Bayesian inference.

Classes

BayesFactor

The pair of Bayes Factors for two models \(M_0\) and \(M_1\).

Functions

savage_dickey_t_test(→ BayesFactor)

Perform a Bayesian t-test on the two arrays.

Module Contents

rwskit.stats.bayesian.ttest.log

rwskit.stats.bayesian.ttest.InferenceData

Data generated during Bayesian inference.

class rwskit.stats.bayesian.ttest.BayesFactor

The pair of Bayes Factors for two models \(M_0\) and \(M_1\). Both directions are provided for convenience even though one can always be derived from the other.

bf_10: float

The Bayes Factor \(BF_{10} = \frac{p(y | M_1)}{p(y | M_0)}\).

inference_data: InferenceData | None = None

The trace of the sampling process, i.e., all the samples generated and used to compute the Bayes Factors.

property bf_01: float

The Bayes Factor \(BF_{01} = \frac{1}{BF_{10}} = \frac{p(y | M_0)}{p(y | M_1)}\).

is_significant(significance_level: float = 10.0) → bool

Test if support for the hypothesis is significant.

Returns True if the Bayes Factor in support of the hypothesis model (i.e., \(BF_{10}\)) is greater than the specified significance_level. See this article for a table of common significance levels.

If bf_10 is inf, this function will always return True.

Parameters:

significance_level (float, default = 10.0) – The Bayes Factor to consider evidence in favor of the hypothesis significant.

Returns:

True if bf_10 is greater than the specified significance_level.

Return type:

bool

plot_trace(combined: bool = True, show: bool = False) → matplotlib.axes.Axes | None

Plot the inference trace using the inference_data.

This method fixes the title formatting of arviz.plot_trace(). Without this fix, the subplot titles on the left graph overlap with the x-axis of the subplot just above it making them illegible.

Parameters:

• combined (bool) – If True, the inference chains will all be plotted using a single series.

• show (bool) – If True render the plot to the current device.

Returns:

The matplotlib axes.

Return type:

Axes

rwskit.stats.bayesian.ttest.savage_dickey_t_test(group_1: numpy.typing.ArrayLike, group_2: numpy.typing.ArrayLike, beta_0_range: tuple[float, float] | None = None, cohen_d_scale: float = 1, sigma_scale: float = 10, n_posterior_samples: int = 1000, n_prior_samples: int = 4000, random_seed: int | None = None, return_trace: bool = False) → BayesFactor

Perform a Bayesian t-test on the two arrays.

The test can be used to estimate the likelihood that the true means of two groups are different. The returned Bayes Factor represents the ratio between the probability of the null-hypothesis \(H_{0}\), that the means are equal, and the alternate hypothesis \(H_{1}\), the means are different.

\(BF_{10} > 1\) indicates support that the means are different. \(BF_{10} > 10\) generally indicates strong support that the means are different (\(H_{1}\) is more than 10 times as likely than \(H_{0}\)). Note, \(BF_{10} = \frac{1}{BF_{01}}\).

The default priors are extremely weakly informative. In practice, this means more data is needed to make strong conclusions about the data. If more information is known about the groups it should be used to refine the priors.

See the description of The Model below for the meaning of the parameters: beta_0, cohen_d, and sigma.

Parameters:

• group_1 (array-like) – The first array.

• group_2 (array_like) – The second array

• beta_0_range (tuple[float, float], optional) – The range of values for the Uniform prior on \(\beta_{1}\). If not provided, defaults to (min(group_1 + group_2), max(group_1 + group_2)).

• cohen_d_scale (float) – The scale to use for the cohen_d prior drawn from a Cauchy distribution.

• sigma_scale (float) – The scale to use for model noise drawn from a half Normal distribution. Using a large value will bias the model towards indicating the groups are the same (because the prior belief is most of the variability is due to random noise).

• n_posterior_samples (int, default = 1000) – The number of samples to draw from the posterior distribution.

• n_prior_samples (int, default = 4000) – The number of samples to draw from the predictive prior distribution.

• random_seed (int, default = None) – The random seed to use when drawing samples.

• return_trace (bool, default = False) – If True return the sampled inference data.

Returns:

The Bayes Factors for \(H_{10}\) and \(H_{01}\). When the groups are very different sometimes the posterior of the fit model does not contain the reference value for comparison (i.e., the value 0). In this case bf_10=np.inf and bf_01=0. If this behavior is not caused by an error in the sampling process, it provides extreme support that the group means are different (but probably not infinite support).

Return type:

BayesFactor

The Model

The test is implemented based on the method proposed by Rouder et al. (2009) and described by Maarten Speekenbrink.

The model is defined as:

\[\begin{split}Y &\sim \mathcal{N}(\mu, \sigma_{y}) \\ \sigma_{y} &\sim \mathcal{N^{+}}(\tau_{\sigma_{y}}) \\ d_{c} &\sim \text{Cauchy}(0, \tau_{d_{c}}) \\ \beta_{1} &= d_{c} \cdot \sigma_{y} \\ \beta_{0} &\sim \text{Uniform}(\tau_{\beta_{0}[0]}, \tau_{\beta_{0}[1]}) \\ \mu &= \beta_{0} \cdot \beta_{1} \cdot C_{h}(X) \\\end{split}\]

\(\sigma_{y}\) is the random noise in the process, i.e., the expected standard deviation of the values from the mean du to chance. The prior is drawn from a half Normal distribution with scale \(\tau_{\sigma_{y}}\).

\(X\) is an indicator variable for our two groups. It is coded using Helmert coding, \(C_{h}(X)\), which compares a group mean to the average of the other group means. In the case of 2 groups, it will cause \(\beta_{1}\) to converge to \(\mu_{1} - \mu_{2}\) and \(\beta_{0}\) to \(\frac{\mu_{1} + \mu_{2}}{2}\).

Because the range of values of \(\mu_{1} - \mu_{2}\) is highly dependent data generating process, it is difficult to pick a prior that is likely to work well for any 2 groups we might encounter. The proposed method is to sample from Cohen’s d, \(d_{c} = \frac{\mu_1 - \mu_{2}}{\sigma_{y}}\), instead. \(\beta_{1}\) can be easily recovered by multiplying \(d_{c}\) by \(\sigma_{y}\).

\(\beta_{0}\) is drawn from a Normal distribution with scale \(\tau_{\beta_{0}}\).

\(\mu\) is the linear model using the indicator variable defined above. The intercept is \(\beta_{0}\) and the slope is \(\beta_{1}\).

\(Y \sim \mathcal{N}(\mu, \sigma_{y})\) is the full model that incorporates the random noise.

Hypothesis Testing

We can use the model for hypothesis testing by setting up our two hypotheses as follows:

\[\begin{split}H_0: \mu_1 - \mu_2 &= \beta_{1} &= 0 \\ H_1: \mu_1 - \mu_2 &= \beta_{1} &\neq 0\end{split}\]

Because the models for \(H_{0}\) and \(H_{1}\) are properly nested [1], we can use the Savage-Dickey density ratio to estimate the Bayes Factor. So that:

\[\text{BF}_{01} = \frac{p(\beta_{1}=0 | D, Y)}{p(\beta_{1}=0 | Y)}\]

Where \(D\) is our observed data and \(Y\) is our (linear) model. In words this essentially says, the Bayes Factor in support of the null-hypothesis can be found by taking the ratio of the probability that \(\beta_{1}\) equals 0 in the posterior distribution to the probability that \(\beta_{1}\) equals 0 in the (predictive) prior distribution.

[1]

Briefly, two models are properly nested if they have the same parameters, but in one model, \(M_{1}\), they are all free to take on any value while in the other (at least) one of the parameters is fixed to a constant value.