rwskit.stats.bayesian.ttest
===========================
.. py:module:: rwskit.stats.bayesian.ttest
.. autoapi-nested-parse::
An implementation of Bayesian t-tests.
Attributes
----------
.. autoapisummary::
rwskit.stats.bayesian.ttest.log
rwskit.stats.bayesian.ttest.InferenceData
Classes
-------
.. autoapisummary::
rwskit.stats.bayesian.ttest.BayesFactor
Functions
---------
.. autoapisummary::
rwskit.stats.bayesian.ttest.savage_dickey_t_test
Module Contents
---------------
.. py:data:: log
.. py:data:: InferenceData
Data generated during Bayesian inference.
.. py:class:: BayesFactor
The pair of Bayes Factors for two models :math:`M_0` and :math:`M_1`.
Both directions are provided for convenience even though one can always be
derived from the other.
.. py:attribute:: bf_10
:type: float
The Bayes Factor :math:`BF_{10} = \frac{p(y | M_1)}{p(y | M_0)}`.
.. py:attribute:: inference_data
:type: Optional[InferenceData]
:value: None
The trace of the sampling process, i.e., all the samples generated and
used to compute the Bayes Factors.
.. py:property:: bf_01
:type: float
The Bayes Factor :math:`BF_{01} = \frac{1}{BF_{10}} = \frac{p(y | M_0)}{p(y | M_1)}`.
.. py:method:: 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., :math:`BF_{10}`) is greater than the specified
``significance_level``. See `this article `__
for a table of common significance levels.
If :data:`bf_10` is ``inf``, this function will always return ``True``.
:param significance_level: The Bayes Factor to consider evidence in favor of the hypothesis
significant.
:type significance_level: float, default = 10.0
:returns: ``True`` if :data:`bf_10` is greater than the specified ``significance_level``.
:rtype: bool
.. py:method:: plot_trace(combined: bool = True, show: bool = False) -> Optional[matplotlib.axes.Axes]
Plot the inference trace using the :data:`inference_data`.
This method fixes the title formatting of :func:`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.
:param combined: If ``True``, the inference chains will all be plotted using a
single series.
:type combined: bool
:param show: If ``True`` render the plot to the current device.
:type show: bool
:returns: The matplotlib axes.
:rtype: Axes
.. py:function:: savage_dickey_t_test(group_1: numpy.typing.ArrayLike, group_2: numpy.typing.ArrayLike, beta_0_range: Optional[tuple[float, float]] = None, cohen_d_scale: float = 1, sigma_scale: float = 10, n_posterior_samples: int = 1000, n_prior_samples: int = 4000, random_seed: Optional[int] = 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 :math:`H_{0}`, that the means are
equal, and the alternate hypothesis :math:`H_{1}`, the means are different.
:math:`BF_{10} > 1` indicates support that the means are different.
:math:`BF_{10} > 10` generally indicates strong support that the means
are different (:math:`H_{1}` is more than 10 times as likely than
:math:`H_{0}`). Note, :math:`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 :ref:`model` below for the meaning of the parameters:
``beta_0``, ``cohen_d``, and ``sigma``.
:param group_1: The first array.
:type group_1: array-like
:param group_2: The second array
:type group_2: array_like
:param beta_0_range: The range of values for the Uniform prior on :math:`\beta_{1}`. If
not provided, defaults to
``(min(group_1 + group_2), max(group_1 + group_2))``.
:type beta_0_range: tuple[float, float], optional
:param cohen_d_scale: The scale to use for the ``cohen_d`` prior drawn from a Cauchy
distribution.
:type cohen_d_scale: float
:param sigma_scale: 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).
:type sigma_scale: float
:param n_posterior_samples: The number of samples to draw from the posterior distribution.
:type n_posterior_samples: int, default = 1000
:param n_prior_samples: The number of samples to draw from the predictive prior
distribution.
:type n_prior_samples: int, default = 4000
:param random_seed: The random seed to use when drawing samples.
:type random_seed: int, default = None
:param return_trace: If ``True`` return the sampled inference data.
:type return_trace: bool, default = False
:returns: The Bayes Factors for :math:`H_{10}` and :math:`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).
:rtype: BayesFactor
.. _model:
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:
.. math::
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) \\
:math:`\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
:math:`\tau_{\sigma_{y}}`.
:math:`X` is an indicator variable for our two groups. It is
coded using Helmert coding, :math:`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 :math:`\beta_{1}` to converge to :math:`\mu_{1} - \mu_{2}`
and :math:`\beta_{0}` to :math:`\frac{\mu_{1} + \mu_{2}}{2}`.
Because the range of values of
:math:`\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 `__,
:math:`d_{c} = \frac{\mu_1 - \mu_{2}}{\sigma_{y}}`, instead.
:math:`\beta_{1}` can be easily recovered by multiplying
:math:`d_{c}` by :math:`\sigma_{y}`.
:math:`\beta_{0}` is drawn from a Normal distribution with scale
:math:`\tau_{\beta_{0}}`.
:math:`\mu` is the linear model using the indicator variable defined above.
The intercept is :math:`\beta_{0}` and the slope is :math:`\beta_{1}`.
:math:`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:
.. math::
H_0: \mu_1 - \mu_2 &= \beta_{1} &= 0 \\
H_1: \mu_1 - \mu_2 &= \beta_{1} &\neq 0
Because the models for :math:`H_{0}` and :math:`H_{1}` are
properly nested [1]_, we can use the
`Savage-Dickey density ratio `__
to estimate the Bayes Factor. So that:
.. math::
\text{BF}_{01} = \frac{p(\beta_{1}=0 | D, Y)}{p(\beta_{1}=0 | Y)}
Where :math:`D` is our observed data and :math:`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
:math:`\beta_{1}` equals 0 in the posterior distribution to the
probability that :math:`\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, :math:`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.