specify different random structures for subsets within the same mixed model?

Question

I would like to do a meta-model using data from different experiments with different blocking structures. For this, I would need to specify the different blocking structure (random effects structure) for the data from each experiment within the same model. Genstat has a function called vrmeta that does this (see here for more info) but I prefer to work in R, and I can't figure out how to do it in R.

For example, one experiment has blocks and main plots, while another has blocks, main plots and split plots. I have tried giving each experiment unique columns for its blocks and plots, and then coding the model as:

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

This does not work and I get:

Error: Invalid grouping factor specification, EXP1main:EXP1block

... presumably because all data for EXP2 has NA values in EXP1main and EXP1block (and vice versa).

It would be great if someone could explain how specifying different structures could be achieved. I am currently using package lme4 but if this is easier in a different package please let me know.

Here is a dput of some fake data as a reproducible example if needed:

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))

score 2 · Accepted Answer · answered Nov 26 '19 at 02:10

Here's a solution using dummy().

First we actually have to replace the NA values with non-NA values; it doesn't matter what they are since they will be multiplied by zero and/or ignored ... (there may be a tidyverse and/or simpler version of this)

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))

Now fit the model with dummy(exp,"level")+0 as the treatment effect for each grouping variable: this effectively multiplies the random variables by an indicator variable for whether the observation is in the focal group or not.

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

The results look sensible: here are the estimated variances.

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

Thank you for this! I have checked it on my full dataset as well as the data in the dput above, and it works fine. I also compared this with putting the same dataset through Genstat's vrmeta and I get the same output for the random effect variance and the fixed effect estimates, so it must be doing the same thing. Just for the sake of consistency with my original question, I would include EXP1block and EXP2block in the model with the `dummy` variables - but this makes no difference to how the code is structured :) — corn_bunting, Nov 28 '19 at 11:53

specify different random structures for subsets within the same mixed model?

1 Answers1