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

[vlagani at ics.forth.gr]

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