EDUC 847 Winter 24

class: center, middle, inverse, title-slide

.title[
# EDUC 847 Winter 24
]
.subtitle[
## Week 8 - BayesFactors
]
.author[
### Eric Brewe <br> Professor of Physics at Drexel University <br>
]
.date[
### 11 March 2024, last update: 2024-03-12
]

---

# Let's start with a complete dataset

Start by getting the data into R

```r
#This loads the csv and saves it as a dataframe titled sci_df

sci_df <- read_csv(here("data",
                        "sci_scores.csv"))
```

---
# Let's remember prior work

### We predicted science scores from motivation scores.

```r
*sci_mod_motiv <- lm(sci ~ motiv, data = sci_df)

tidy(summary(sci_mod_motiv))
```

```
## # A tibble: 2 × 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)   405.      179.        2.26  0.0242
## 2 motiv           5.46      2.42      2.25  0.0247
```

---
# Let's decide if it is a good model.

```r
tidy(glance(sci_mod_motiv)) %>% select(column, mean)
```

```
## # A tibble: 12 × 2
##    column            mean
##    <chr>            <dbl>
##  1 r.squared      1.01e-2
##  2 adj.r.squared  8.10e-3
##  3 sigma          6.64e+2
##  4 statistic      5.07e+0
##  5 p.value        2.47e-2
##  6 df             1   e+0
##  7 logLik        -3.96e+3
##  8 AIC            7.92e+3
##  9 BIC            7.93e+3
## 10 deviance       2.20e+8
## 11 df.residual    4.98e+2
## 12 nobs           5   e+2
```

---
# Let's compare to the null model

## What would the null model be?

```r
*sci_mod_null <- lm(sci ~ 1, data = sci_df)

tidy(summary(sci_mod_null))
```

```
## # A tibble: 1 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)     804.      29.8      27.0 1.87e-99
```

---
# Let's get the AIC for the null

```r
tidy(glance(sci_mod_null)) %>% select(column, mean)
```

```
## # A tibble: 12 × 2
##    column              mean
##    <chr>              <dbl>
##  1 r.squared             0 
##  2 adj.r.squared         0 
##  3 sigma               667.
##  4 statistic           NaN 
##  5 p.value             NaN 
##  6 df                  NaN 
##  7 logLik            -3960.
##  8 AIC                7924.
##  9 BIC                7933.
## 10 deviance      221799764.
## 11 df.residual         499 
## 12 nobs                500
```
AIC went down.

---
# Let's use ANOVA

```r
anova(sci_mod_null,sci_mod_motiv)
```

```
## Analysis of Variance Table
## 
## Model 1: sci ~ 1
## Model 2: sci ~ motiv
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
## 1    499 221799764                              
## 2    498 219562843  1   2236921 5.0737 0.02473 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

This suggests we can reject the null hypothesis.

---
# Let's try another variable

## What if we wanted to investagte group differences

```r
*sci_mod_gr <- lm(sci ~ gr, data = sci_df)

tidy(summary(sci_mod_gr))
```

```
## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    713.       69.9     10.2  2.58e-22
## 2 gr              30.4      21.2      1.43 1.52e- 1
```

---
# Let's decide if it is a good model.

```r
tidy(glance(sci_mod_gr)) %>% select(column, mean)
```

```
## # A tibble: 12 × 2
##    column            mean
##    <chr>            <dbl>
##  1 r.squared      4.12e-3
##  2 adj.r.squared  2.12e-3
##  3 sigma          6.66e+2
##  4 statistic      2.06e+0
##  5 p.value        1.52e-1
##  6 df             1   e+0
##  7 logLik        -3.96e+3
##  8 AIC            7.92e+3
##  9 BIC            7.94e+3
## 10 deviance       2.21e+8
## 11 df.residual    4.98e+2
## 12 nobs           5   e+2
```

---
# Let's use ANOVA

```r
anova(sci_mod_null,sci_mod_gr)
```

```
## Analysis of Variance Table
## 
## Model 1: sci ~ 1
## Model 2: sci ~ gr
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1    499 221799764                           
## 2    498 220886978  1    912786 2.0579  0.152
```

Doesn't look like we can reject the null. So...?

---
# Let's check out Bayes Factors

## Bayes factors draw on Bayes' Theorem

The book, Learning Statistics with R has a great description of Bayes' rule, but one of the keys is that in Bayesian thinking is updating.

### https://learningstatisticswithr.com/book/bayes.html

Here is how it works:
1. Start with some set of beliefs (called the Prior or Prior Belief)
2. Choose a likelihood function (sometimes this is the Conditional Probability)
3. Multiply the Prior with the Likelihood function to update the beliefs. 
4. Divide by the probability of the data to give the posterior probability (typically called the Posterior)

The equation for the Posterior probability is:

`$P(h|d) = \frac{P(d|h)*P(h)}{P(d)}$`

---
# Let's use Bayes Theory in Hypothesis Tests

## For the null hypothesis we can write the posterior probability as:

`$P(h_{0}|d) = \frac{P(d|h_{0})*P(h_{0})}{P(d)}$`

## For the alternative hypothesis we can write the posterior probability as:

`$P(h_{1}|d) = \frac{P(d|h_{1})*P(h_{1})}{P(d)}$`

---
# Let's take an odds ratio

## In principle, we can compare these models by taking the odds ratio of the probablity of Model 1 to the Null Model.

`$\frac{Model_{1}}{Model_{0}}$`

## What this tells us is how much more probable Model 1 is than Model 0.

---
# Let's do some math

## When we actually do this division...

`$\frac{Model_{1}}{Model_{0}} = \frac{P(h_{1}|d)}{P(d|h_{0})} = \frac{\frac{P(d|h_{1})*P(h_{1})}{P(d)}}{\frac{P(d|h_{0})*P(h_{0})}{P(d)}}$`

Which is great, because a number of things divide out, and we end up with

`$Posterior Odds = Bayes Factor * Prior Odds$`

`$\frac{P(h_{1}|d)}{P(h_{0}|d)} = \frac{P(d|h_{1})}{P(d|h_{0})} * \frac{P(h_{1})}{P(h_{0})}$`

---
# Let's reflect.

## Here are some really great features of this approach.

- What is awesome about this Bayes Factor approach is that it works for any pair of models.  This means we can do this with t-tests or regression models.  
- The results are sort of easy to interpret, how much more probable is one model than the other. 
- The ratio can be inverted too - so if we want to know how much more probable the null model is than the alternative, just do the division the other way - which is really powerful when you are doing equity work.

---
# Let's give it a shot.

## We can start with the model for motivation

First we have to estimate the model.

```r
regressionBF(
  formula = sci ~ motiv,
  data = sci_df
)
```

```
## Bayes factor analysis
## --------------
## [1] motiv : 1.158637 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS
```

Hmm, our Bayes Factor is only 1.15, so not a ton better than the null model.

---
# Let's try a model with groups

## Remember we know we couldn't reject the null.

```r
regressionBF(
  formula = sci ~ gr,
  data = sci_df
)
```

```
## Bayes factor analysis
## --------------
## [1] gr : 0.2697961 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS
```
So, this one is less than 1. To get the confidence in the null model over the model with group, we can just take the inverse, so `$1/0.2697961 = 3.706503$`

---
# Let's see if there is a better model

```r
models <- regressionBF(
  formula = sci ~ .,
  data = sci_df
)
```

---
# Let's see if there is a better model

```r
topModels<-head(models,5)

topModels
```

```
## Bayes factor analysis
## --------------
## [1] motiv + pre_sc + math_sc           : 1.057403e+379 ±0%
## [2] id + motiv + pre_sc + math_sc      : 2.164667e+377 ±0%
## [3] gr + motiv + pre_sc + math_sc      : 1.85498e+377  ±0%
## [4] motiv + pre_sc + math_sc + pass    : 1.803902e+377 ±0%
## [5] id + gr + motiv + pre_sc + math_sc : 4.266269e+375 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS
```

A Bayes Factor of 1.057403e+379!  Is that good?

---
# Let's see the heuristics

## There isn't a clear rule about what counts as significant evidence

|Bayes factor	| Interpretation|
|--------|------|
|1 - 3	|Negligible evidence|
|3 - 20	|Positive evidence|
|20 - 150	|Strong evidence|
|$>$150	|Very strong evidence|

So, yeah a Bayes factor of 1.057403e+379 is really strong evidence.

---
# Let's compare models

## One of the great features of Bayes Factors is that you can use the same framework to compare models.

So is the top model (motiv + pre_sc + math_sc) much better than the third best (gr + motiv + pre_sc + math_sc) model?

```r
topModels[1]/topModels[3]
```

```
## Bayes factor analysis
## --------------
## [1] motiv + pre_sc + math_sc : 57.00344 ±0%
## 
## Against denominator:
##   sci ~ gr + motiv + pre_sc + math_sc 
## ---
## Bayes factor type: BFlinearModel, JZS
```
That tells us how confident to be in Model 1 vs. Model 3.

---
# Let't consider dropping variables

### To identify variables to drop, we can use the whichModels function

```r
regressionBF( 
 formula = sci ~ motiv + pre_sc + math_sc + gr,
 data = sci_df,
*whichModels = "top"
)
```

```
## Bayes factor top-down analysis
## --------------
## When effect is omitted from motiv + pre_sc + math_sc + gr , BF is...
## [1] Omit gr      : 57.00344      ±0%
## [2] Omit math_sc : 7.07287e-60   ±0%
## [3] Omit pre_sc  : 5.822091e-374 ±0.01%
## [4] Omit motiv   : 1.645125e-08  ±0%
## 
## Against denominator:
##   sci ~ motiv + pre_sc + math_sc + gr 
## ---
## Bayes factor type: BFlinearModel, JZS
```

---
# Let't decide what to drop

### The Bayes factor for dropping gr is ~ 57, so yeah get rid of it.

### The Bayes factor for dropping the others are really small. So motiv, pre_sc, and math_sc should all be in the model

---
# Let's do one more

```r
BestModels <- generalTestBF(
  formula = sci ~ motiv + math_sc*pre_sc,
  data = sci_df)
BestModels
```

```
## Bayes factor analysis
## --------------
## [1] motiv                                     : 1.158637       ±0%
## [2] math_sc                                   : 3703.986       ±0%
## [3] motiv + math_sc                           : 18803.81       ±0%
## [4] pre_sc                                    : 3.557574e+318  ±0.01%
## [5] motiv + pre_sc                            : 4.611555e+319  ±0%
## [6] math_sc + pre_sc                          : 2.007702e+371  ±0%
## [7] motiv + math_sc + pre_sc                  : 1.057403e+379  ±0%
## [8] math_sc + pre_sc + math_sc:pre_sc         : 5.133498e+620  ±0%
## [9] motiv + math_sc + pre_sc + math_sc:pre_sc : 9.857432e+1227 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS
```

---
# Let's do one more

```r
BestModels[8]
```

```
## Bayes factor analysis
## --------------
## [1] math_sc + pre_sc + math_sc:pre_sc : 5.133498e+620 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS
```

## Hey! That is the way I specified the data to be simulated, so well done!
---
# Let's review

We spent today looking at how to use Bayes factords

We can...
 - calculate Bayes factors
 - use Bayesian regression to compare models
 - interpret Bayes factors
 - use Bayes factors to determine what variable to include (and those to drop)
 
---
# Let's reflect

--
.font150[.center[YOU DID IT!]]

--
.font150[.center[Congratulations, this wasn't easy!]]

--
.font150[.center[I appreciate your effort, patience, and committment.]]

---
# Let's reflect

.font190[.center[Thank you!]]