MCMCglmm multinomial model in R

Question

I'm trying to create a model using the MCMCglmm package in R.

The data are structured as follows, where dyad, focal, other are all random effects, predict1-2 are predictor variables, and response 1-5 are outcome variables that capture # of observed behaviors of different subtypes:

 dyad focal other r    present  village  resp1 resp2 resp3 resp4 resp5 
 1    10101 14302 0.5  3        1        0     0     4     0     5
 2    10405 11301 0.0  5        0        0     0     1     0     1
 …

So a model with only one outcome (teaching) is as follows:

 prior_overdisp_i <- list(R=list(V=diag(2),nu=0.08,fix=2), 
 G=list(G1=list(V=1,nu=0.08), G2=list(V=1,nu=0.08), G3=list(V=1,nu=0.08), G4=list(V=1,nu=0.08)))

 m1 <- MCMCglmm(teaching ~ trait-1 + at.level(trait,1):r + at.level(trait,1):present, 
 random= ~idh(at.level(trait,1)):focal + idh(at.level(trait,1)):other + 
 idh(at.level(trait,1)):X + idh(at.level(trait,1)):village, 
 rcov=~idh(trait):units, family = "zipoisson", prior=prior_overdisp_i, 
 data = data, nitt = nitt.1, thin = 50, burnin = 15000, pr = TRUE, pl = TRUE, verbose = TRUE, DIC  = TRUE)

Hadfield's course notes (Ch 5) give an example of a multinomial model that uses only a single outcome variable with 3 levels (sheep horns of 3 types). Similar treatment can be found here: http://hlplab.wordpress.com/2009/05/07/multinomial-random-effects-models-in-r/ This is not quite right for what I'm doing, but contains helpful background info.

Another reference (Hadfield 2010) gives an example of a multi-response MCMCglmm that follows the same format but uses cbind() to predict a vector of responses, rather than a single outcome. The same model with multiple responses would look like this:

 m1 <- MCMCglmm(cbind(resp1, resp2, resp3, resp4, resp5) ~ trait-1 + 
 at.level(trait,1):r + at.level(trait,1):present, 
 random= ~idh(at.level(trait,1)):focal + idh(at.level(trait,1)):other + 
 idh(at.level(trait,1)):X + idh(at.level(trait,1)):village, 
 rcov=~idh(trait):units, 
 family = cbind("zipoisson","zipoisson","zipoisson","zipoisson","zipoisson"), 
 prior=prior_overdisp_i, 
 data = data, nitt = nitt.1, thin = 50, burnin = 15000, pr = TRUE, pl = TRUE, verbose = TRUE, DIC  = TRUE)

I have two programming questions here:

1. How do I specify a prior for this model? I've looked at the materials mentioned in this post but just can't figure it out.

2. I've run a similar version with only two response variables, but I only get one slope - where I thought I should get a different slope for each resp variable. Where am I going wrong, or having I misunderstood the model?

Mosley

Did you check, whether fix = 2 in R=list(V=diag(2),nu=0.08,fix=2) really makes sense? In my understanding of the MCMCglmm's prior specification fix should be read like a boolean value: fix = 0 is the default value for not fixing the variance to V, and fix = 1 means "fix the variance at the value of V". So fix = 2 (or similar) imo should have no meaning at all. (But on page 103 of his course nots Hadfield use this specification: cran.r-project.org/pub/R/web/packages/MCMCglmm/vignettes/…) – Steersman

@Steersman I am coming back to these data after a couple of years and looking at these models again. My understanding is that the "fix" component has to do with which part of the model the prior is for... since there is a categorical component (predicting zeros) and a continuous component (predicting non-zero values). This is specific to zipoisson models which are technically multinomial in their own right. Caveat: I may be confused! – Mosley

Answer 1

Answer to my first question, based on the HLP post and some help from a colleage/stats consultant:

# values for prior
k <- 5 # originally: length(levels(dative$SemanticClass)), so k = # of outcomes for SemanticClass     aka categorical outcomes 
I <- diag(k-1) #should make matrix of 0's with diagonal of 1's, dimensions k-1 rows and k-1 columns
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1)) # should make k-1 x k-1 matrix of 1's 

And for my model, using the multinomial5 family and 5 outcome variables, the prior is:

prior = list(
             R = list(fix=1, V=0.5 * (I + J), n = 4),
             G = list(
               G1 = list(V = diag(4), n = 4))

For my second question, I need to add an interaction term to the fixed effects in this model:

 m <- MCMCglmm(cbind(Resp1, Resp2...) ~ -1 + trait*predictorvariable,
 ...

The result gives both main effects for the Response variables and posterior estimates for the Response/Predictor interaction (the effect of the predictor variable on each response variable).

Mosley