I want to estimate a multilevel multinomial logit model but I am struggling with the terminology and notation used by the R-package
MCMCglmm. There is documentation available in form of a J Stat Softw article and course notes.
According to the article, "the [MCMCglmm] syntax used to specify the model closely follows that used by asreml" (Hadfield 2010: 6). Consequently, some abbreviations like
idh="independent heterogeneous" are never stated explicitly in the MCMCglmm documentation, although they are well presented in formulas. Unfortunately, I am only acquainted with the notation for hierarchical models in the social sciences and I have no experience with ASReml.
The examples of multinomial multilevel models in the MCMCglmm course notes are subsumed under multiresponse models and the course notes (p. 91) show that
trait is used as identifier when the dataset is reshaped for estimating multiresponse models.
- Is it necessary to use
traitin the fixed (or random) formula also for models with a nominal outcome variable like employment status or gender (in contrast to a real multiresponse variable)?
- Does the inclusion of
traitchange the estimated formula or is it only a signifier for MCMCglmm to expect reshaped data?
- Can I define the error structure without using
trait, e.g. by specifying the error structure for coefficients from the fixed part?
Here is an (almost working) example code using a typical social science survey dataset:
library(nnet) library(MCMCglmm) # Load dataset from TraMineR-package data(mvad, package="TraMineR") ### Create variable "region" for 2nd level attach(mvad) mvad$region[Belfast=="yes"] <- "Belfast" mvad$region[N.Eastern=="yes"] <- "N.Eastern" mvad$region[Southern=="yes"] <- "Southern" mvad$region[S.Eastern=="yes"] <- "S.Eastern" mvad$region[Western=="yes"] <- "Western" detach(mvad) # Add identifier for complete cases mvad$complete <- complete.cases(mvad[, c("Jul.93", "male", "Grammar", "gcse5eq", "funemp", "region")]) # Multinomial logit model without random effects using nnet-package multinom(Jul.93 ~ male + Grammar + gcse5eq + funemp + region, data = mvad[mvad$complete==TRUE, ]) # Multinomial logit models with random effects using MCMCglmm mod1 <- MCMCglmm(Jul.93 ~ trait + male + Grammar + gcse5eq + funemp, random = ~ region, rcov = ~ us(trait):units, family = "categorical", verbose=FALSE, data=mvad[mvad$complete==TRUE, ]) # Including trait and defining rcov gives an estimable command mod2 <- MCMCglmm(Jul.93 ~ male + Grammar + gcse5eq + funemp, random = ~ region, rcov = ~ us(trait):units, family = "categorical", verbose=FALSE, data=mvad[mvad$complete==TRUE, ]) # Command works as well, but is trait now only used for the error terms? mod3 <- MCMCglmm(Jul.93 ~ male + Grammar + gcse5eq + funemp, random = ~ region, family = "categorical", verbose=FALSE, data=mvad[mvad$complete==TRUE, ]) # Command does not work as definition of error structure using e.g. "us(trait):units" is needed.
Hadfield, J.D. (2010) 'MCMC Methods for Multi-Response Generalized Linear Mixed Models: The MCMCglmm R Package', Journal of Statistical Software 33(2): 1–22.