Comparison of ordbetareg and glmmTMB versions of the ordered beta regression model

compare_ordbeta_glmmTMB

# test glmmTMB's MLE implementation of the ordered beta regression model
# and compare it to the ordbetareg package, which uses
# Bayesian estimation (Stan/brms)

library(glmmTMB)
library(ordbetareg)
library(dplyr)

data("pew")

model_data <- select(pew,therm,age="F_AGECAT_FINAL",
                     sex="F_SEX_FINAL",
                     income="F_INCOME_FINAL",
                     ideology="F_IDEO_FINAL",
                     race="F_RACETHN_RECRUITMENT",
                     education="F_EDUCCAT2_FINAL",
                     region="F_CREGION_FINAL",
                     approval="POL1DT_W28",
                     born_again="F_BORN_FINAL",
                     relig="F_RELIG_FINAL",
                     news="NEWS_PLATFORMA_W28") %>%
  mutate_at(c("race","ideology","income","approval","sex","education","born_again","relig"), function(c) {
    factor(c, exclude=levels(c)[length(levels(c))])
  }) %>%
  # need to make these ordered factors for BRMS
  # normalize them for glmmTMB
  mutate(education=ordered(education),
         income=ordered(income),
         therm=(therm - min(therm))/(max(therm) - min(therm)))

ord_fit <- ordbetareg(therm ~ approval + (1|region),
                      data=model_data,
                      chains=1,
                      iter=1000,
                      backend="cmdstanr")

# NOTE: until this bug is fixed in their package, need to specify
# the starting values to the start option for the
# cutpoints, which they term psi

glmm_tmb_fit <- glmmTMB(therm ~ approval + (1|region),
                        data=model_data,
                        family=ordbeta(),
                        start=list(psi = c(-1, 1)))

# how close are they? pretty close

print(ord_fit,digits=5)

summary(glmm_tmb_fit)

# at least two significant digits... not bad...
# SEs are same to one significant digit

# check random effects

ranef(ord_fit)
ranef(glmm_tmb_fit)

# check cutpoints/phi

glmm_tmb_fit$fit$parfull
print(ord_fit,digits=5)

# lower cutpoint very close, upper cutpoint a bit off
# phi very close, although different sign (this is probably just an estimation difference)

# make nice table with modelsummary

library(modelsummary)

modelsummary(models=list(ordbetareg=ord_fit,
                         glmmmTMB=glmm_tmb_fit),
             statistic = "conf.int",metrics="RMSE",
             coef_rename=c("b_Intercept"="Intercept",
                           "(Intercept)"="Intercept",
                           "b_approvalDisapprove"="Disapprove",
                           "approvalDisapprove"="Disapprove",
                           'sd_region__Intercept'="Region SD",
                           "SD (Intercept)"="Region SD"))


# on the whole seems like a fantastic variant of ordbetareg,
# especially for getting estimates of predictors very quickly.
# at least if you are not overly concerned with the uncertainty/exact values of the random effects
# of course, does not have all brms features (missing data/latent variables etc.)