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

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，但如果使用不同的包更容易，请告诉我。

如果需要，这里是一些假数据的 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() 的解决方案。

• 首先我们实际上必须用 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