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

发布于 2025-02-11 11:59:37 字数 8 浏览 1 评论 0原文

continue

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

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

遮了一弯 2025-02-18 11:59:37

对于 lme4 安装的模型test1，呼叫coef（test1）是在内部进行lme4 ::: coef.mermod（test1）lme4 ::: /代码>。这是一种用户友好的例程，可以将固定效应系数和随机效应系数（条件模式）一起添加在一起。以下是此不错功能的源代码。

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

但是，在 ordinal 中，coef没有此功能。相反，它只需提取拟合模型的$系数即可。因此，coef（test2）与test2 $系数一样。如果您阅读？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）。

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.

回复收藏 0 原文

~没有更多了~