Performing a Bayesian T-Test in Excel
Objectives
In this post we will compare two samples for differences in their means. The catch is that we want to implement this in a spreadsheet programme (MS Excel), so that it can be easily used by people with no background in R/Python. Microsoft have included many functions allowing for standard frequentist analyses to be performed in Excel, however we will use these functions in order to perform bayesian inference.
We will be guided by the following principles when deciding how to perform the analyses:
- Only use base excel (No VBA). The reason for this is that macro enabled worksheets can be blocked by antivirus software and email programs, which makes it slightly more of a hassle to share.
- Minimise the amount of customisation that needs to be done for each analysis.
- Minimise the amount of data that needs to be stored in the worksheet.
Problem Setup
We have two groups, Target and Control, of size \(n_t\) and \(n_c\) respectively. For these groups we have each of their sample means (\(\bar{x}_t\) and \(\bar{x}_c\)) and sample variances (\(s^2_t\) and \(s^2_c\)).
We will define \(\gamma\) (the standard effect size) as the difference between the true means of the Target and Control group, divided by the standard deviation of the Target group. We parameterise it in this manner so that we do not have to adjust the analysis based on the magnitude of the values.
All the results below will be presented in terms of \(\gamma\). If we wish to return to the original scale we will need to multiply by the standard deviation of the target group, which can be approximated by \(s_t\). Doing so does slightly underestimate the uncertainty (as we are not considering the distribution of the standard deviation around \(s_t\)), however this effect is typically small.
Normality
Assuming that the underlying distribution of values has both a mean and variance, the Central Limit Theorem allows us to say that the sample mean will converge to a Normal (Gaussian) distribution. Basing our calculations on this immediately grants us several advantages:
- The distribution does not have to be changed for different analyses
- The Normal distribution is part of the Exponential Family. One of the properties of this family is that all of the information contained in the data is contained in the sufficient statistics. This means that we do not need to import all of the data into our spreadsheet, we can instead simply input the sufficient statistics. For the normal distribution these are the sample mean and sample variance
- A significant amount of theory around the sampling distributions of the Normal has already been done, and these analyses are common enough that important functions are included in Excel
Note that this does mean that the analysis could be inaccurate for small sample sizes.
Prior Distribution
It is necessary to define a prior distribution for \(\gamma\). For the credibility intervals and parameter estimation we can use an improper uniform prior (ie, all real numbers are equally likely). However we cannot do this for the hypothesis tests, as this will result in our alternative hypothesis having zero likelihood. I have chosen to use the Uniform distribution with a maximum of \(m\) and minimum of \(-m\) as this allows the math to work out simply, meaning that we will not need to use Monte Carlo methods to simulate the posterior. We want to choose a reasonably wide distribution so that we do not ignore large effect sizes. A value of around 2 should be wide enough for most use cases, but this can be adjusted as needed. Note that the wider the prior distribution is the more the hypothesis tests will favour the null hypothesis.
If we wish to use a one-sided prior (if we believe that the effect cannot be negative), then we simply set the minimum to zero.
Conditional Distributions
Conditional on \(\gamma\), we have that \(\bar{x}_t - \bar{x}_c\) will have an approximately Normal distribution (from the Central Limit Theorem). Since we do not know the variances we can convert this into Student’s t-Distribution. We will not assume that the variances of the two groups are equal but this does mean that we need to make some approximations. We will divide the difference in sample means by: \(s_\Delta = \sqrt{\frac{s_t^2}{n_t} + \frac{s_c^2}{n_c}}\). The degrees of freedom are approximately: \(\nu = \frac{s_\Delta^4}{\frac{(s_t^2/n_t)^2}{n_t - 1} +\frac{(s_c^2/n_c)^2}{n_c - 1}}\). This statistic (\(t = \frac{\bar{x}_t - \bar{x}_c}{s_\Delta}\)) has a non-standard t-distribution with a mean of approximately \(\gamma\frac{s_t}{s_\Delta}\).
Posterior Distributions
Since the t-distribution is symmetric, if we assume an improper uniform prior for \(\gamma\), then \(\gamma\frac{s_t}{s_\Delta} - t\) has an approximate t-distribution with \(\nu\) degrees of freedom (conditional on the observed data). If we wish to have a one-sided prior the distribution is proportional to the t-distribution but only where \(\gamma\) is positive, otherwise it is zero.
Inference
With the requisite setup out of the way we can now move on to performing inference. Here I will present the bayesian analgues of common frequentist inferential techniques. ##Hypothesis Testing Firstly we set up our 2 hypotheses: \[H_0: \gamma = 0 \\ H_1: \gamma \sim \mbox{Unif}(-m,m)\]
Or, in the one sided case: \[H_1: \gamma \sim \mbox{Unif}(0,m)\]
For a Bayesian Hypothesis test we calculate the odds ratio: \(K = \frac{Pr(T=t|H_1)}{Pr(T=t|H_0)}\)
This shows how much the data favour the \(H_1\) compared to \(H_0\). The denominator is simply the density of the t-distibution with \(\nu\) degrees of freedom at the observed value of \(t\).
The numerator is given by \(\int_{-m}^m \frac{1}{2m}f_\nu(t-\gamma\frac{s_t}{s_\Delta}) d\gamma\), where \(f_\nu\) is the pdf of a standard t distribution with \(\nu\) degrees of freedom. This simplifies to \(\frac{s_\Delta}{2ms_t}(F_\nu(t+m\frac{s_t}{s_\Delta})-F_\nu(t-m\frac{s_t}{s_\Delta}))\). In the one sided case it instead simplifies to : \(\frac{s_\Delta}{ms_t}(F_\nu(t)-F_\nu(t-m\frac{s_t}{s_\Delta}))\)
Excel contains functions for both the pdf and cdf of the t-distribution (the T.Dist function, which has a boolean parameter ‘cumulative’).
Credibility Intervals
Recall that with an improper uniform prior \(\gamma\frac{s_t}{s_\Delta} - t\) has an approximate t-distribution with \(\nu\) degrees of freedom. Excel has a function for the inverse cumulative t-distribution (T.Inv) so we can obtain a percentile \(\alpha\) of \(\gamma\) as \((F_\nu^{-1}(\alpha)+t)\frac{s_\Delta}{s_t}\). We can create credibility intervals using this. As the t-distribution is symmetric, a symmetric interval will also be the minimum width (highest probability density) interval.
With the one sided prior, the posterior distribution of \(\gamma\frac{s_t}{s_\Delta} - t\) is proportional to \(t_\nu\), except that the density below \(-t\) is zero. We can determine the percentile \(\alpha\) of \(\gamma\) as follows: \[(F_\nu^{-1}\left(1-(1-\alpha)F_\nu(t)\right)+t)\frac{s_\Delta}{s_t}\] In this case the symmetric interval (which has equal probability above and below) is not the same as the minimum width interval. This interval will be symmetric about the mode (which is \(\frac{ts_\Delta}{s_t}\)), unless this results in the lower bound becoming negative, in which case the lower bound is zero and the upper bound will be the percentile corresponding to the required credibility level. To summarise, the lower bound of a \(1-\alpha\) credibility interval is \(\max\left(F_\nu^{-1}\left(\frac{\alpha}{2}F_\nu(t)\right)+t,0\right)\frac{s_\Delta}{s_t}\). The upper bound is either given by the percentile formula presented earlier (when the lower bound is zero), or \(\left(F_\nu^{-1}\left((1-\frac{\alpha}{2})F_\nu(t)\right)+t\right)\frac{s_\Delta}{s_t}\).
Parameter Estimation
In the two-sided case we can determine the mean of \(\gamma\) as \(\frac{ts_\Delta}{s_t}\).
The one-sided case is more complicated. I will present the mean of \(\gamma\) here and leave the derivation as an exercise for the reader:
\[\mbox{E}[\gamma|t,\nu] = \frac{s_\Delta}{s_t}\left(t+\frac{\sqrt{\nu}\Gamma\left({\frac{\nu+1}{2}}\right)}{\sqrt{\pi}(\nu-1)F_\nu(t)\Gamma\left(\frac{\nu}{2}\right)}\left(1+\frac{t^2}{\nu}\right)^{-\frac{\nu-1}{2}}\right) \] There is no gamma function in excel, however there is a GammaLn function which is the natural logarithm of the gamma function. We can therefore obtain the ratio of two gamma functions by taking the difference of the corresponding gammaln functions and exponentiating the result.
Extensions
The above analyses cover the majority of simple use cases. They can be extended to apply to arbitrary cost functions for parameter estimation.
The larger limitation is the restriction on the prior. The uniform prior does not give us the nice ‘regression’ properties that can be observed in bayesian analyses with informative priors. This can be somewhat replicated by including a probability mass at zero in the prior. Another solution is to use the triangular distribution instead of the uniform, as this is still simple enough to implement in Excel.