bayesian.Rmd

---
title: "Bayesian Analysis"
output:
  html_document:
    code_folding: show
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  error = TRUE,
  comment = "",
  class.source = "fold-show")
```

Adapted from `brms` workshop by Paul-Christian Bürkner

# Preamble

## Install Libraries

```{r}
#install.packages("remotes")
#remotes::install_github("DevPsyLab/petersenlab")
```

## Load Libraries

```{r, message = FALSE, warning = FALSE}
library("lme4")
library("rstan")
library("brms")
library("bayestestR")
library("mice")
```

## Simulate Data

```{r, class.source = "fold-hide"}
set.seed(52242)

sampleSize <- 1000

id <- rep(1:100, each = 10)
X <- rnorm(sampleSize)
M <- 0.5*X + rnorm(sampleSize)
Y <- 0.7*M + rnorm(sampleSize)

X[sample(1:length(X), size = 10)] <- NA
M[sample(1:length(M), size = 10)] <- NA
Y[sample(1:length(Y), size = 10)] <- NA

mydata <- data.frame(
  id = id,
  X = X,
  Y = Y,
  M = M)
```

## Load Data

```{r}
data("sleepstudy", package = "lme4")
```

## Prepare Data

```{r}
conditions <- make_conditions(sleepstudy, "Subject")
```

# brms

## Post-Processing Methods

```{r, class.source = "fold-hide"}
methods(class = "brmsfit")
```

# Multilevel Models

## Fit the Models

### Complete Pooling

```{r, results = "hide"}
fit_sleep1 <- brm(
  Reaction ~ 1 + Days,
  data = sleepstudy,
  seed  = 52242)
```

### Random Intercepts

```{r, results = "hide"}
fit_sleep2 <- brm(
  Reaction ~ 1 + Days + (1 | Subject), 
  data = sleepstudy,
  seed  = 52242)
```

### Random Intercepts and Slopes

```{r, results = "hide"}
fit_sleep3 <- brm(
  Reaction ~ 1 + Days + (1 + Days | Subject), 
  data = sleepstudy,
  seed  = 52242
)
```

## Summarize Results

For convergence, `Rhat` values should not be above 1.00.

```{r}
summary(fit_sleep3)
```

## Model Priors

```{r}
prior_summary(fit_sleep3)
```

## Model Parameters

```{r}
variables(fit_sleep3)
```

## Model Coefficients

```{r}
coef(fit_sleep3)
```

## Plots

### Trace Plots

```{r}
plot(fit_sleep3, ask = FALSE)
```

### Visualize Predictions

#### Sample-Level

```{r}
plot(conditional_effects(fit_sleep1), points = TRUE)
```

#### Person-Level

```{r}
# re_formula = NULL ensures that group-level effects are included
ce2 <- conditional_effects(
  fit_sleep3,
  conditions = conditions,
  re_formula = NULL)

plot(ce2, ncol = 6, points = TRUE)
```

### Check Model Fit

#### Posterior Predictive Check

Evaluate how closely the posterior predictions match the observed values.
If they do not match the general pattern of the observed values, a different response distribution may be necessary.

```{r}
pp_check(fit_sleep3)
pp_check(fit_sleep3, type = "dens_overlay")
pp_check(fit_sleep3, "error_scatter_avg")
```

## Fitted Values

```{r}
fitted(fit_sleep3)
```

## Residuals

```{r}
residuals(fit_sleep3)
```

## Compare Models

`elpd` values: higher is better
`looic` values: lower is better

`elpd_diff` values that are greater than ~2 standard errors of the `elpd_diff` values indicate a significantly better model (i.e., if `elpd_diff` value is greater than 2 times the `se_diff` value).

```{r}
loo(fit_sleep1, fit_sleep2, fit_sleep3)
print(loo(fit_sleep1, fit_sleep2, fit_sleep3), simplify = FALSE)
```

### Compute Model Weights

```{r}
model_weights(fit_sleep1, fit_sleep2, fit_sleep3, weights = "loo")

round(model_weights(fit_sleep1, fit_sleep2, fit_sleep3, weights = "loo"))
```

# Multilevel Mediation {#multilevelMediation}

The syntax below estimates random intercepts (which allows each participant to have a different intercept) to account for nested data within the same participant.

```{r}
bayesianMediationSyntax <-
  bf(M ~ X + (1 |i| id)) +
  bf(Y ~ X + M + (1 |i| id)) +
  set_rescor(FALSE) # don't add a residual correlation between M and Y

bayesianMediationModel <- brm(
  bayesianMediationSyntax,
  data  = mydata,
  seed  = 52242
)

summary(bayesianMediationModel)

hypothesis(
  bayesianMediationModel,
  "b_M_X * b_Y_M = 0", # indirect effect = a path * b path
  class = NULL,
  seed  =  52242
)

mediation(bayesianMediationModel)
```

# Setting Priors

## Parameters for Which to Set Priors

```{r}
get_prior(
  Reaction ~ 1 + Days + (1 + Days | Subject), 
  data = sleepstudy)
```

## Define Priors

```{r}
bprior <- c(
  set_prior("normal(5, 5)", coef = "Days"),
  set_prior("cauchy(0, 10)", class = "sd"),
  set_prior("lkj(2)", class = "cor"))

bprior
```

## Fit the Model

Fit the model with these priors, and sample from these priors:

```{r, results = "hide"}
fit_sleep4 <- brm(
  Reaction ~ 1 + Days + (1 + Days | Subject), 
  data = sleepstudy,
  prior = bprior, 
  sample_prior = TRUE,
  seed  = 52242
)
```

# Evaluate a Hypothesis

```{r}
# Evid.Ratio is the ratio of P(Days > 7) / P(Days <= 7)
(hyp1 <- hypothesis(fit_sleep4, "Days < 7"))
plot(hyp1)
```

```{r}
# Evid.Ratio is the Bayes Factor of the posterior
# vs the prior that Days = 10 is TRUE (Savage-Dickey Ratio)
(hyp2 <- hypothesis(fit_sleep4, "Days = 10"))
plot(hyp2)
```

# Parallel Processing

## Between-Chain Parallelization

```{r, results = "hide"}
fit_sleep4 <- brm(
  Reaction ~ 1 + Days + (1 + Days | Subject), 
  data = sleepstudy,
  prior = bprior, 
  sample_prior = TRUE,
  cores = 4,
  seed  = 52242
)
```

## Within-Chain Parallelization

https://paul-buerkner.github.io/brms/articles/brms_threading.html (archived at https://perma.cc/NCG3-KV4G)

# Missing Data Handling

## Multiply Imputed Datasets from `mice`

```{r}
?brm_multiple
```

```{r}
imp <- mice::mice(
  mydata,
  m = 5,
  print = FALSE)

fit_imp <- brm_multiple(
  bayesianMediationSyntax,
  data = imp,
  chains = 2)
```

## Imputation on the Fly During Model Fitting

https://paul-buerkner.github.io/brms/articles/brms_missings.html (archived at https://perma.cc/4Y9L-USQR)

```{r}
?mi
```

```{r}
bayesianRegressionImputationSyntax <-
  bf(X | mi() ~ (1 |i| id)) +
  bf(M | mi() ~ mi(X) + (1 |i| id)) +
  bf(Y | mi() ~ mi(X) + mi(M) + (1 |i| id)) +
  set_rescor(FALSE) # don't add a residual correlation between X, M, and Y

bayesianRegressionModel <- brm(
  bayesianRegressionImputationSyntax,
  data = mydata,
  seed = 52242
)

summary(bayesianRegressionModel)

hypothesis(
  bayesianRegressionModel,
  "bsp_M_miX * bsp_Y_miM = 0", # indirect effect = a path * b path
  class = NULL,
  seed = 52242
)
```

# Session Info

```{r, class.source = "fold-hide"}
sessionInfo()
```