r-在有序CLMM模型中查看每个级别的系数

Question

概述

我想使用R中的ordinal::clmm函数来访问多级有序响应模型中每个级别的截距和系数。

我可以很容易地通过使用lme4::lmer和调用coef函数估计的多级线性模型来实现这一点。我似乎不知道如何使用序数模型来做到这一点。

示例

复制数据集

下面是一个随机生成的复制数据集。我创建了一个自变量("indv")、一个因变量("depv")和因变量的有序版本("depv2")以及一些级别("级别")。

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")))

运行模型

首先，我估计线性模型：

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

现在，我估计有序响应模型：

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

如何访问级别系数

我可以很容易地访问线性模型的水平截距和系数：

> 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

在顺序响应模型上这样做不会产生相同的结果：

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

Answer 1

对于lme4安装的型号test1，调用coef(test1)就是在内部执行lme4:::coef.merMod(test1)。这是一个用户友好的例程，将固定效果系数和随机效果系数(条件模式)添加在一起。下面是这个漂亮函数的源代码。

## 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
}

但是，在序数中，coef不具有此功能。相反，它只是提取拟合模型的$coefficients。所以coef(test2)和test2$coefficients是一样的。如果您阅读?clmm，您将看到该向量收集alpha(阈值参数)、beta(固定效应系数)和tau(如果存在)。因此，要获得类似lme4提供的输出，我们需要自己定义以下函数。

## 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
}

现在，您可以调用coef_ordinal(test2)来获得所需的输出。