[R-sig-ME] MCMCglmm for more than two levels

[David Atkins]

 >Dear All,
 >
 >I am trying to estimating a three-level multinomial logit model using 
 >MCMCglmm (I am not aware of other options for random effect
 >multinomial logit models). My question is: can MCMCglmm handle more
 >than two levels of nesting? I have been trying to find answers from
 >the documents without luck, and would really appreciate for some
 > examples if it is feasible at all. Many thanks.
 >
 >Best,
 >Shige

Yes, MCMCglmm can fit multiple levels, though outside of social 
sciences, I don't think mixed-model folks think in terms of "levels" 
(quite possibly with good reason, due to cross-classification, etc.).

I've included an example below that uses data from a marital therapy 
study (which I used in the following article):

Atkins, D. C. (2005). Using multilevel models to analyze marital and 
family treatment data: Basic and advanced issues. Journal of Family 
Psychology, 19, 98-110.

The data has 4 repeated measures, on 268 individuals, in 134 couples; 
thus, a 3-level model in Raudenbush and Bryk jargon.

The outcome is the "dyadic adjustment scale", which is reasonably 
normally distributed; however, I think the code for MCMCglmm should be 
similar, with the exception of an additional "family" argument.

Hope it helps.

cheers, Dave



df <- read.table(file = 
"http://depts.washington.edu/cshrb/newweb/stats%20documents/Atkins%202005%20data.txt",
					header = TRUE)
head(df)					
names(df) <- casefold(names(df))			

### create unique individual ID
df$ind.id <- with(df, interaction(id, sex))

### fit "3 level" model with lmer
library(lme4)		

out.lmer <- lmer(das ~ therapy*time + (1 | id) + (1 | ind.id),
					data = df, verbose = TRUE)
summary(out.lmer)					

### fit random slopes across couples
out.lmer2 <- lmer(das ~ therapy*time + (time | id) + (1 | ind.id),
					data = df, verbose = TRUE)
summary(out.lmer2)					

### fit same models using MCMCglmm
library(MCMCglmm)

out.mcmc <- MCMCglmm(das ~ therapy*time,
					random = ~ id + ind.id,
					data = df, verbose = TRUE)
summary(out.mcmc$Sol)
summary(out.mcmc$VCV)
### results similar to out.lmer
### NOTE: this is functionally random-intercepts

prior <- list(R = list(V = 1, nu = 1),
				G = list(G1 = list(V = diag(2), nu = 2),
						  G2 = list(V = 1, nu = 1)))

out.mcmc2 <- MCMCglmm(das ~ therapy*time,
					random = ~ us(1 + time):id + ind.id,
					rcov = ~ units,
					data = df, verbose = TRUE,
					prior = prior)
summary(out.mcmc2$Sol)
summary(out.mcmc2$VCV)
### NOTE: this fits random intercept and slope for couples and
### intercept only for individuals
### results largely comparable to out.lmer2;
### using iterations, burn-in, thinning "out of the box"

-- 
Dave Atkins, PhD
Research Associate Professor
Department of Psychiatry and Behavioral Science
University of Washington
datkins at u.washington.edu

Center for the Study of Health and Risk Behaviors (CSHRB)		
1100 NE 45th Street, Suite 300 	
Seattle, WA  98105 	
206-616-3879 	
http://depts.washington.edu/cshrb/
(Mon-Wed)	

Center for Healthcare Improvement, for Addictions, Mental Illness,
   Medically Vulnerable Populations (CHAMMP)
325 9th Avenue, 2HH-15
Box 359911
Seattle, WA 98104?
206-897-4210
http://www.chammp.org
(Thurs)