Hereâs an example of how you could approach sum-score contrasts
# Required packages
Pkgs <- c("tidyverse"
, "tidybayes"
, "brms"
, "parallel" )
## Load packages
lapply(Pkgs, require, c = T)
## Set computing options
options(mc.cores = detectCores())
ncores = detectCores()
# Make data
## Ordinal function
get_probs <- function(eta,cuts){
p = rep(NA,(length(cuts)+1))
p[1] = 1-plogis(eta-cuts[1])
for(i in 2:length(cuts)){
p[i] = plogis(eta-cuts[i-1]) - plogis(eta-cuts[i])
}
p[length(p)] = plogis(eta-cuts[i])
plot = ggplot(
data = data.frame(x=1:length(p),y=p)
, mapping = aes(x=x,y=y)
) +
geom_bar(stat='identity',position='identity')+
scale_y_continuous(limits=c(0,1),expand=c(0,0))
print(plot)
return(p)
}
## Create data for two groups
p1 <- get_probs(-1, seq(-3, 3, length.out = 6))
p2 <- get_probs(1, seq(-3, 3, length.out = 6))
## Make a data frame
set.seed(1)
data <- data.frame(
sex = rep(c('male','female'), each = 30)
, response = c(
sample(1:7,30, prob = p1, replace = TRUE)
, sample(1:7,30, prob = p2, replace = TRUE)
)
) %>%
mutate(sex = factor(sex))
## Look at the simulated data
data %>%
ggplot(aes(x = response)) +
geom_histogram(aes(fill = sex)
, bins = 7
, position = position_dodge(width = 1)) +
scale_x_continuous(breaks = 1:7) +
theme_classic()
# Model
## Formula
fmla <- bf(response ~ sex)
## Priors
priors <- c(set_prior("normal(0,5)"
, class = "Intercept"
)
, set_prior("normal(0,1)"
, class = "b"
)
)
## Fit the model
fit <- brm(fmla
, data
, family = cumulative()
, prior = priors
, iter = 2000
, warmup = 1000
, chains = 4
, cores = ncores/2
, backend = "cmdstanr"
, threads = threading(2)
, normalize = FALSE
, save_pars = save_pars(all = TRUE)
, stan_model_args=list(stanc_options = list("O1"))
, control = list(adapt_delta = 0.99
, max_treedepth = 14)
)
# Analysis
## Posterior predictive check
fit$data %>%
add_epred_draws(fit
, ndraws = 100) %>%
group_by(sex, .category, .draw) %>%
summarise(count = sum(.epred)) %>%
group_by(sex, .category) %>%
mean_qi(count) %>%
rename(response = .category) %>%
ggplot(aes(x = as.numeric(response))) +
geom_bar(data = data,
fill = "grey70") +
facet_wrap(~sex) +
geom_pointrange(aes(y = count, ymin = .lower, ymax = .upper),
size = 1/6)
## Get sum-score intervals
intervals <- data %>%
add_epred_draws(fit, re_formula = NULL, ndraws = 100) %>%
ungroup() %>%
mutate(product = as.double(.category) * .epred
) %>%
group_by(sex, .row, .draw) %>%
summarise(sum_score = sum(product)) %>%
ungroup() %>%
group_by(sex) %>%
point_interval(
sum_score
, .width = 0.89
, .point = mean
, .interval = hdci
, .simple_names = TRUE
, na.rm = TRUE
, .exclude = c(".width", ".point", ".interval")
) %>%
ungroup()
## Plot sum-score intervals
intervals %>%
select(-.width, -.point, -.interval) %>%
ggplot(aes(sex, sum_score)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = .lower, ymax = .upper),
width = 0.1) +
scale_x_discrete("Sex") +
scale_y_continuous("Rating (Sum-Score)",
limits = c(1,7),
breaks = c(1:7)) +
theme_classic()
## Calculate sum-score difference
contrast <- data %>%
add_epred_draws(fit, re_formula = NULL, ndraws = 100) %>%
ungroup() %>%
mutate(product = as.double(.category) * .epred
) %>%
group_by(sex, .row, .draw) %>%
summarise(sum_score = sum(product)) %>%
ungroup() %>%
select(sex, sum_score, .draw) %>%
dplyr::group_by(sex) %>%
mutate(row = row_number()) %>%
group_by(row, .draw) %>%
compare_levels(sum_score
, by = "sex") %>%
ungroup() %>%
group_by(sex) %>%
point_interval(
sum_score
, .width = 0.89
, .point = mean
, .interval = hdci
, .simple_names = TRUE
, na.rm = TRUE
, .exclude = c(".width", ".point", ".interval")
) %>%
ungroup()
## Plot sum-score difference
contrast %>%
select(-.width, -.point, -.interval) %>%
ggplot(aes(sex, sum_score)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = .lower, ymax = .upper),
width = 0.1) +
scale_x_discrete("Male - Female") +
scale_y_continuous("Rating (Sum-Score) Difference",
limits = c(-3,3),
breaks = c(-3:3)) +
geom_hline(yintercept = 0, linetype = 2) +
theme_classic()