Overview

I want to access the intercepts and coefficients for each level in a multilevel ordinal response model using the ordinal::clmm function in R.

I can easily do this with multilevel linear models estimated using lme4::lmer and calling the coef function. I can't seem to figure out how to do so with the ordinal model.

Example

Replication dataset

Below is a randomly-generated replication dataset. I create an independent variable ("indv"), a dependent variable ("depv") and an ordinal version of the dependent variable ("depv2") as well as some levels ("level").

test <- data.frame(depv = sample(1:4, 250, replace = TRUE),
                   indv = runif(250),
                   level = sample(1:4, 250, replace = TRUE)) %>%
  mutate(depv2 = factor(depv, levels = 1:4, labels = c("bad", "okay", "good", "great")),
         level = factor(level, levels = 1:4, labels = c("USA", "France", "China", "Brazil")))

Running the models

First, I estimate the linear model:

test1 <- lmer(depv ~ indv + (1 + indv | level), data = test)

Now, I estimate the ordinal response model:

test2 <- clmm(depv2 ~ indv + (1 + indv | level), data = test)

How to access the level coefficients?

I can easily access the level intercepts and coefficients for the linear model:

> coef(test1)
$level
       (Intercept)       indv
USA       2.239171  0.6238428
France    2.766888 -0.4173206
China     1.910860  1.2715852
Brazil    2.839156 -0.5599012

Doing this on the ordinal response model does not produce the same result:

> coef(test2)
   bad|okay   okay|good  good|great        indv 
-1.13105544  0.09101709  1.32240904  0.37157688 

For model test1 fitted by lme4, calling coef(test1) is internally doing lme4:::coef.merMod(test1). This is a user-friendly routine that adds fixed-effect coefficients and random-effect coefficients (conditional mode) together. Below is the source code of this nice function.

## source code of `lme4:::coef.merMod`
function (object, ...) 
{
    if (length(list(...))) 
        warning("arguments named \"", paste(names(list(...)), 
            collapse = ", "), "\" ignored")
    fef <- data.frame(rbind(fixef(object)), check.names = FALSE)  ## fixed-effect coefficients
    ref <- ranef(object, condVar = FALSE)  ## random-effect coefficients
    refnames <- unlist(lapply(ref, colnames))
    nmiss <- length(missnames <- setdiff(refnames, names(fef)))
    if (nmiss > 0) {
        fillvars <- setNames(data.frame(rbind(rep(0, nmiss))), 
            missnames)
        fef <- cbind(fillvars, fef)
    }
    val <- lapply(ref, function(x) fef[rep.int(1L, nrow(x)), 
        , drop = FALSE])
    for (i in seq(a = val)) {
        refi <- ref[[i]]
        row.names(val[[i]]) <- row.names(refi)
        nmsi <- colnames(refi)
        if (!all(nmsi %in% names(fef))) 
            stop("unable to align random and fixed effects")
        for (nm in nmsi) val[[i]][[nm]] <- val[[i]][[nm]] + refi[, 
            nm]
    }
    class(val) <- "coef.mer"
    val
}

In ordinal however, coef does not have this functionality. Instead, it simply extracts the $coefficients of the fitted model. So coef(test2) is as same as test2$coefficients. If you read ?clmm, you will see that this vector collects alpha (threshold parameters), beta (fixed-effect coefficients) and tau (if exists). Therefore, to get an output similar to what lme4 provides, we need to define the following function ourselves.

## inspired by `lme4:::coef.merMod`
coef_ordinal <- function (object, ...) 
{
    if (length(list(...))) 
        warning("arguments named \"", paste(names(list(...)), 
            collapse = ", "), "\" ignored")
    fef <- data.frame(rbind(object$beta), check.names = FALSE)  ## adapted
    ref <- ordinal::ranef(object, condVar = FALSE)  ## adapted
    refnames <- unlist(lapply(ref, colnames))
    nmiss <- length(missnames <- setdiff(refnames, names(fef)))
    if (nmiss > 0) {
        fillvars <- setNames(data.frame(rbind(rep(0, nmiss))), 
            missnames)
        fef <- cbind(fillvars, fef)
    }
    val <- lapply(ref, function(x) fef[rep.int(1L, nrow(x)), 
        , drop = FALSE])
    for (i in seq(a = val)) {
        refi <- ref[[i]]
        row.names(val[[i]]) <- row.names(refi)
        nmsi <- colnames(refi)
        if (!all(nmsi %in% names(fef))) 
            stop("unable to align random and fixed effects")
        for (nm in nmsi) val[[i]][[nm]] <- val[[i]][[nm]] + refi[, 
            nm]
    }
    #class(val) <- "coef.mer"  ## removed
    val
}

You can now call coef_ordinal(test2) for your desired output.