为同一混合模型内的子集指定不同的随机结构吗？

Question

我想使用来自具有不同阻断结构的不同实验的数据来做一个元模型。为此，我需要为同一模型中每个实验的数据指定不同的阻止结构（随机效应结构）。 Genstat具有一个称为vrmeta的函数来执行此操作（有关更多信息，请参见here），但是我更喜欢在R中工作，但我不知道如何在R中进行操作。

例如，一个实验有块状图和主图，而另一个实验有块状图，主图和分割图。我尝试为每个实验的块和图赋予唯一的列，然后将模型编码为：

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

这不起作用，我得到：

错误：无效的分组因子说明，EXP1main：EXP1block

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

如果有人能解释如何指定不同的结构，那就太好了。我当前正在使用软件包lme4，但如果在其他软件包中更方便使用，请告诉我。

这里是一些虚假数据的副本，如果需要，可以作为可重现的示例：

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()的解决方案。

• 首先，我们实际上必须用non-NA值替换NA值；无关紧要，因为它们将乘以零和/或被忽略...（可能有一个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