[R-sig-ME] Jointly test multiple terms in a MCMCglmm model to be zero

[Paul Johnson]

Hi Vincenzo,

I think best the way to choose between the two models with MCMCglmm is to compare DICs. There's a section in the MCMCglmm vignette on DIC.

set.seed(1234)

model.mcmc1 <- MCMCglmm(weight ~ Time + Diet + Time : Diet, data = BodyWeight, random= ~ Time, verbose=FALSE)

model.mcmc1$DIC
[1] 1750.806

model.mcmc2 <- MCMCglmm(weight ~ Time + Time : Diet, data = BodyWeight, random= ~ Time, verbose=FALSE)

model.mcmc2$DIC
[1] 1987.564

model.mcmc2$DIC - model.mcmc1$DIC
[1] 236.7585

The model with the main effect of Diet has much lower DIC so is much the better model. Although it seems odd to test for the Diet terms being zero while allowing the interaction between Diet and Time, which is effectively what you're doing here and in the lme and lmer examples.

For more guidance see Spiegelhalter et al. "Bayesian measures of model complexity and fit" (ref copied below):

"9.2.4. What is an important difference in DIC?
Burnham and Anderson (1998) suggested models receiving AIC within 1-2 of the 'best' deserve
consideration, and 3-7 have considerably less support: these rules of thumb appear to work
reasonably well for DIC. Certainly we would like to ensure that differences are not due to
Monte Carlo error: although this is straightforward for ¯D, Zhu and Carlin (2000) have explored
the difficulty of assessing the Monte Carlo error on DIC."

Paul

@article {RSSB:RSSB353,
author = {Spiegelhalter, David J. and Best, Nicola G. and Carlin, Bradley P. and Van Der Linde, Angelika},
title = {Bayesian measures of model complexity and fit},
journal = {Journal of the Royal Statistical Society: Series B (Statistical Methodology)},
volume = {64},
number = {4},
publisher = {Blackwell Publishers},
issn = {1467-9868},
url = {http://dx.doi.org/10.1111/1467-9868.00353},
doi = {10.1111/1467-9868.00353},
pages = {583--639},
keywords = {Bayesian model comparison, Decision theory, Deviance information criterion, Effective number of parameters, Hierarchical models, Information theory, Leverage, Markov chain Monte Carlo methods, Model dimension},
year = {2002},
}

-----Original Message-----
From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of vlagani at ics.forth.gr
Sent: 15 June 2011 13:54
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] Jointly test multiple terms in a MCMCglmm model to be zero

Hello listmembers,

I am interested in testing whether multiple coefficients of a MCMCglmm
model are all zero. In particular, I have a categorical term in the
model with four different values, and I want to check if this term is
significant or not (that means, I want to check if *any* of the three
coefficients encoding this term is statistically different from zero).

Please accept my apologizes if this question has already been
answered, but I spent two days looking for an answer across mailing
lists, MCMCglmm help pages and vignettes, and I have not found
anything (maybe I am just not so good in searching!)

So, here a small example code:

#loading the libraries
require(MCMCglmm)
require(nlme)
require(lme4)

#setting the contrasts
contrasts(BodyWeight$Diet) = contr.sum(3)

#fitting the models: I am interested in the Diet term, that is a categorical
#variable with 3 values
model.lme <- lme(weight ~ Time * Diet, data = BodyWeight, random= ~ Time)
model.lmer <- lmer(weight ~ Time * Diet +(Time|Rat) , data = BodyWeight)
model.mcmc <- MCMCglmm(weight ~ Time * Diet, data = BodyWeight,
random= ~ Time, verbose=FALSE)

#checking the significance of the terms: the anova function will provide a
#p-value for the Diet term in the lme model, and an F value for the lme4 model
#How to obtain a similar statistic/pvalue for the MCMC model?
anova(model.lme, type = 'm')
anova(model.lmer, type = 'm')

Thanks in advance,

Vincenzo

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


The University of Glasgow, charity number SC004401