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

[Jarrod Hadfield]

Hi,

The chain is getting "stuck" at values of zero for the variance, hence  
the spike. This should happen less with proper priors, and parameter  
expansion in particular increases mixing under these conditions. Using:

prior=list(R=list(V=1, nu=0.002), G=list(G1=list(V=1, nu=1,  
alpha.mu=0, alpha.V=1000),G2=list(V=1, nu=1, alpha.mu=0, alpha.V=1000)))

as your prior should help.

Note that the variance in the data due to the second set of random  
effects is not equal to the variance component because the diagonal  
elements of ZZ' are on average ~0.07. Depending on model assumptions  
you might want to square root your weights.  The diagonal elements of  
ZZ' are then all one and the variance component can be interpreted  
directly.  I can't see how the weights enter into the lme4a code, but  
presumably they do?

Jarrod



Quoting Malcolm Fairbrother <m.fairbrother at bristol.ac.uk> on Tue, 28  
Jun 2011 17:04:03 +0100:

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



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.