Step 1: Specify priors for hypotheses

We have an experiment where 20 rats were randomised to one of four doses of the antidepressant fluoxetine, given in the drinking water. The time that the rats spent immobile in the Forced Swim Test (FST) was recorded. The FST is a standard behavioural test of “depression” in rodents.

First, we specify the prior probability of three hypotheses: that fluoxetine increases (H>), decreases (H<), or has no effect (H0) on immobility time. These three hypotheses are exhaustive (include all possible outcomes) and mutually exclusive (only one can be true). Although fluoxetine is known to decrease immobility time, we will specify an equality probability of 1/3 for each hypothesis to illustrate the approach.

Step 2: Specify a prior for the effect size

Next, we need to specify a prior for the effect size (we define the effect size in the Step 3). For now we will use the default prior, which is calculated from the data. It is a normal prior, centred at zero, with the width calculated such that the 99% confidence interval (CI) of the prior matches the 99% CI of the data distribution. The results can be sensitive to the choice of prior, especially if the width is large relative to the data. The default prior is suitable for most situations where you have no prior information to include, and the results are insensitive to small changes near the default value (see below).

Step 3: Calculate effect size and standard error

The data for this example are in the labstats package (available on CRAN) and plotted in Figure 2. To calculate the effect size (and the prior in the previous step) we need to define the analysis. Here, dose is treated as a continuous variable (see reference [5]) and so a linear regression quantifies the relationship between fluoxetine and immobility time. No effect corresponds to a flat line (slope = 0) in Figure 2.

library(labstats)
par(las=1)
plot(time.immob ~ dose, data=fluoxetine, col="royalblue",
    ylab="Time immobile (s)", ylim=c(0, 250), xlim=c(0, 250),
    xlab="Dose of fluoxetine (mg/L in drinking water)")
abline(lm(time.immob ~ dose, data=fluoxetine), col="royalblue") # add reg. line

Effect of fluoxetine (Prozac) on rats in the Forced Swim Test. Data are from Lazic [4]

The code below calculates the effect size and standard error. These are the required inputs to calculate the posterior probabilities and are returned from lm(), glm(), and related functions (e.g. from ANOVAs, t-tests, or regressions with Gaussian, Poisson, or binomial outcomes).

summary(lm(time.immob ~ dose, data=fluoxetine))$coef

##             Estimate Std. Error  t value    Pr(>|t|)
## (Intercept)  170.440 14.8368131 11.48764 1.01541e-09
## dose          -0.252  0.0991326 -2.54205 2.04365e-02

From the above output we see that the estimated slope is -0.252, with a standard error of 0.099, and a p-value of 0.020. We now have all the information to calculate the probability that fluoxetine increases, decreases, or has no effect on immobility time.¹

Step 4: Calculate Bayes factors

Bayes factors are used as an intermediate step in calculating the posterior probabilities of each hypothesis. The functions in the BayesCombo package calculate these automatically, but we take a brief digression from the fluoxetine example to illustrate how they are calculated to provide some insight into the method. Figure 3 shows a standard Bayesian updating of a prior distribution to a posterior distribution based on the data (likelihood). The prior in the top panel is normal with a mean of zero and standard deviation of 1.29. The middle panel shows the data (likelihood) distribution, which has a mean of 0.75 and standard deviation of 1. The bottom panel shows the posterior, which is shifted to positive values relative to the prior and reflects the influence of the data. The prior and likelihood correspond to steps 2 and 3 in Figure 1. The BFs for the three hypotheses are then calculated as the ratio of posterior to prior areas or heights of the distributions:

\[ BF_{H<0} = \dfrac{\mathrm{area~of~d}}{\mathrm{area~of~a}} \]

\[ BF_{H=0} = \dfrac{\mathrm{height~of~point~e}}{\mathrm{height~of~point~b}} \]

\[ BF_{H>0} = \dfrac{\mathrm{area~of~f}}{\mathrm{area~of~c}} \]

For example, 50% of the prior distribution is above 0 (region c), as is 72% of the posterior (region f). The interpretation is that the data have increased the plausibility of hypothesis H> from 50% to 72%. The ratio of these values is the Bayes factor and is equal to 0.72/0.5 = 1.4. Thus, we can say that H> is 1.4 times more likely. More generally, a BF > 1 means that the data support a hypothesis, whereas a BF < 1 means that data do not support a hypothesis.

Prior, Likelihood (data), and posterior distributions for an experiment. Bayes factors for the three hypotheses are calculated as the ratio of posterior to prior areas or heights.

Step 5: Use the BFs to update the prior probability of each hypothesis

The final step is to update the prior probability for each hypothesis with the BFs to get the posterior probability for each hypothesis. The equation below shows the calculation for the hypothesis that the effect is greater than zero (H>), and other probabilities are calculated in the same way, just substituting other BFs in the numerator. \(Pr()\) are the prior probabilities for each hypothesis, and they cancel out from the equation when they are all equal (e.g. if they are all 1/3).

\[ P(H>) = \dfrac{Pr(H>)BF_{H>}}{Pr(H<)BF_{H<} + Pr(H0)BF_{H0} + Pr(H>)BF_{H>}}\]

All of the above steps can be conveniently calculated using the pph() function (Posterior Probability of a Hypothesis). Returning to the fluoxetine example, we can calculate the probability that the slope is negative, positive, or zero. Below, we specify the slope (beta = -0.252) and its standard error (se.beta = 0.099) that we obtained previously from the output of the lm() function. The default settings are used for all other options. The output below shows that the probability of a negative slope (H<) is 0.9120, a positive slope (H>) is 0.0106, and a zero slope (H0, corresponding to no effect) is 0.0774. Unlike a p-value, which is the probability of the data given a hypothesis, these probabilities have a direct and intuitive interpretation as the probability of a hypothesis given the data.

x <- pph(beta = -0.252, se.beta = 0.099)
summary(x)

##     H<     H0     H> 
## 0.9120 0.0774 0.0106

Plotting the output of the pph() function returns the likelihood, prior, and posterior distributions (Fig. 4). The mean of the likelihood (data distribution; dotted line) is centred on -0.252 and the standard error (0.099) determines the width of this distribution. The prior (thick grey line) was automatically calculated and we can see how it encompasses the likelihood distribution and is centred at zero. The posterior distribution (thin black line) represents the combination of prior and likelihood. The BFs and posterior probabilities of the three hypotheses are then calculated from these prior and posterior distributions (analogous to the distributions in Fig. 3).

par(las=1)
plot(x, leg.loc = "topright", xlab="Effect size (slope of regression line)")

Prior, likelihood (data), and posterior distributions for the fluoxetine data.

The above default prior may be more informative than desired (note how the posterior is pulled towards the prior in Fig. 4) and the easiest way to decrease the influence of the prior is to make it wider by specifying a multiplier. The code below doubles the previous standard error (se.mult = 2), and the posterior is now much closer to the data distribution (Fig. 5). However, the posterior hypothesis probabilities are similar to the previous analysis; rounded to two decimal places, the probability that the slope is negative is still 0.91.

x2 <- pph(beta = -0.252, se.beta = 0.099, se.mult = 2)
summary(x2)

##     H<     H0     H> 
## 0.9051 0.0887 0.0062

par(las=1)
plot(x2, leg.loc = "topright", xlab="Effect size (slope of regression line)")

Prior, likelihood (data), and posterior distributions for the fluoxetine data. Same as Fig. 4 but with a wider prior.

A final example illustrates other options. The prior for the slope is directly specified as having a mean of 0 (beta0 = 0) and standard error of 1.2 (se0 = 1.2). In the previous analyses H0 was defined as exactly equal to 0, but here we define H0 as a range of values close to zero using H0 = c(-0.05, 0.05). Finally, the priors on the hypotheses are also given as an argument to H.priors. The values indicate that the prior probability of the slope being negative, zero, and positive are 0.495, 0.495, 0.01, respectively. The interpretation is that we expect that fluoxetine either decreases immobility time or has no effect, but it is unlikely to increase immobility time.

x3 <- pph(beta = -0.252, se.beta = 0.099, beta0=0, se0=1.2,
          H0 = c(-0.05, 0.05), H.priors=c(0.495, 0.495, 0.01))
summary(x3)

##    H<    H0    H> 
## 0.984 0.016 0.000

par(las=1)
plot(x3, leg.loc = "topright", xlab="Effect size (slope of regression line)")

Prior, likelihood (data), and posterior distributions for the fluoxetine data. Same as Fig. 4 but with alternative options specified.

With these options, the probability that fluoxetine decreases immobility time is now 0.98, and the distributions are shown in Figure 6.