I am using MCMCglmm models in R, with hierarchically nested data. The basic structure of the data is as follows - each dyad is a unique combination of focal/other:

 dyad  focal  other  r    present  behavior  btype1  btype2  btype3 ...
 1     10101  10201  .25  3        4         1       1       0
 2     10101  10402  0    1        0         0       0       0      
 3     13201  12811  0    2        7         1       1       1

The dependent/response variable is behavior, and r and present are predictor/independent variables. The complication is that there are different subcategories within behavior, which are coded with btype 1-5 (dummy variables).

Rather than subdivide the data into separate samples and rerun the same models, I have instead used these dummy variables to create interaction terms within the MCMCglmm. When I do so, the interaction terms seem interpretable and behave in expected (and theoretically interesting) ways. But I have a few problems/questions about this, that are making me question whether this is an appropriate way to use the data/create the interaction terms.

First, the model automatically adds the basic btype terms into the model when I specify the interaction terms. Since the btypes are coding for presence/absence of a subcategory of behavior, they are of course highly significant predictors with large effect sizes. I understand this is because of the circular nature of the prediction, so I'm not attempting to interpret them in any useful way. Also note that they are not mutually exclusive, since a dyad's total behavior will often include more than one subcategory or btype.My first and primary question is whether this is an appropriate way to use interaction terms.

There are also some practical problems that arise from these terms:

1) The p values (pMCMC) drastically increase for the other terms in the models, even those with very low p-values in models without the interaction terms. (P-values aren't everything of course, but still they are something). The only terms with p <.05 in these models are btype terms, which aren't at all useful to me. P values greater than .5 (<-not a typo) for most of the other values make it difficult to argue that any of the main effects or interaction effects matter much at all.

2) Including these terms changes the DIC value drastically, and because of the circular nature of the terms, I can't use DIC values to make a meaningful comparison to similar models that don't use interaction terms. Is there some other way to draw a comparison?

I've read through Hadfield's Course Notes and Hadfield (2010), as well as some other materials on MCMCglmms and related modeling, which is where I came up with this crazy model in the first place. However, I'm a beginner in both R and in Bayesian statistics, so your help would be much appreciated: can I use this kind of interaction term? Are there solutions to problems #1 and #2 ?

  • $\begingroup$ Out of curiosity, how are you modeling the response variable? In R, this would be what you pass to the family argument. $\endgroup$ Oct 2 '14 at 1:52
  • $\begingroup$ In other words, what is $behavior$ describing? Also, what is the goal of fitting this model? $\endgroup$ Oct 2 '14 at 1:54
  • $\begingroup$ It's a zero-inflated poisson distribution, because they are count data with many "false" zeros. The goal of the model is to predict the prevalence of teaching behavior among dyads, based on relatedness and some other predictors. So behavior is total # of teaching behaviors observed. $\endgroup$
    – M.A.Kline
    Oct 2 '14 at 17:05

I've had a discussion with a stats consultant on this topic - he says I should use a multinomial model instead of interaction terms. If you're having a similar issue there are examples in Hadfield's course notes (chapter 5): http://cran.r-project.org/web/packages/MCMCglmm/vignettes/CourseNotes.pdf

I'll update further when I figure out the multinomial model.

  • $\begingroup$ Further update: a multinomial model is appropriate only if your frequencies add up to 1 across the outcomes. In my case, the 5 outcomes are not mutually exclusive, so I need to use a multiresponse model, but not a multinomial model. Hadfield's course notes don't appear to cover this. $\endgroup$
    – M.A.Kline
    May 1 '18 at 20:26

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.