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)