Jointly test multiple terms in a MCMCglmm model to be zero

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

Hi Paul,

thanks a lot for your answer and for the very useful reference. Yes,  
you are right when you says that is odd to test whether the Diet  
coefficients are all zeros while allowing the interaction between Diet  
and Time.

However, my question was more general: how to jointly test whether  
multiple coefficients in a MCMCglmm model are all zero. For example,  
given a model with coefficients Beta0, Beta1, ..., Beta10, I may want  
to test this null hypothesis on a given subset of coefficients:

H0: Beta1 = Beta4 = Beta5 = Beta6 = 0

Let me give another, more practical example. In the model "weight ~  
Time + Diet + Time:Diet" (weight and Time continuous, Diet categorical  
with three values), there are 6 different coefficients: intercept,  
B_time, B_diet1, B_diet2, B_(time:diet1), B_(time:diet2). If one wants  
to assess whether the time has any effect on the body weight in AT  
LEAST ONE diet, the following null hypothesis should be tested:

H0: B_time = B_(time:diet1) = B_(time:diet2) = 0

The anova.lme function provide the option "Terms" for performing such  
type of tests. After reading your email, I guess that the easiest way  
to perform a similar test in MCMCglmm models is the following:

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

#setting the contrasts and the seed
contrasts(BodyWeight$Diet) = contr.sum(3)
set.seed(1234)

#fitting the models
model.lme <- lme(weight ~ Time * Diet, data = BodyWeight, random= ~ Time)
model.mcmc <- MCMCglmm(weight ~ Time * Diet, data = BodyWeight,  
random= ~ Time, verbose=FALSE)
model.mcmc.noTime <- MCMCglmm(weight ~ Diet, data = BodyWeight,  
random= ~ Time, verbose=FALSE)

#checking the significance of the terms: the anova function with the "Terms"
#option will provide a p-value for checking if ANY coefficients  
related to Time
#is different from zero
anova(model.lme, type = 'm', Terms = c('Time', 'Time:Diet'))

#DIC- bases "test"
model.mcmc$DIC - model.mcmc.noTime$DIC

Both the DIC and the anova tests confirm that the couple of terms  
<Time, Time:Diet> is relevant, thus I conclude that at least one of  
the three coefficients is different from zero.

My final intention is to employ both MCMC and nlme models in order to  
strengthen my analysis. In particular, my idea is to compare nlme  
p-values with MCMC statistics; since there is not general agreement  
about how to calculate p-values  for mixed models, I prefer to double  
check my finding.

Ciao,

Vincenzo


Paul Johnson <Paul.Johnson at glasgow.ac.uk> ha scritto:

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

-- 
Vincenzo Lagani
Visiting Researcher
BioInformatics Laboratory
Institute of Computer Science
Foundation for Research and Technology - Hellas

Vincenzo Lagani
Visiting Researcher
BioInformatics Laboratory
Institute of Computer Science
Foundation for Research and Technology - Hellas