MCMCglmm multidimensional model in R

Question

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

The data is structured as follows: where dyad, focal, others are all random effects, prediction 1-2 are predictor variables, and answer 1-5 are result variables that record # 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 one result (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) 

The Hadfield Course Notes (Ch 5) provide an example of a multinomial model that uses only one result variable with three levels (sheep’s horns of 3 types). A similar reference can be found here: http://hlplab.wordpress.com/2009/05/07/multinomial-random-effects-models-in-r/ This is not entirely correct for what I am doing, but it contains useful background information .

Another link (Hadfield 2010) provides an example of MCMCglmm with multiple answers, which follows the same format but uses cbind () to predict the response vector, and not for a single result. The same model with multiple answers will 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:

• How to specify the previous model? I looked at the materials mentioned in this post, but just can't figure it out.

• I run a similar version with only two response variables, but I get only one slope - where I thought there should be a different slope for each resp variable. Where am I mistaken, or did I misunderstand the model?

Answer 1

Answer my first question based on an HLP post and some help from a peer / specialist 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 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 

And for my model, using the multinomial5 variable and 5 source variables, the previous one:

 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 the main effects for the Response variables and the rear estimates of the Response / Predictor interaction (the effect of the predictor variable on each response variable).