6
$\begingroup$

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 us="unstructured" and 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.

  1. Is it necessary to use trait in 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)?
  2. Does the inclusion of trait change the estimated formula or is it only a signifier for MCMCglmm to expect reshaped data?
  3. 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. 

Reference:
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.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.