An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20110624/0a407033/attachment.pl>
code for multiple membership models?
5 messages · Douglas Bates, Malcolm Fairbrother, Jarrod Hadfield
Hi,
MCMCglmm solution below, I think. It's a bit ugly:
library(foreign); lips <-
read.dta("http://www.bristol.ac.uk/cmm/media/runmlwin/lips1.dta")[,c(1,3,5,9:30)]
neigh.factors<-unique(unlist(lips[,4:14]))
neigh.factors<-neigh.factors[-which(neigh.factors==0)]
# vector of all potential neighbours (0 really means missing)
for(i in 4:14){
lips[,i]<-factor(lips[,i], levels=neigh.factors)
# make sure neigh1, neigh2 ... all have the same potential levels
lips[,i][which(is.na(lips[,i]))]<-neigh.factors[1]
# change missing values to arbitrary neighbour (makes no difference
since weight is zero). I should probably allow missing values!
}
m1<-MCMCglmm(obs~1,
random=~idv(mult.memb(~weight1:neigh1+weight2:neigh2+weight3:neigh3+weight4:neigh4+weight5:neigh5+weight6:neigh6+weight7:neigh7+weight8:neigh8+ weight9:neigh9+weight10:neigh10+weight11:neigh11)),
data=lips)
# It might make more sense to square root the weights depending on
model assumptions.
Cheers,
Jarrod
Quoting Malcolm Fairbrother <m.fairbrother at bristol.ac.uk> on Fri, 24
Jun 2011 13:18:05 +0100:
Dear Doug, Jarrod, and/or perhaps others, A couple posters have previously asked about multiple membership models, and Doug has said he could provide code to fit such models (https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q2/003860.html), if given some sample data. Jarrod has similarly suggested that MCMCglmm can handle multiple membership (https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q2/006195.html), though I can't find any reference to multiple membership models in the MCMCglmm documentation. (Maybe I missed it, however.) Here's is a simple sample dataset (Scottish Lip Cancer): library(foreign); lips <- read.dta("http://www.bristol.ac.uk/cmm/media/runmlwin/lips1.dta")[,c(1,3,5,9:30)] I would be very grateful for any indication of how to fit a multiple membership model to such data. In this instance, the outcome ("obs") is Poisson-distributed (number of instances of lip cancer in an area), and "perc_aff" is a covariate (percentage of the population working in outdoor activities). Each unit ("area") is observed only once, but a mixed model could be useful for modelling the spatial dependence in the data, where each area is taken to be nested in a set of neighbours (identified in columns 4-14). Each area has anywhere from 1 to 11 neighbours, so each neighbour's share of the overall "neighbour" classification would have to be weighted somewhere from 1/11 to 1. Weights are given in columns 14-25. So a starting model is: summary(M1 <- glm(obs ~ 1, lips, family=poisson)) Then: library(lme4) (M2 <- lmer(obs ~ 1 + (1 | neigh1), lips, family=poisson)) However, the random effect in M2 should be replaced by a random effect for "neighbours" (not just "neigh1") where each of the 1 to 11 neighbours "counts equally" towards the classification, not just the area's first neighbour. Can this indeed be done in lme4 and/or MCMCglmm? If so, can you please show how? A similar analysis of this dataset, using MLwiN, is described in chapter 17 of http://www.bristol.ac.uk/cmm/software/mlwin/download/mcmc-09.pdf. Thanks, Malcolm [[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
3 days later
My apologies for the delay on responding to this question and to a similar question off-list by Jeroen Ooms. I have been immersed in another development effort which ultimately will feed back into lme4 and, unfortunately, I don't multiplex well. The basic approach to multiple membership models is to build the model from one of the factors representing the membership, such as area here, but not fit it - just return the numeric representation. Specifying doFit=FALSE in a call to glmer in the lme4a package provides the raw structure without doing the optimization. Inside that structure is the sparse matrix Z. You then add the sparse matrices from the other factors representing the membership, recalculate the Cholesky factor L (because the structure of Z has changed) and go on with the optimization. I enclose a first cut at this. I didn't use weights so the results are not directly comparable to those in the manual from the Multilevel Modeling group. On Fri, Jun 24, 2011 at 7:18 AM, Malcolm Fairbrother
<m.fairbrother at bristol.ac.uk> wrote:
Dear Doug, Jarrod, and/or perhaps others, A couple posters have previously asked about multiple membership models, and Doug has said he could provide code to fit such models (https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q2/003860.html), if given some sample data. Jarrod has similarly suggested that MCMCglmm can handle multiple membership (https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q2/006195.html), though I can't find any reference to multiple membership models in the MCMCglmm documentation. (Maybe I missed it, however.) Here's is a simple sample dataset (Scottish Lip Cancer): library(foreign); lips <- read.dta("http://www.bristol.ac.uk/cmm/media/runmlwin/lips1.dta")[,c(1,3,5,9:30)] I would be very grateful for any indication of how to fit a multiple membership model to such data. In this instance, the outcome ("obs") is Poisson-distributed (number of instances of lip cancer in an area), and "perc_aff" is a covariate (percentage of the population working in outdoor activities). Each unit ("area") is observed only once, but a mixed model could be useful for modelling the spatial dependence in the data, where each area is taken to be nested in a set of neighbours (identified in columns 4-14). Each area has anywhere from 1 to 11 neighbours, so each neighbour's share of the overall "neighbour" classification would have to be weighted somewhere from 1/11 to 1. Weights are given in columns 14-25. So a starting model is: summary(M1 <- glm(obs ~ 1, lips, family=poisson)) Then: library(lme4) (M2 <- lmer(obs ~ 1 + (1 | neigh1), lips, family=poisson)) However, the random effect in M2 should be replaced by a random effect for "neighbours" (not just "neigh1") where each of the 1 to 11 neighbours "counts equally" towards the classification, not just the area's first neighbour. Can this indeed be done in lme4 and/or MCMCglmm? If so, can you please show how? A similar analysis of this dataset, using MLwiN, is described in chapter 17 of http://www.bristol.ac.uk/cmm/software/mlwin/download/mcmc-09.pdf. Thanks, Malcolm ? ? ? ?[[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-------------- next part -------------- R version 2.14.0 Under development (unstable) (2011-06-24 r56210) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: x86_64-unknown-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R.
library(lme4a)
Loading required package: Matrix
Loading required package: lattice
Attaching package: ?Matrix?
The following object(s) are masked from ?package:base?:
det
Loading required package: minqa
Loading required package: Rcpp
Loading required package: MatrixModels
Attaching package: ?MatrixModels?
The following object(s) are masked from ?package:stats?:
getCall
library(foreign)
## convert the neighXX columns to factors with a consistent set of levels
## drop the weights columns
lips <- within(read.dta("http://www.bristol.ac.uk/cmm/media/runmlwin/lips1.dta")[,c(1,3,5,9:30)],
+ {
+ area <- factor(area)
+ neigh1 <- factor(neigh1, levels=1:56)
+ neigh2 <- factor(neigh2, levels=1:56)
+ neigh3 <- factor(neigh3, levels=1:56)
+ neigh4 <- factor(neigh4, levels=1:56)
+ neigh5 <- factor(neigh5, levels=1:56)
+ neigh6 <- factor(neigh6, levels=1:56)
+ neigh7 <- factor(neigh7, levels=1:56)
+ neigh8 <- factor(neigh8, levels=1:56)
+ neigh9 <- factor(neigh9, levels=1:56)
+ neigh10 <- factor(neigh10, levels=1:56)
+ neigh11 <- factor(neigh11, levels=1:56)
+ })[ , 1:14]
## truncate the names
names(lips)[4:14] <- paste(sprintf("n%02d", 1:11))
str(lips)
'data.frame': 56 obs. of 14 variables: $ area : Factor w/ 56 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ... $ obs : int 9 39 11 9 15 8 26 7 6 20 ... $ perc_aff: int 16 16 10 24 10 24 10 7 7 16 ... $ n01 : Factor w/ 56 levels "1","2","3","4",..: 5 7 6 18 1 3 2 6 1 2 ... $ n02 : Factor w/ 56 levels "1","2","3","4",..: 9 10 12 20 11 8 10 NA 11 7 ... $ n03 : Factor w/ 56 levels "1","2","3","4",..: 11 NA NA 28 12 NA 13 NA 17 16 ... $ n04 : Factor w/ 56 levels "1","2","3","4",..: 19 NA NA NA 13 NA 16 NA 19 22 ... $ n05 : Factor w/ 56 levels "1","2","3","4",..: NA NA NA NA 19 NA 17 NA 23 NA ... $ n06 : Factor w/ 56 levels "1","2","3","4",..: NA NA NA NA NA NA NA NA 29 NA ... $ n07 : Factor w/ 56 levels "1","2","3","4",..: NA NA NA NA NA NA NA NA NA NA ... $ n08 : Factor w/ 56 levels "1","2","3","4",..: NA NA NA NA NA NA NA NA NA NA ... $ n09 : Factor w/ 56 levels "1","2","3","4",..: NA NA NA NA NA NA NA NA NA NA ... $ n10 : Factor w/ 56 levels "1","2","3","4",..: NA NA NA NA NA NA NA NA NA NA ... $ n11 : Factor w/ 56 levels "1","2","3","4",..: NA NA NA NA NA NA NA NA NA NA ...
## Convert each nXX column to a sparse model matrix. We can almost do ## this with as(lips$n01, "sparseMatrix") but that drops unused levels ## and we must keep them. Instead we use the hidden function behind ## that coercion function. Zl <- lapply(names(lips[4:14]), # list of sparse model matrices, one for each column
+ function(nm) Matrix:::fac2sparse(lips[[nm]], "d", drop=FALSE))
ZZ <- Reduce("+", Zl[-1], Zl[[1]])
str(ZZ)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:264] 4 8 10 18 6 9 5 11 17 19 ... ..@ p : int [1:57] 0 4 6 8 11 16 18 23 24 30 ... ..@ Dim : int [1:2] 56 56 ..@ Dimnames:List of 2 .. ..$ : chr [1:56] "1" "2" "3" "4" ... .. ..$ : NULL ..@ x : num [1:264] 1 1 1 1 1 1 1 1 1 1 ... ..@ factors : list()
## generate, but do not fit, the model structure mstr <- glmer(obs ~ perc_aff + (1|area), lips, poisson, doFit=FALSE) str(mstr at re)
Formal class 'reTrms' [package "lme4a"] with 9 slots ..@ flist :List of 1 .. ..$ area: Factor w/ 56 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ... .. ..- attr(*, "assign")= int 1 ..@ cnms :List of 1 .. ..$ area: chr "(Intercept)" ..@ L :Formal class 'dCHMsimpl' [package "Matrix"] with 10 slots .. .. ..@ x : num [1:56] 1.41 1.41 1.41 1.41 1.41 ... .. .. ..@ p : int [1:57] 0 1 2 3 4 5 6 7 8 9 ... .. .. ..@ i : int [1:56] 0 1 2 3 4 5 6 7 8 9 ... .. .. ..@ nz : int [1:56] 1 1 1 1 1 1 1 1 1 1 ... .. .. ..@ nxt : int [1:58] 1 2 3 4 5 6 7 8 9 10 ... .. .. ..@ prv : int [1:58] 57 0 1 2 3 4 5 6 7 8 ... .. .. ..@ colcount: int [1:56] 1 1 1 1 1 1 1 1 1 1 ... .. .. ..@ perm : int [1:56] 0 1 2 3 4 5 6 7 8 9 ... .. .. ..@ type : int [1:4] 2 1 0 1 .. .. ..@ Dim : int [1:2] 56 56 ..@ Lambda:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots .. .. ..@ i : int [1:56] 0 1 2 3 4 5 6 7 8 9 ... .. .. ..@ p : int [1:57] 0 1 2 3 4 5 6 7 8 9 ... .. .. ..@ Dim : int [1:2] 56 56 .. .. ..@ Dimnames:List of 2 .. .. .. ..$ : NULL .. .. .. ..$ : NULL .. .. ..@ x : num [1:56] 1 1 1 1 1 1 1 1 1 1 ... .. .. ..@ factors : list() ..@ Lind : int [1:56] 1 1 1 1 1 1 1 1 1 1 ... ..@ Zt :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots .. .. ..@ i : int [1:56] 0 1 2 3 4 5 6 7 8 9 ... .. .. ..@ p : int [1:57] 0 1 2 3 4 5 6 7 8 9 ... .. .. ..@ Dim : int [1:2] 56 56 .. .. ..@ Dimnames:List of 2 .. .. .. ..$ : chr [1:56] "1" "2" "3" "4" ... .. .. .. ..$ : NULL .. .. ..@ x : num [1:56] 1 1 1 1 1 1 1 1 1 1 ... .. .. ..@ factors : list() ..@ lower : num 0 ..@ theta : num 1 ..@ u : num [1:56] 0 0 0 0 0 0 0 0 0 0 ...
# replace the Zt matrix and L factor
re <- mstr at re
re at Zt <- ZZ
re at L <- Cholesky(tcrossprod(ZZ), LDL=FALSE, Imult=1)
mstr at re <- re
# fit and summarize the model
summary(ans <- lme4a:::PIRLSest(mstr, verbose=2L, control=list(), nAGQ=1L))
npt = 3 , n = 1
rhobeg = 0.2 , rhoend = 2e-07
0.020: 5: 164.386;0.600000
0.0020: 12: 159.497;0.420054
0.00020: 17: 159.485;0.427453
2.0e-05: 19: 159.485;0.427058
2.0e-06: 20: 159.485;0.427026
2.0e-07: 22: 159.485;0.427026
At return
25: 159.48455: 0.427025
npt = 5 , n = 3
rhobeg = 0.2 , rhoend = 2e-07
0.020: 6: 159.484;0.427025 2.45348 -0.00103278
0.0020: 9: 159.484;0.427025 2.45348 -0.00103278
0.00020: 17: 159.462;0.428080 2.43929 -0.00209994
2.0e-05: 51: 159.456;0.426677 2.42444 -0.00112120
2.0e-06: 77: 159.456;0.426957 2.42220 -0.00100983
2.0e-07: 94: 159.456;0.426955 2.42214 -0.00100766
At return
110: 159.45636: 0.426955 2.42214 -0.00100739
Generalized linear mixed model fit by maximum likelihood ['summary.mer']
Family: poisson
Formula: obs ~ perc_aff + (1 | area)
Data: lips
AIC BIC logLik deviance
165.4564 171.5324 -79.7282 159.4564
Random effects:
Groups Name Variance Std.Dev.
area (Intercept) 0.1823 0.427
Number of obs: 56, groups: area, 56
Fixed effects:
Estimate Std. Error z value
(Intercept) 2.422140 0.244911 9.890
perc_aff -0.001007 0.015505 -0.065
Correlation of Fixed Effects:
(Intr)
perc_aff -0.646
proc.time()
user system elapsed 10.730 0.110 11.611
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)
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.