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.