lmer 模型的标准化系数答案

【问题标题】：Standardized coefficients for lmer modellmer 模型的标准化系数
【发布时间】：2014-08-05 15:39:11
【问题描述】：

我曾经使用下面的代码来计算lmer 模型的标准化系数。但是，随着新版本的 lme 返回对象的结构发生了变化。

如何调整函数stdCoef.lmer 使其与新的lme4 版本一起使用？

# Install old version of lme 4
install.packages("lme4.0", type="both",
                 repos=c("http://lme4.r-forge.r-project.org/repos",
                         getOption("repos")[["CRAN"]]))

# Load package
detach("package:lme4", unload=TRUE)
library(lme4.0)

# Define function to get standardized coefficients from an lmer
# See: https://github.com/jebyrnes/ext-meta/blob/master/r/lmerMetaPrep.R
stdCoef.lmer <- function(object) {
  sdy <- sd(attr(object, "y"))
  sdx <- apply(attr(object, "X"), 2, sd)
  sc <- fixef(object)*sdx/sdy
  #mimic se.ranef from pacakge "arm"
  se.fixef <- function(obj) attr(summary(obj), "coefs")[,2]
  se <- se.fixef(object)*sdx/sdy
  return(list(stdcoef=sc, stdse=se))
}

# Run model
fm0 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

# Get standardized coefficients
stdCoef.lmer(fm0)

# Comparison model with prescaled variables
fm0.comparison <- lmer(scale(Reaction) ~ scale(Days) + (scale(Days) | Subject), sleepstudy)

【问题讨论】：

标签： multi-level mixed-models lmer coefficients

【解决方案1】：

@LeonardoBergamini 的答案有效，但这个更紧凑和易于理解，并且仅使用标准访问器 - 如果/当 summary() 输出的结构或拟合的内部结构时，将来不太可能中断模型，变化。

stdCoef.merMod <- function(object) {
  sdy <- sd(getME(object,"y"))
  sdx <- apply(getME(object,"X"), 2, sd)
  sc <- fixef(object)*sdx/sdy
  se.fixef <- coef(summary(object))[,"Std. Error"]
  se <- se.fixef*sdx/sdy
  return(data.frame(stdcoef=sc, stdse=se))
}
library("lme4")
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
fixef(fm1)
## (Intercept)        Days 
##   251.40510    10.46729
stdCoef.merMod(fm1)
##               stdcoef      stdse
## (Intercept) 0.0000000 0.00000000
## Days        0.5352302 0.07904178

（这确实给出了与stdCoef.lmer 相同的结果 @LeonardoBergamini 的回答...）

【讨论】：

【解决方案2】：

这应该可行：

stdCoef.lmer <- function(object) {
  sdy <- sd(attr(object, "resp")$y) # the y values are now in the 'y' slot 
  ###                                 of the resp attribute
  sdx <- apply(attr(object, "pp")$X, 2, sd) # And the X matriz is in the 'X' slot of the pp attr
  sc <- fixef(object)*sdx/sdy
  #mimic se.ranef from pacakge "arm"
  se.fixef <- function(obj) as.data.frame(summary(obj)[10])[,2] # last change - extracting 
  ##             the standard errors from the summary
  se <- se.fixef(object)*sdx/sdy
  return(data.frame(stdcoef=sc, stdse=se))
}

【讨论】：