时间戳

首页 > 解决方案 > 获取 glmer 模型的标准化系数？

r - 获取 glmer 模型的标准化系数？

问题描述

我被要求为glmer模型提供标准化系数，但不知道如何获得它们。不幸的是，该beta功能不适用于glmer模型：

Error in UseMethod("beta") : 
  no applicable method for 'beta' applied to an object of class "c('glmerMod', 'merMod')"

还有其他功能我可以使用，还是我必须自己写一个？

另一个问题是该模型包含几个连续预测变量（它们在相似的尺度上运行）和 2 个分类预测变量（一个有 4 个级别，一个有 6 个级别）。使用标准化系数的目的是将分类预测变量的影响与连续变量的影响进行比较，我不确定标准化系数是否适合这样做。标准化系数是一种可接受的方法吗？

模型如下：

model=glmer(cbind(nr_corr,maximum-nr_corr) ~ (condition|SUBJECT) + categorical_1 + categorical_2 + continuous_1 + continuous_2 + continuous_3 + continuous_4 + categorical_1:categorical_2 + categorical_1:continuous_3, data, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000)), family = binomial)

标签： rlme4

解决方案

reghelper::beta简单地标准化我们数据集中的数值变量。因此，假设您的分类变量是factors 而不是数字虚拟变量或其他对比编码，我们可以相当简单地将数据集中的数字变量标准化

vars <- grep('^continuous(.*)?', all.vars(formula(model)))
f <- function(var, data)
   scale(data[[var]])
data[, vars] <- lapply(vars, f, data = data)
update(model, data = data)

现在对于更一般的情况，我们可以或多或少地创建我们自己的beta.merMod函数。但是，我们需要考虑标准化是否有意义y。例如，如果我们有一个poisson模型，那么只有正整数值才有意义。另外一个问题变成了是否要缩放随机斜率效应，以及首先问这个问题是否有意义。在其中我假设分类变量被编码为characterorfactor而不是numericor integer。

beta.merMod <- function(model, 
                        x = TRUE, 
                        y = !family(model) %in% c('binomial', 'poisson'), 
                        ran_eff = FALSE, 
                        skip = NULL, 
                        ...){
  # Extract all names from the model formula
  vars <- all.vars(form <- formula(model))
  lhs <- all.vars(form[[2]])
  # Get random effects from the 
  ranef <- names(ranef(model))
  # Remove ranef and lhs from vars
  rhs <- vars[!vars %in% c(lhs, ranef)]
  # extract the data used for the model
  env <- environment(form)
  call <- getCall(model)
  data <- get(dname <- as.character(call$data), envir = env)
  # standardize the dataset
  vars <- character()
  if(isTRUE(x))
    vars <- c(vars, rhs)
  if(isTRUE(y))
    vars <- c(vars, lhs)
  if(isTRUE(ran_eff))
    vars <- c(vars, ranef)
  data[, vars] <- lapply(vars, function(var){
    if(is.numeric(data[[var]]))
      data[[var]] <- scale(data[[var]])
    data[[var]]
  })
  # Update the model and change the data into the new data.
  update(model, data = data)
}

该函数适用于线性和广义线性混合效应模型（未针对非线性模型进行测试），并且与其他 beta 函数一样使用reghelper

library(reghelper)
library(lme4)
# Linear mixed effect model
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
fm2 <- beta(fm1)
fixef(fm1) - fixef(fm2)
(Intercept)        Days 
  -47.10279   -19.68157 

# Generalized mixed effect model
data(cbpp)
# create numeric variable correlated with period
cbpp$nv <- 
  rnorm(nrow(cbpp), mean = as.numeric(levels(cbpp$period))[as.numeric(cbpp$period)])
gm1 <- glmer(cbind(incidence, size - incidence) ~ nv + (1 | herd),
              family = binomial, data = cbpp)
gm2 <- beta(gm1)
fixef(gm1) - fixef(gm2)
(Intercept)          nv 
  0.5946322   0.1401114

但是请注意，与beta该函数不同的是，该函数返回更新后的模型，而不是模型的摘要。

另一个问题是该模型包含几个连续预测变量（它们在相似的尺度上运行）和 2 个分类预测变量（一个有 4 个级别，一个有 6 个级别）。使用标准化系数的目的是将分类预测变量的影响与连续变量的影响进行比较，我不确定标准化系数是否适合这样做。标准化系数是一种可接受的方法吗？

现在这是一个很好的问题，而且更适合stats.stackexchange，但我不能确定答案。

推荐阅读