Ordered Beta Regression

ordbetareg_example.R

# Title: Ordered Beta Regression 

# Description: This code is to demonstrate using ordered beta regression. 
# This method can be particularly helpful for continuous data with upper/lower bounds. 

# For more information recommend reading Kubinec (2022)

# Setup ----

#... Packages ----

library(tidyverse) # ol faithful
library(brms) # fitting bayesian models 
library(tidybayes) # working with outputs and pltos
library(ordbetareg) # ordered beta regression
library(marginaleffects) # for computing marginal effects with model 

#... Custom Theme ----

#... Setting Theme ----

custom_theme <- function() {
  theme_minimal() %+replace% # basing this on the theme_minimal() function but editing some of the components
    theme(
      plot.title = element_text(hjust = 0.5, family = "Jost", face = "bold", size = 26),
      plot.subtitle = element_text(hjust = 0.5, size = 24),
      text = element_text(family = "Jost", size = 24)
    ) 
}

# Simulating Data ----

# - Sample size of 250
# - Binary treatment
# - Continuous outcome with bounds
# - One continuous confounder

ss = 250 # arbitrarily choosing a sample size of 250

df <- data.frame(
  z = rnorm(n = ss, mean = 5, sd = 2) # confounder
) %>% 
  dplyr::mutate(
    trt = rbinom(n = ss, size = 1, prob = plogis(0.03*z)),
    outcome = rordbeta(n = ss, mu = 0.1 + 0.2*trt + 0.05*z)
  )

# Plotting the Outcome ----

ggplot(data = df, 
       mapping = aes(x = outcome)
       ) + 
  tidybayes::stat_halfeye(fill = "#FFBE00", point_interval = NULL) +
  custom_theme() +
  labs(x = "Outcome", y = "Density", title = "Continuous Outcome with Bounds")

# Fitting OrdBetaReg ----

ordmod <- ordbetareg::ordbetareg(formula = outcome ~ trt + z, 
                                 data = df)

#... Posterior Predictive Check ----

p <- ordbetareg::pp_check_ordbeta(ordmod) 

p$discrete + custom_theme()
p$continuous + custom_theme()

# Results ----

# When interpretting the results, it can be useful to use the marginaleffects package. 
# This is especially useful when the model becomes more complicated. 

marginaleffects::avg_comparisons(ordmod)

# This can be compared to the above value of 0.2
# (When defining the outcome, mu = 0.2*trt + 0.05 *z)