Simulate Ordinal Mixed Model Coefficients and Make Posterior Predictions

1-clmm-sim-example.R

# Let's load the packages and the data. -----
# See this instead: http://svmiller.com/blog/2020/04/summarizing-ordinal-models-with-simulation-multivariate-normal/
library(stevemisc)
library(post8000r)
library(tidyverse)


TV16 %>%
  filter(racef == "White") %>%
  filter(state %in% c("Indiana","Ohio","Pennsylvania","Wisconsin","Michigan")) %>%
  mutate(whiteadvf = ordered(whiteadv)) -> Data


Data %>%
  mutate_at(vars("age", "famincr","pid7na","ideo","lcograc","lemprac"),
            list(z = ~r2sd(.))) %>%
  rename_at(vars(contains("_z")),
            ~paste("z", gsub("_z", "", .), sep = "_") ) -> Data


# Simple ordinal model with state random effects. ----
require(ordinal)
Ord1 <- clmm(whiteadvf ~ z_age + z_famincr + z_pid7na + female + collegeed  + z_ideo + (1 | state),
             data=Data)

# Extract some model parameters and make sims ----

coefOrd1 <- coef(Ord1)
vcovOrd1 <- vcov(Ord1) # vcov without state random effect
vcovOrd1 <- vcovOrd1[-nrow(vcovOrd1), -ncol(vcovOrd1)]

set.seed(8675309)
simOrd1 <- MASS::mvrnorm(1000, coefOrd1, vcovOrd1) %>% tbl_df() %>%
  mutate(sim = seq(1:1000)) %>% select(sim, everything())


# Get hypothetical data.
require(modelr)
Data %>%
  # Let's keep it simple and obvious. The min of z_ideo is "very liberal."
  # The max of z_ideo is "very conservative."
  # We should expect to see magnitude differences.
  data_grid(.model = Ord1, z_ideo = c(min(z_ideo, na.rm=T), max(z_ideo, na.rm=T))) %>%
  # because we got 1000 sims, we need to repeat this 1000 times
  dplyr::slice(rep(row_number(), 1000)) ->  newdatOrd1


simOrd1 %>%
  # repeat it twice because we have two values of z_ideo
  dplyr::slice(rep(row_number(), 2)) %>%
  # arrange by simulation number after we repeated the data twice
  arrange(sim) %>%
  # rename these to be clear they're simulated coefficients
  rename_at(vars("z_age", "z_famincr", "z_pid7na", "female", "collegeed", "z_ideo"),
            ~paste0("coef", .)) %>%
  bind_cols(., newdatOrd1) %>%
  # Manually calculate the logits
  mutate(logit1 = `1|2` - ((z_age*coefz_age) + (z_famincr*coefz_famincr) + (z_pid7na*coefz_pid7na) +
           (female*coeffemale) + (collegeed*coefcollegeed) + (z_ideo*coefz_ideo)),
         logit2 = `2|3` - ((z_age*coefz_age) + (z_famincr*coefz_famincr) + (z_pid7na*coefz_pid7na) +
           (female*coeffemale) + (collegeed*coefcollegeed) + (z_ideo*coefz_ideo)),
         logit3 = `3|4` - ((z_age*coefz_age) + (z_famincr*coefz_famincr) + (z_pid7na*coefz_pid7na) +
           (female*coeffemale) + (collegeed*coefcollegeed) + (z_ideo*coefz_ideo)),
         logit4 = `4|5` - ((z_age*coefz_age) + (z_famincr*coefz_famincr) + (z_pid7na*coefz_pid7na) +
           (female*coeffemale) + (collegeed*coefcollegeed) + (z_ideo*coefz_ideo)))  %>%
  mutate_at(vars(contains("logit")), list(p = ~plogis(.))) %>%
  mutate(p1 = logit1_p,
         p2 = logit2_p - logit1_p,
         p3 = logit3_p - logit2_p,
         p4 = logit4_p - logit3_p,
         p5 = 1 - logit4_p,
         sump = p1 + p2 + p3 + p4 + p5) -> simOrd1


simOrd1 %>%
  group_by(z_ideo) %>%
  summarize_at(vars("p1","p2","p3","p4","p5"),
               list(mean = ~mean(.),
                    lwr = ~quantile(., .025),
                    upr = ~quantile(., .975)))  %>%
  select(z_ideo, p5_mean, p5_lwr, p5_upr)