具有固定和随机效应的库（lme4）中的 MLM 设计矩阵答案

【问题标题】：Design matrix for MLM from library(lme4) with fixed and random effects具有固定和随机效应的库（lme4）中的 MLM 设计矩阵
【发布时间】：2014-10-21 16:33:43
【问题描述】：

应用上下文
我有一个随机斜率和截距的模型。有许多级别的随机效应。新数据（待预测）可能具有也可能不具有所有这些级别。

为了更具体地说明这一点，我正在处理专辑级别的音乐收入 (title)。每张专辑可能有多种类型format2（CD、黑胶唱片、电子音频等）。我对每种专辑的每张专辑的收入都有衡量标准。模型指定为：

lmer(physical~ format2+ (0+format2|title))

问题是未来的数据可能没有title 或format2 的所有级别。对于随机拦截，这很容易通过predict(..., allow.new.levels= TRUE) 解决。但是对于固定效应和随机斜率是有问题的。因此，我正在尝试编写一个函数来灵活预测merMod 对象，类似于lme4::predict.merMod；但这将处理训练数据和预测数据之间的差异。出于对lme4::predict.merMod 的确切细节的无知，这是一个与其他问题一样多的问题。

问题描述
问题的症结在于获得正确的model.matrix()，具有固定和随机效应来计算预测和 SE。 merMod 类的 S3 方法仅返回固定效果。

stats::model.matrix() 基本函数的文档非常有限。不幸的是，我不拥有Statistical Models in S 或Software for Data Analysis，它们似乎拥有这些功能背后的细节。

model.matrix() 应该采用模型公式和新数据框并生成设计矩阵。但我遇到了一个错误。您能提供的任何帮助将不胜感激。

示例数据

dat1 <- structure(list(dt_scale = c(16, 16, 16, 16, 16, 16, 16, 16, 16, 
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16), title = c("Bahia", 
"Jazz Moods: Brazilian Romance", "Quintessence", "Amadeus: The Complete Soundtrack Recording (Bicentennial Edition)", 
"Live In Europe", "We'll Play The Blues For You", "The Complete Village Vanguard Recordings, 1961", 
"The Isaac Hayes Movement", "Jazz Moods: Jazz At Week's End", 
"Blue In Green: The Concert In Canada", "The English Patient - Original Motion Picture Soundtrack", 
"The Unique Thelonious Monk", "Since We Met", "You're Gonna Hear From Me", 
"The Colors Of Latin Jazz: Cubop!", "The Colors Of Latin Jazz: Samba!", 
"Homecoming", "Consecration: The Final Recordings Part 2 - Live At Keystone Korner, September 1980", "More Creedence Gold", "The Stardust Session"), format2 = c("CD", "CD", 
"CD", "CD", "CD", "CD", "CD", "SuperAudio", "SuperAudio", "CD", "E Audio", "CD", 
"Vinyl", "CD", "E Audio", "CD", "CD", "CD", "CD", "CD"), mf_day = c(TRUE, 
TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, 
TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE), xmas = c(FALSE, 
FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, 
FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, 
FALSE), vday = c(FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, 
FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, 
FALSE, FALSE, FALSE, FALSE), yr_since_rel = c(16.9050969937038, 
8.41815617876864, 9.2991404674865, 25.0870296783559, 39.1267038232812, 
27.9156764326061, 9.11596751812513, 23.3052837112449, 14.3123922258974, 
30.5208152866414, 5.83025071417496, 21.3090003877291, 7.75022155568392, 
11.3601605287827, 0.849006673421519, 31.9918631305662, 13.8861905547041, 
12.8342695062012, 29.6916671402534, 13.5912612705038), physical = c(1327.17849171096, 
-110.2265302258, -795.37376268564, 355.06192702004, -1357.3492884345, 
-1254.93442612023, -816.713683621225, 881.201935773452, -3092.02845691036, 
-2268.6304275652, 907.347941142021, -699.130275178185, 377.867849132077, 
-1047.50531157311, 1460.25978951805, 1376.84579069304, 3619.03629114089, 
962.888173535704, 2514.77880599199, 2539.14958588771)), .Names = c("dt_scale", 
"title", "format2", "mf_day", "xmas", "vday", "yr_since_rel", 
"physical"), row.names = c(1L, 2L, 5L, 6L, 7L, 8L, 9L, 11L, 12L, 
13L, 14L, 15L, 20L, 22L, 23L, 25L, 27L, 32L, 35L, 36L), class = "data.frame")

公式：

f1 <- as.formula(~1 + dt_scale + yr_since_rel + format2 + (0 + format2 + mf_day + 
xmas + vday | title))

执行/错误

library(lme4)
model.matrix(f1, data= dat1)
Error in 0 + format2 : non-numeric argument to binary operator

注意我也用Orthodont 数据试过这个；但是，我得到一个不同的错误。

library(lme4)
data("Orthodont",package="MEMSS")
fm1 <- lmer(formula = distance ~ age*Sex + (1+age|Subject), data = Orthodont)
newdat <- expand.grid(
  age=c(8,10,12,14)
  , Sex=c("Male","Female")
  , distance = 0
  , Subject= c("F01", "F02")
)


f1 <- formula(fm1)[-2] # simpler code via Ben Bolker below
mm <- model.matrix(f1, newdat) # attempt to use model.matrix
Warning message
In Ops.factor(1 + age, Subject) : | not meaningful for factors

# use lme4:::mkNewReTrms as suggested in comments
mm <- lme4:::mkNewReTrms(f1, newdat) 
Error in lme4:::mkNewReTrms(f1, newdat) : object 'ReTrms' not found
In addition: Warning message:
In Ops.factor(1 + age, Subject) : | not meaningful for factors

# check if different syntax would fix this
mm <- lme4::mkNewReTrms(f1, newdat)
Error: 'mkNewReTrms' is not an exported object from 'namespace:lme4'
mm <- mkNewReTrms(f1, newdat)
Error: could not find function "mkNewReTrms"

【问题讨论】：

我有几个问题/cmets。（1）您包含的样本数据只有一个值format2，因此指定的模型不起作用（大概您的真实数据有更多）。 (2) 您的示例不可重复地继续； b1 和 a 是什么？ (3) 你的f1 公式看起来很可疑；所有这些影响是否在title 的范围内变化，您是否有足够的数据来估计所有这些影响的（相关）标题之间的可变性？ (4) 用于预测的新数据中固定或随机效应水平的缺失水平不成问题；麻烦的是额外的关卡
如果你想构建一个新的随机效应模型矩阵，你可以使用lme4:::mkNewReTrms(object,newdata,re.form)，其中object是一个公式；然后提取并转置结果对象的$Zt 组件
formula(fm1)[-2] 会让您更轻松地获得公式的 RHS。最后，您能否展示一个（可重现的）示例，说明您的预测数据没有以您关注的方式工作？正如我在评论 #1 第 4 部分中所说，我不相信/不明白使用内置 predict 函数做你想做的事情是否有任何困难。
@BenBolker -- 更新 (1) 和 (2)。再（3），数据不平衡。因此，标题之间的随机斜率项存在差异，但不一定每个标题都有大量案例。 RE (4)，我无法提供完整的数据（专有）。但是我所有多余的级别都已清理干净。我会调查lme4:::mkNewReTrms()。
“如何为模型的随机效应部分构建模型矩阵（通常表示为 Z）？”是一个合理的问题，但我仍然不相信您提供了任何实际需要这样做的示例。理想情况下，您会生成一个小/假的可重现示例（例如，使用Orthodont 数据）显示您正在尝试做什么（您的最终目标，即预测新的和/或缺失的随机和/或固定效应水平）跨度>

标签： r lme4 s

【解决方案1】：

编辑于 2015 年 8 月 12 日：见 changes on Github 和 GitHub Repo

已编辑，2014 年 10 月 15 日：这个答案还不完美。仍然有几个有错误的用例（见下面的评论链）。但它在大多数情况下都有效。我会在某个时候完成它。

我相信这个功能将解决更重要的问题，准确预测 merMod 对象。 Bolker 博士，这里仍然存在一些问题（例如稀疏性和效率）；但我相信该方法有效：

data("Orthodont",package="MEMSS")
fm1 <- lmer(formula = distance ~ age*Sex + (1+age|Subject), data = Orthodont)
newdat <- expand.grid(
  age=c(8,10,12,14)
  , Sex=c("Male","Female")
  , distance = 0
  , Subject= c("F01", "F02")
)

predict.merMod2 <- function(object, newdat=NULL) {
# 01. get formula and build model matrix
  # current problem--model matrix is not sparse, as would be ideal
  f1 <- formula(object)[-2]
  z.fe <- model.matrix(terms(object), newdat)
  z.re <- t(lme4:::mkReTrms(findbars(f1), newdat)$Zt)
  mm <- cbind(z.fe, 
              matrix(z.re, nrow= dim(z.re)[1], ncol= dim(z.re)[2],
                     dimnames= dimnames(z.re)))

  # 02. extract random effect coefficients needed for the new data
  # (a) - determine number of coef
  len <- length(ranef(object)) 
  re.grp.len <- vector(mode= "integer", length= len) 
  for (i in 1:len) { # for each random group
    re.grp.len[i] <- dim(ranef(object)[[i]])[2] # number of columns (slope and intercept terms)
  }

  # (b) - create beta vector
  fe.names <- unique(colnames(mm)[1:length(fixef(object)) - 1])
  re.names <- unique(colnames(mm)[-c(1:length(fixef(object)) - 1)]) 
  beta.re <- as.vector(rep(NA, length= sum(re.grp.len) * length(re.names)), mode= "numeric")
  for (i in 1:len) {
    re.beta  <- ranef(object)[[i]][rownames(ranef(object)[[i]]) %in% re.names,] 
    ind.i <- sum(!is.na(beta.re)) + 1; ind.j <- length(as.vector(t(re.beta))) 
    beta.re[ind.i:ind.j]  <- as.vector(t(re.beta)) 
  }
  beta <- c(fixef(object)[names(fixef(object)) %in% fe.names], beta.re)
  # 03. execute prediction
  return(mm %*% beta)
}

predict.merMod2(fm1, newdat)

【讨论】：

我在我的 5 个模型上对音乐数据进行了测试（如最初所述）。它适用于 5 个中的 4 个（即识别错误）。我会尽快更新这个答案以解决这个问题，可能是明天。
我想我已经解决了 lme4 中的问题（参见 github.com/lme4/lme4/commit/… ）。您愿意尝试从 Github (devtools::install_github("lme4","lme4")) 重新安装并尝试您的测试用例吗？你需要编译器等——如果这是一个问题，请联系我，我会为你构建一个二进制文件......
@BenBolker 谢谢——我下周假期回来时会检查一下。
@Alex : 你已经答应@BenBolker 检查lme4 的开发版本（可从 github 获得）......但在你假期后可能忘记了这一点？
@MartinMächler 开发版本（当时的1.1.8）无法满足我的需求，这就是为什么我分叉了 lme4/lme4 存储库并根据我们之前的电子邮件更新了predict.merMod。