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

[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.