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

Posted 2023-03-12

【中文标题】获取 glmer 模型的标准化系数？

【英文标题】：Getting standardized coefficients for a glmer model?

【发布时间】：2021-01-12 11:49:21

【问题描述】：

有人要求我为glmer 模型提供标准化系数，但我不确定如何获得它们。不幸的是，beta 函数不适用于glmer 模型：

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

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

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

型号如下：

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)

【参考方案1】：

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 模型，那么只有正整数值才有意义。另外一个问题变成了是否要缩放随机斜率效应，以及首先问这个问题是否有意义。在其中我假设分类变量被编码为character 或factor 而不是numeric 或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)

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

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 个级别，一个有六个级别）。使用标准化系数的目的是将分类预测变量的影响与连续变量的影响进行比较，我不确定标准化系数是否适合这样做。标准化系数是一种可接受的方法吗？

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

【讨论】：

非常感谢奥利弗！我仍在努力解决这一切。请问'cbpp$nv

你做对了。我使用的示例数据cbpp 没有数字变量，我手头也没有示例。但它确实有一个因素period（值：“1”到“4”）。 period 在?glmer 的示例中用作固定效应，因此我通过生成一些平均值在 1 和 4 之间（相当于周期）的随机数据来为该因子创建了一个数字代理。简而言之，它是rnorm(nrow(cbpp), mean = (1:4)[match with period as numeric])。 rnorm(n, mean = mu) 只是生成随机的 n 数字，其手段为 mu（mu 可能是更多的数字或长度为 n 的向量）

我也可以使用cbpp$nv <- rnorm(nrow(cbpp), mean = as.numeric(cbpp$period))。

【参考方案2】：

再次感谢你，奥利弗！对于任何对我的问题最后一部分的答案感兴趣的人，

另一个问题是模型包含几个连续的 预测变量（在相似的尺度上运行）和 2 个分类 预测变量（一个有 4 个级别，一个有六个级别）。的目的 使用标准化系数将比较 分类预测器对那些连续预测器，我是 不确定标准化系数是否合适 所以。标准化系数是一种可接受的方法吗？

你可以找到答案here。 tl;博士是，无论如何，使用标准化回归系数并不是混合模型的最佳方法，更不用说像我这样的模型了......