[R-sig-ME] code for multiple membership models?

[Malcolm Fairbrother]

Dear Jarrod and Doug,

Thanks very much to both of you for the code. I think I understand more or less what is going on, in both instances, and I would guess that this will be useful content for other people via the archives.

Pushing this further, I tried to simulate some data--slightly more complex--and use both MCMCglmm and lme4a to try to recover the parameters (see code below). The results are encouraging (the two functions return similar results), but not entirely so, and was hoping I might abuse your patience just a bit more on this topic.

What I've done here is to simulate grouped data, where each level-1 unit is uniquely nested in a higher-level unit AND is treated as a member of a set of neighbouring higher-level units--where the number of neighbouring units can vary. So in lmer terms, the model looks like y ~ X1 + X2 + (1 | L1) + (1 | L2), where the X's are covariates, L1 is a garden variety grouping factor, and L2 is the set of neighbours. The complexity (need for a multiple membership approach) comes from L2. The outcome here is Normally distributed, unlike my earlier example.

(A) I think I was able to adapt Doug's lme4a code to this situation (see middle section of the code below), but I'm not 100% sure for a couple reasons:
(1) When I specify but do not fit the model, it returns an error/warning. It doesn't sound too bad, but I don't know for sure: "Error in if (length(e1) <= 1 && any(.Generic == c("*", "/", "+")) && (e1 >  :   missing value where TRUE/FALSE needed".
(2) The switch from a Poisson to Normally distributed outcome seems to require new code in ways that dig into the guts of lme4a (i.e., different optimizer). Does the code below get this right?
(3) The results from an otherwise identical model missing the L2 random term are--except for the absence of an estimate of the Variance/SD for L2--otherwise identical. In other words, the beta coefficients, and the Variance/SD estimates for L1 and the Residuals, are identical for the models with and without L2. That seems suspicious to me.
(4) I'm not sure I'm simulating data in a way that makes perfect sense, but obviously that's my issue not yours and I'll not worry about that for now.

(B) For MCMCglmm (see bottom section of the code below), the code seems straightforward, but I'm finding that the chain for the variance of L2 tends to experience occasional massive spikes, such that the difference between the mean and median for that chain is substantial, the density plot doesn't look good, and the confidence interval is very wide. This problem also occurs for some real data to which I've tried applying the method. The median is reasonably close to the result from lme4a, but the huge spikes in the chain don't leave me very confident.

If this makes sense, what do you think about these results? Are they cause for concern? Can they be addressed, such as by resolving problems I'm overlooking in the code below? Many thanks once again for your help.

- Malcolm





# SIMULATE DATA
set.seed(1234)
N <- 30 # set the number of lower-level units per higher level unit
N2 <- 60 # set the number of higher level units
dat <- data.frame(X1=rnorm(N*N2, 2, 2), X2=rep(rnorm(N2, 2, 2), each=N), L1=rep(seq(from=1, to=N2), each=N)) # create data
L1 <- matrix(0, nrow=nrow(dat), ncol=N2) # this line and the next create an N*N2 by N2 matrix of group membership
for (i in 1:nrow(L1)) { L1[i,(dat$L1[i])] <- 1 }
dat$U_j1 <- L1 %*% rnorm(N2, 0, 4) # group-level random intercept
L2 <- matrix(rbinom(N2^2, 1, 0.25), N2, N2)*upper.tri(matrix(1, N2, N2)) # random connectivity matrix, sized N2 by N2
L2 <- L2 + t(L2) # made symmetrical
L2 <- L1 %*% L2 # expanded to N*N2 by N2
L2 <- L2/apply(L2, 1, sum) # row-standardised (each row sums to 1)
dat$U_j2 <- L2 %*% rnorm(N2, 0, 4) # group-level random intercept for neighbours
dat$y <- dat$X1 + dat$X2 + dat$U_j1 + dat$U_j2 + rnorm(nrow(dat), 1, 4) # y data including Intercept (=1) and random error
L2 <- apply(L2, 1, function(x) which(x>0))
L22 <- data.frame(matrix(NA, nrow=length(L2), ncol=max(unlist(lapply(L2, length))))) # new kind of connectivity matrix
for (i in 1:nrow(L22)) L22[i,1:length(L2[[i]])] <- L2[[i]]
names(L22) <- paste(sprintf("n%02d", 1:ncol(L22)))
w <- data.frame(t(apply(L22, 1, function(x) as.numeric(!is.na(x))/sum(!is.na(x))))) # weights matrix
for (i in 1:ncol(L22)) L22[,i] <- factor(L22[,i], levels=unique(dat$L1))
names(w) <- paste(sprintf("w%02d", 1:ncol(w)))
dat <- data.frame(dat, L22, w)

# ESTIMATION WITH lme4a
library(lme4a)
Zl <- lapply(names(dat)[7:(6+ncol(L22))], function(nm) Matrix:::fac2sparse(dat[[nm]], "d", drop=FALSE)) # list of sparse model matrices
ZZ <- Reduce("+", Zl[-1], Zl[[1]])
dat$L2 <- dat$L1 # an "artificial" grouping factor, used only in setting up the model (in the next line)
mstr <- lmer(y ~ X1 + X2 + (1 | L1) + (1 | L2), dat, doFit=FALSE) # throws an error, but doesn't seem to matter
re <- mstr at re
re at Zt[(1+nrow(re at Zt)/2):nrow(re at Zt),] <- ZZ # replace the Zt
re at L <- Cholesky(tcrossprod(re at Zt), LDL=FALSE, Imult=1) # replace the L factor
mstr at re <- re
opt <- bobyqa(mstr at re@theta, mkdevfun(mstr), mstr at re@lower, control = list()) # optimize
summary(mstr <- lme4a:::updateMod(mstr, opt$par, opt$fval)) # summarize the model
(mstr2 <- lmer(y ~ X1 + X2 + (1 | L1), dat)) # for comparison

# ESTIMATION WITH MCMCglmm
library(MCMCglmm)
for(i in 7:(6+ncol(L22))) { dat[,i][which(is.na(dat[,i]))] <- levels(dat$n01)[1] }
m1 <- MCMCglmm(y ~ X1 + X2, random = ~ L1 + idv(mult.memb(~ w01:n01 + w02:n02 + w03:n03 + w04:n04 + w05:n05 + w06:n06 + w07:n07 + w08:n08 + w09:n09 + w10:n10 + w11:n11 + w12:n12 + w13:n13 + w14:n14 + w15:n15 + w16:n16 + w17:n17 + w18:n18 + w19:n19 + w20:n20 + w21:n21)), data=dat)
rbind(as.numeric(c(apply(m1$VCV, 2, mean), apply(m1$Sol, 2, mean))), as.numeric(c(apply(m1$VCV, 2, median), apply(m1$Sol, 2, median))))
plot(m1)