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

Jarrod Hadfield j.hadfield at ed.ac.uk
Fri Jun 24 15:20:29 CEST 2011


MCMCglmm solution below, I think. It's a bit ugly:

library(foreign); lips <-  


# 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
# change missing values to arbitrary neighbour (makes no difference  
since weight is zero). I should probably allow missing values!

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

# It might make more sense to square root the weights depending on  
model assumptions.



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.

More information about the R-sig-mixed-models mailing list