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

[Shige Song]

Dear Dave,

This is very helpful, thank you very much!

Best,
Shige

On Thu, Apr 22, 2010 at 12:07 AM, David Atkins <datkins at u.washington.edu> wrote:
>>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)
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>