---
title: "Getting started with speedyBBT"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with speedyBBT}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup, warning = FALSE}
library(speedyBBT)
library(coda)
```

`speedyBBT` is a package for Bayesian Bradley-Terry modelling. It provides functions for fitting the Bradley-Terry model using Markov Chain Monte Carlo (MCMC) methods, allowing for inference on player abilities and hyperparameters. It uses latent Pólya-Gamma variables to facilitate efficient sampling.

To fit a Bradley-Terry model using `speedyBBT`, you can first install the package from CRAN or install the development version from GitHub. 

```{r install_pkg}
# Install from CRAN
# install.packages("speedyBBT")
# Install the development version from GitHub
# remotes::install_github("rowlandseymour/speedyBBT")
```

You can learn more about the package functions and how to use them from the documentation which can be accessed using `?function_name` in R. For example, to learn about the `speedyBBTm` function, you can run:

```{r help_function}
?speedyBBTm
```

The package includes some example datasets that you can use to practice fitting models. These datasets are automatically loaded when you install the package. To view the `forcedMarriage` dataset, you can run:
````{r load_data}
data("forcedMarriage")
````

This dataset was collected from 12 judges in Nottinghamshire who were asked to judge which of a pair of 76 wards had a higher risk of forced marriage. The dataset contains the results of the 1,846 pairwise comparisons made by the judges. An adjacency matrix is constructed based on proximity, assuming that geographically close wards are more highly correlated. This serves as a covariance matrix for the multivariate normal prior distribution of the item qualities. In this case, the "items" are wards in Nottinghamshire and their "quality" is the risk of forced marriage in that ward. 

The `speedyBBTm` function can be used to fit a Bradley-Terry model to the data. To start with, you need to construct the prior distribution for the item qualities. This needs to be a Normal distribution. For this example, you can use the adjacency matrix provided with the package data and then use it to create a covariance matrix for the prior distribution of the item qualities. For your own data set, you will want to create an adjacency matrix with values that denote the similarity of different items to one another. For this example, the adjacency matrix ranks wards that are closer to one another as more similar. You can specify the outcome variable, the items involved in each comparison, and any prior information you have about the items' rating.  In this example, the items are the individual wards. The function will return an MCMC object containing samples from the posterior distribution of the item qualities. 

```{r covariance}
# Construct covariance matrix
expA <- expm::expm(forcedMarriage$adjacencyMatrix)
prior.var <- diag(diag(expA)^-0.5) %*% expA %*% diag(diag(expA)^-0.5)
```

Next, you can fit the model using the `speedyBBTm()` function. You will need to specify the outcome variable, the items involved in each comparison, and the prior covariance matrix for the item qualities. You can also specify the number of MCMC iterations to run. This example uses the default values for the hyperparameters of the prior distribution, but you can also specify your own values if you have prior knowledge about the items' rating using the `hyperparameter` argument. To match the plots in this vignette, you will need to use the same argument to `set.seed` as the one provided here. 

```{r fit_model}
# Fit model
set.seed(432)
forcedMarriageModel <- speedyBBTm(
  outcome = rep(1, length(forcedMarriage$comparisons$win)),
  player1 = forcedMarriage$comparisons$win,
  player2 = forcedMarriage$comparisons$lost,
  player.prior.var = prior.var,
  n.iter = 1000,
  burn.in = 100
)
```

To view some preliminary information about the results, you can calculate the means of the item quality parameters from the MCMC samples and summarise them. This will give you an idea of the estimated qualities of the items based on the comparisons in the dataset.

```{r view_output}
forcedMarriageModelMeans <- colMeans(forcedMarriageModel[, grep("lambda", varnames(forcedMarriageModel))])
forcedMarriageModelMeans
summary(forcedMarriageModel)
```

You can run some MCMC diagnostics to check convergence. The `speedyBBTm()` function returns an ["mcmc"](https://cran.r-project.org/package=coda) object that can be used with the `coda` package for diagnostics. You can use the `effectiveSize()` function to calculate the estimated effective sample size for convergence. Another example diagnostic illustrated here is the Geweke diagnostic, calculated using the `gewecke.diag()` function from `coda`. This diagnostic compares the difference between the two sample means in the first 10% and the last 50% of iterations in the Markov chain to assess convergence. For more information about available plots and diagnostics, you can refer to the [**coda** package documentation](https://cran.r-project.org/package=coda).

```{r some_diags}
effectiveSize(forcedMarriageModel)
geweke.diag(forcedMarriageModel)
``` 

You can use the `plot` function in `coda` or install `bayesplot` and `ggplot2` to get additional, colourful plots. For example, you can use the `mcmc_intervals()` function to create interval plots of the parameter posterior distributions and the `mcmc_dens()` function to create density plots of the MCMC samples. For more information, see the [**bayesplot** package documentation](https://mc-stan.org/bayesplot/).

```{r bayesplot_diagnostics, fig.width=8, fig.height=6, warning = FALSE, message = FALSE}
if (requireNamespace("bayesplot", quietly = TRUE)) {
  # Install bayesplot if not already installed
  # install.packages("bayesplot")
  # install.packages("ggplot2")
  library(bayesplot)
  library(ggplot2)
  # Trace plots
  bayesplot::mcmc_intervals(forcedMarriageModel, pars = c("lambda[10]", "lambda[20]", "lambda[30]", "lambda[40]")) +
    ggtitle("Posterior distributions with 95% credible intervals") +
    ylab(expression(lambda[i])) +
    xlab("Density")
  # Density plots
  bayesplot::mcmc_dens(forcedMarriageModel, pars = c("lambda[10]", "lambda[20]", "lambda[30]", "lambda[40]", "alpha.sq"))
} else {
  # Trace plots
  plot(forcedMarriageModel[, paste0("lambda[", c(10, 20, 30, 40), "]")], main = "Trace plots for item quality parameters")
}
```