Ordinal Bayesian Power

gistfile1.txt

# Load packages
library(tidyverse)
library(brms)
library(broom)

# We need to simulate data from an ordinal distribution. This is hard.
# Unlike many distributions, treatment effects can affect the mean *and*
# the sd of the latent y* distribution that underlies the ordinal distribution
# (see Liddell & Kruschke (2018) and Bürkner & Vuorre (2019))

# Create function to simulate 5-point Likert scale distributions
# t = thresholds
# b_T = treatment effect on latent mean
# d_t = treatment effect on latent sd

r_likert <- function(seed, n, t = c(-.5, 0, .5, 1), b_T = c(0, .8), d_T = c(0, .5)) {
  
  # Set random seed
  set.seed(seed)
  
  # Simulate random variables
  d <- tibble(treatment = rep(c("control", "treatment"), each = n) %>% factor())
  
  # Compute latent-y mean and sd parameters for each individual
  d <- 
    d %>% 
    mutate(
      mean = b_T[treatment],
      sd = 1/exp(d_T[treatment])
    )
  
  # Calculate probabilities of each response for each individual
  d <- 
    d %>% 
    mutate(
      p1 = pnorm(t[1], mean, sd),
      p2 = pnorm(t[2], mean, sd) - pnorm(t[1], mean, sd),
      p3 = pnorm(t[3], mean, sd) - pnorm(t[2], mean, sd),
      p4 = pnorm(t[4], mean, sd) - pnorm(t[3], mean, sd),
      p5 = 1 - pnorm(t[4], mean, sd)
    )
    
  
  # For each individual, sample a single response category
  d <- d %>% mutate(y = NA)
  for (i in 1:nrow(d)) {
    d$y[i] <- sample(
      x = c(1:5), 
      size = 1, 
      prob = c(d$p1[i], d$p2[i], d$p3[i], d$p4[i], d$p5[i])
    )
  }

  # Return data frame
  d %>% select(treatment, y)
  
}

# Simulate data
test <- r_likert(seed = 1,
                 n = 50,
                 t = c(-.5, 0, .5, 1),
                 b_T = c(0, .8),
                 d_T = c(0, .5))

# Fit model
fit <- 
  brm(formula = 
        bf(y ~ treatment) +
        lf(disc ~ 0 + treatment, cmc = F),
      family = cumulative(link = "probit",
                          link_disc = "log"),
      prior = 
        prior(normal(0, 1.5), class = "Intercept") + 
        prior(normal(0, .5), class = "b") +
        prior(normal(0, 1), class = "b", dpar = "disc"),
      data = test,
      iter = 2e3,
      chains = 2,
      cores = 2,
      seed = 1)
      
      
# Now that we've fit the model, we can do the same power calculations as
# any other type of model. Though note, we have estimates for *two*
# treatment effects

# Create function to simulate Likert data and fit model
sim_and_fit <- function(seed, n){
  
  # Simulate data
  d <- r_likert(seed = seed, n = n)
  
  # Fit model
  update(fit,
         newdata = d, 
         seed = seed,
         iter = 2e3,
         chains = 4,
         cores = 4,
         refresh = 0) %>% 
    tidy(prob = .95) %>% 
    rename(error = std.error) %>% 
    filter(term %in% c("b_treatmenttreatment", "b_disc_treatmenttreatment"))
}

# Simulate 100 times (this can be memory intensive)
start <- Sys.time()
s1 <-
  tibble(seed = 1:100) %>% 
  mutate(tidy = map(seed, sim_and_fit, n = 100)) %>% 
  unnest(tidy)
end <- Sys.time()
end - start

# Power Calculations - latent mean
s1 %>% 
  filter(term == "b_treatmenttreatment") %>% 
  ggplot(aes(x = seed, y = estimate, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = c(0, .8), color = "grey") +
  geom_pointrange(fatten = 1/2) +
  labs(x = "seed (i.e., simulation index)",
       y = expression(beta[1])) +
  theme_minimal() +
  theme(panel.grid = element_blank())

s1 %>% 
  filter(term == "b_treatmenttreatment") %>% 
  mutate(check = ifelse(lower > 0, 1, 0)) %>% 
  summarise(power = mean(check))
  
# Power calculations - latent sd
s1 %>% 
  filter(term == "b_disc_treatmenttreatment") %>% 
  ggplot(aes(x = seed, y = estimate, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = c(0, .5), color = "grey") +
  geom_pointrange(fatten = 1/2) +
  labs(x = "seed (i.e., simulation index)",
       y = expression(beta[1])) +
  theme_minimal() +
  theme(panel.grid = element_blank())

s1 %>% 
  filter(term == "b_disc_treatmenttreatment") %>% 
  mutate(check = ifelse(lower > 0, 1, 0)) %>% 
  summarise(power = mean(check))