---
title: "Bayesian Beta Regression in JAGS"
author: "Ahmed Ali"
output:
  html_document:
    df_print: paged
header-includes: \usepackage{amsmath}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(dev = 'pdf')
```


## Notation:
$y_t$ = repsonse variable for observation t=1,..,N - should be in the range (0, 1).

$x_t$ = explanatory variable for obs t.

$\alpha, \beta$ = intercept and slope parameters to be estimated.

$\sigma$ = residual standard deviation.

## The model

\begin{align}

\ y_t \sim beta(a_t,b_t) \\

\mu_t = \frac{a_t}{a_t+b_t} \\

a_t = \mu_t*\phi \\

b_t = (1-\mu_t)*\phi \\

\end{align}

## Apply to dataset using JAGS

```{r}
# Load in
library(R2jags)
library(boot)
library(datasets) # Package which consists of a variety of datasets
head(attenu) # We will be applying the model to a dataset of peak acceleration measured at various observations for 33 earthquakes in California
```

```{r}
#Set up the data
acc=with(attenu,list(y=attenu$accel
                       ,x=attenu$dist
                       ,T=nrow(attenu)))
```

```{r}
plot(attenu$dist,attenu$accel)
```

```{r}
# Set up jags model

# Jags code to set up the model

model_code = '
model
{
  # Likelihood
  for (t in 1:T) {
    y[t] ~ dbeta(a[t], b[t])
    a[t] <- mu[t] * phi
    b[t] <- (1 - mu[t]) * phi
    logit(mu[t]) <- alpha + beta * x[t]
  }

  # Priors
  alpha ~ dnorm(0, 10^-2)
  beta ~ dnorm(0, 10^-2)
  phi ~ dunif(0, 10)
}
'
# Choose the parameters to watch
model_parameters =  c("alpha","beta","phi")

# Jags code to fit the model onto the dataset
jags_model=jags(acc,
                parameters.to.save = model_parameters
                ,model.file = textConnection(model_code),
                n.chains=4,
                n.iter=1000,
                n.burnin=200,
                n.thin=2)
```

```{r}
# Plot the jags output
print(jags_model)
```

