Beta regression using R and Stan

beta_regression.Rmd

---
title: "Beta regression"
output: html_notebook
---

## Setup

```{r setup}
library(ggplot2)
library(betareg)
library(cmdstanr)
library(posterior)
library(bayesplot)
```

### Data generation

$$
Y \sim \mathrm{Beta}(\alpha, \beta) \\
\mu = \frac{\alpha}{\alpha + \beta} = \mathrm{logit}^{-1}(-4.5 + 0.5 X) \\
\kappa = \alpha + \beta
$$

```{r}
set.seed(1234)
inv_logit <- function(x) {
  1 / (1 + exp(-x))
}
N <- 100
X <- runif(N, 0, 10)
logit_mu <- -4.5 + 0.7 * X
mu <-inv_logit(logit_mu)
kappa <- 6
alpha <- mu * kappa
beta <- (1 - mu) * kappa
Y <- rbeta(N, alpha, beta)
sim_data <- data.frame(X = X, Y = Y)
```

### View data

```{r}
p0 <- ggplot(sim_data, aes(x = X, y = Y)) +
  geom_point()
print(p0)
```

## Models

### Linear regression

```{r}
fit1 <- lm(Y ~ X, data = sim_data)
summary(fit1)
```

Plot results

```{r}
p0 +
  geom_abline(intercept = coef(fit1)[1],
              slope = coef(fit1)[2])
```

## Beta regression using the betareg package

```{r}
fit2 <- betareg(Y ~ X, sim_data, link = "logit")
summary(fit2)
```

Plot results

```{r}
p0 +
  geom_function(fun = function(x)
    inv_logit(coef(fit2)[1] + coef(fit2)[2] * x))
```

## Beta regression using Stan

```{r}
model_file <- file.path("models", "beta_regression.stan")
model <- cmdstan_model(model_file)
stan_data <- list(N = N, X = X, Y = Y)
fit3 <- model$sample(data = stan_data)
fit3$summary(variables = c("beta", "kappa"))
```

Plot results

```{r}
beta_mean <- fit3$summary("beta")$mean
p0 +
  geom_function(fun = function(x)
    inv_logit(beta_mean[1] + beta_mean[2] * x))
```

Posterior predictive check

```{r}
yrep <- fit3$draws("yrep") |>
  as_draws_matrix()
ppc_dens_overlay(y = Y, yrep = yrep[1:100, ])
```


## References

- Imad Ali, Jonah Gabry and Ben Goodrich (2020) Modeling Rates/Proportions using Beta Regression with rstanarm. https://mc-stan.org/rstanarm/articles/betareg.html

beta_regression.stan

data {
  int<lower=0> N;               // number of data points
  vector[N] X;                  // explanatory variable
  vector<lower=0,upper=1>[N] Y; // objective variable
}

parameters {
  array[2] real beta;  // intercept and slope (logit scale)
  real<lower=0> kappa; // precision parameter
}

transformed parameters {
  vector<lower=0,upper=1>[N] mu = inv_logit(beta[1] + beta[2] * X);
}

model {
  Y ~ beta_proportion(mu, kappa);

  // priors
  beta ~ normal(0, 10);
}

generated quantities {
  vector<lower=0,upper=1>[N] yrep;
  
  for (n in 1:N)
    yrep[n] = beta_proportion_rng(mu[n], kappa);
}