r - 为同一混合模型中的子集指定不同的随机结构？

Question

我想使用来自不同阻塞结构的不同实验的数据来做一个元模型。为此，我需要为同一模型中每个实验的数据指定不同的分块结构(随机效应结构)。 Genstat 有一个名为 vrmeta 的函数这样做(有关更多信息，请参阅 here)，但我更喜欢在 R 中工作，而且我不知道如何在 R 中进行。

例如，一个实验有块和主图，而另一个实验有块、主图和 split 图。我尝试为每个实验的块和图提供唯一的列，然后将模型编码为:

model <- lmer(response<-treatment1*treatment2*exp+
               (1|EXP1block/EXP1main)+
               (1|EXP2block/EXP2main/EXP2split),
           data=df)

这不起作用，我得到:

Error: Invalid grouping factor specification, EXP1main:EXP1block

...大概是因为 EXP2 的所有数据在 EXP1main 和 EXP1block 中都有 NA 值(反之亦然)。

如果有人能解释如何实现指定不同的结构，那就太好了。我目前正在使用包 lme4但如果这在不同的包中更容易，请告诉我。

如果需要，这是一些假数据的 dput 作为可重现的示例:

df<-structure(list(exp = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("EXP1", "EXP2"
), class = "factor"), treatment1 = structure(c(2L, 2L, 1L, 1L, 
2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L), .Label = c("N", 
"Y"), class = "factor"), treatment2 = c(40L, 60L, 40L, 60L, 40L, 
60L, 40L, 60L, 40L, 60L, 40L, 60L, 40L, 60L, 40L, 60L), response = c(780L, 
786L, 784L, 778L, 869L, 844L, 734L, 784L, 963L, 715L, 591L, 703L, 
925L, 720L, 642L, 678L), EXP1block = structure(c(1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, NA, NA, NA, NA, NA, NA, NA, NA), .Label = c("A", 
"B"), class = "factor"), EXP1main = c(1L, 2L, 3L, 4L, 1L, 2L, 
3L, 4L, NA, NA, NA, NA, NA, NA, NA, NA), EXP2block = structure(c(NA, 
NA, NA, NA, NA, NA, NA, NA, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L), .Label = c("A", 
"B"), class = "factor"), EXP2main = c(NA, NA, NA, NA, NA, NA, 
NA, NA, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L), EXP2split = structure(c(NA, 
NA, NA, NA, NA, NA, NA, NA, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("a", 
"b"), class = "factor")), class = "data.frame", row.names = c(NA, 
-16L))

Answer 1

这是使用 dummy() 的解决方案.

首先我们实际上必须替换NA值与 非北美 值(value)观；它们是什么并不重要，因为它们将被乘以零和/或被忽略......(可能有一个 tidyverse 和/或更简单的版本)

rep_nafac <- function(x,rval="other") {
    if (!any(is.na(x))) return(x)
    w <- which(is.na(x))
    old_lev <- levels(x)
    x <- as.character(x)
    x[is.na(x)] <- rval
    x <- factor(x,levels=c(old_lev,rval))
    return(x)
}
df_nona <- lapply(df,
                  function(x) if (!is.factor(x))
                                  replace(x,which(is.na(x)),1)
                  else rep_nafac(x))

现在用 dummy(exp,"level")+0 拟合模型作为每个分组变量的处理效果:这有效地将随机变量乘以指示变量，以判断观察是否在焦点组中。

library(lme4)
model <- lmer(response ~ treatment1*treatment2*exp+
               (dummy(exp,"EXP1")+0|EXP1main)+
               (dummy(exp,"EXP2")+0|EXP2main/EXP2split),
           data=df_nona)

结果看起来很合理:这是估计的方差。

Random effects:
 Groups             Name               Std.Dev.
 EXP2split:EXP2main dummy(exp, "EXP2") 33.361  
 EXP1main           dummy(exp, "EXP1")  7.706  
 EXP2main           dummy(exp, "EXP2") 33.271  
 Residual                              34.018