In brief

Reporting our results a is key part of the research cycle. We need to aim to report findings fairly, and set them in a meaningful context. Our goal is not to make results as impressive as possible: you should aim for a balanced summary of the data which another researcher could and would replicate from your published report.

An important task is to quantify how much evidence the model provides. In recent years, null-hypothesis significance tests (NHST) have fallen out of favour among (although there is a strong inertia in research practice). Courses at Plymouth explain NHST to enable you to read the existing literature. But we teach Bayes Factors and other Bayesian techniques which enable you to quantify the strength of the evidence in favour, or against, a hypothesis of interests.

Before you start

How good?

We said we could ask two questions about how good our regression models are:

1. How much variance (variability) in the outcome do the predictors explain?
2. How much more probable is the model which includes our predictors, than a model without them?

To answer question 1 we can check what proportion of variance is explained byt the model using \(R^2\) (see the previous session).

The second question asks how probable different models are. By ‘probable’ we mean:

How likely is it that our model is ‘correct’, given the data we have observed.

For various reasons, it’s very hard to actually calculate that probability for an individual model. We also know that ALL models are a simplification of reality, so their absolute probability is not always that important.

However, it’s much easier to compare two different models, and ask the slightly different question:

How much more likely is my model than some alternative model?

If we choose suitable alternative models for comparison this can be a powerful tool.

For example, we can compare models with and without certain predictors included, and this gives us a sense of how much evidence we have that they help predict the outcome.

This kind of comparison is expressed using a Bayes Factor.

Bayes Factors

Was worth running our model at all? Maybe our predictors are no good and don’t really help us explain the outcome? A Bayes Factor quantifies how likely that is.

A Bayes Factor can compare the probabilities of different hypotheses if those hypotheses are represented by statistical models, for example:

‘No effects’ hypothesis, \(H_0\)

The predictors are NOT related to the outcome (except by chance), v.s.

Experimental hypothesis, \(H_1\)

The predictors are related to the outcome

library(BayesFactor)

lmBF(grade ~ work_hours * female, data = studyhabits)

Bayes factor analysis
--------------
[1] work_hours * female : 1.084058e+13 ±0%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

lmBF(grade ~ work_hours, data = studyhabits)

Bayes factor analysis
--------------
[1] work_hours : 726159159273 ±0.01%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

Testing interactions

The BF we just calculated only tested whether work hours was a useful predictor. It didn’t include gender as a predictor.

We can also add female using lmBF:

lmBF(grade ~ work_hours * female, data = studyhabits)

Bayes factor analysis
--------------
[1] work_hours * female : 1.084058e+13 ±0%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

Explanation of the code: We added-back female and the interaction of work_hours and female (see the previous session for a reminder about adding multiple predictors)

Explanation of the output:

• When we added the interaction the Bayes Factor got bigger: 1.0840583^{13}, vs 1.6874571^{11} before.

• This suggests that the model with female and the interaction is more likely than the model without

We want to compare these two hypotheses.

• \(H_0\): The model including work_hours, and a difference for female, is the best
• \(H_1\): The model including work_hours and difference for female and their interaction is the best

A neat think about Bayes Factors is that you are not restricted to comparing against the intercept only model. You can compare any 2 models (which use the same outcome and dataset) by dividing their Bayes Factors.

We do this by simply dividing the Bayes Factors for each model

• BF₁₀ = \(\dfrac{BF(H_1)}{BF(H_0)}\) (> 3 favours \(H_1\))

• BF₀₁ = \(\dfrac{BF(H_0)}{BF(H_1)}\) (> 3 favours \(H_0\))

So in this case:

• BF₁₀ = \(\dfrac{1.084058e+13}{1.687457e+11} = 64.24\)

And

• BF₀₁ = \(\dfrac{1.687457e+11}{1.084058e+13} = 0.02\)

Explanation of the result: The BF₁₀ is 64.24, which means we have very strong evidence in favour of the interaction.

lmBF(grade ~ work_consistently + female, data = studyhabits)

Bayes factor analysis
--------------
[1] work_consistently + female : 2838750 ±0%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

lmBF(grade ~ work_consistently * female, data = studyhabits)

Bayes factor analysis
--------------
[1] work_consistently * female : 554911.9 ±0.01%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

Reversing a Bayes Factor

Sometimes, you’ll see Bayes Factor written as BF₁₀, which means the same thing as just plain ‘BF’.

Occasionally you will see BF₀₁, which is the same idea but flipped. BF₀₁ is the inverse of BF₁₀: it provides evidence in favour of \(H_0\).

BF₀₁ > 3 means substantial evidence for \(H_0\) (where \(H_0\) is often the ‘null’ hypothesis of no-difference, or no-relationship).

Flipping a BF from a BF₀₁ to a BF₁₀ (or the reverse) is easy. Just calculate:

\(\dfrac{1}{BF}\)

So, if your BF₀₁ = 0.2, then your BF₁₀ would just be \(1/0.2 = 5\)

Among us

Imagine we had data where students who reported how much Among Us they played each day. We could plot the effect on their achievement like this:

As we have seen before, the lines on this plot can be represented by the coefficients in a linear model.

Our model would include a continuous variable to measure hours playing (gameplay), a categorical variable (gender), and the interaction of these two variables:

Call:
lm(formula = score ~ gameplay * gender, data = amongus)

Coefficients:
          (Intercept)               gameplay           genderFemale  gameplay:genderFemale  
               25.403                 -1.703                  2.260                 -2.617  

Explanation of the output:

• In previous examples we included a variable for gender which was called female and coded 0/1 or FALSE/TRUE for men/women respectively.
• In this example we have included a variable called gender which has the text "Male" or "Female" instead.
• The coefficients returned are the same as they would be with a 0/1 coding, but the names of the coefficients are different
• For example genderFemale now relates to the differences in intercepts between men and women, and gameplay:genderFemale now relates to the difference in slopes.

h0 <- lmBF( score ~ gameplay + gender, data=amongus)
h1 <- lmBF( score ~ gameplay * gender, data=amongus)

# calculate the BF10
h1/h0

Bayes factor analysis
--------------
[1] gameplay * gender : 1.917127 ±1.28%

Against denominator:
  score ~ gameplay + gender 
---
Bayes factor type: BFlinearModel, JZS

# calculate the BF01
h0/h1

Bayes factor analysis
--------------
[1] gameplay + gender : 0.5216139 ±1.28%

Against denominator:
  score ~ gameplay * gender 
---
Bayes factor type: BFlinearModel, JZS