LMER/GLMER How to generate random effects design matrix?

Question

Just to illustrate my situation, imagine this is my data:

y = rbinom(25, 1, 0.5)
names1 = rep(1:5, each = 5)
names2 = rep(1:5 , 5)
x= rnorm(25)
df = data.frame(cbind(names1, names2, y, x))

That is, it looks like the following:

names1	names2	y	x
1	2	0	-1.10152733
1	3	1	1.49929490
1	4	0	0.96740978
1	5	1	-0.52383757
2	1	0	-1.19182503
2	3	1	1.00492382
2	4	1	0.39198228
2	5	1	-0.68029102

I want to estimate y here using glmer, ie:

glmer(y ~ x + (1|Random), family = binomial, data=df)

For my problem, Random effects structure should be based on names, independent of them appearing on names1 or names 2 column. That is, in the general formulation of the model $y=X\beta + Zb +\u$, the matrix Z here is n*p, where in this exercise n=20 (nrow(df)), and p=5 (length(unique(df$names1)). $Z_{ij}=1$ if the name $j$ appear in row $i$, that is names1=j OR names2=j

My question is how can I use this matrix Z as my random effects grouping:

glmer(y ~ x + (1|Z), family = binomial, data=df)

So far, I have generated Z as:

for (i in 1:length(unique(df$names1))) {
     df[paste0("Z", i)] <- ifelse(df$names1 == i | df$names2 == i, 1, 0)
 }
df = df[df$names1 != df$names2,]

But i am at lost in how to incorporate Z's in glmer. Any help would be appreciated.

Answer 1

I don't have enough reputation to comment or edit, so hopefully someone can merge this into Ben's answer:

lmerMultiMember is an R package written specifically to handle this type of multiple membership model. Using it is pretty straightforward, the only reason Ben's example doesn't work is that unlike lme4, lmerMultiMember requires you pass the data in through the data argument of the lmer() or glmer() function call. (This is a technical limitation, unfortunately.)

That means you can't just use variables the model formula when they exist in the top-level R environment, they have to be in the dataframe you pass in through the data argument.

The example below works just fine, and does the same thing as Ben's example:

set.seed(101)
y <- rbinom(25, 1, 0.5)
names1 <- rep(1:5, each = 5)
names2 <- rep(1:5, 5)
x <- rnorm(25)
df <- data.frame(names1, names2, y, x) ## better not to cbind() unnecessarily

## if necessary: remotes::install_github("jvparidon/lmerMultiMember")
library(lmerMultiMember)

## create weight matrix for multiple membership model
W <- weights_from_vector(paste(df$names1, df$names2, sep = ","))

## fit model and inspect summary
m <- lmerMultiMember::lmer(y ~ x + (1 | names), memberships = list(names = W), data = df)
print(summary(m))

The vignettes at https://jvparidon.github.io/lmerMultiMember/ have some more usage examples that might be of help.

Answer 2

The name for this setup is a multi-membership model.

This can be hacked in lme4 following the recipe here, which I illustrate below. Alternatively, there is an lmerMultiMember package on Github which is supposed to fit these models (but I couldn't get it to work, so I show how to do it the hacky way).

set up data

set.seed(101)
y = rbinom(25, 1, 0.5)
names1 = rep(1:5, each = 5)
names2 = rep(1:5 , 5)
x= rnorm(25)
df = data.frame(names1, names2, y, x) ## better not to cbind() unnecessarily

## if necessary: remotes::install_github("jvparidon/lmerMultiMember")
library(lmerMultiMember)
M <- weights_from_vector(paste(df$names1, df$names2, sep = ","))
try(lmerMultiMember::lmer(y ~ x + (1 | names),
                      memberships = list(names = M)))

This throws an error for reasons I don't understand. I'll use M (which is the transpose of the Z matrix) below instead:

library(lme4)
## a fake factor with the right levels, actual values don't matter
fake <- rep(unique(c(df$names1, df$names2)),length.out=nrow(df))
## set up model structures
lmod <- lFormula(y~x+(1|fake), data = df)
## substitute the membership matrix for Z-transpose
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- M
## set up and fit the model
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m1 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)

You can also fit multi-membership models in the MCMCglmm package, if you want to go Bayesian ...