Re: [S] Varcomp for comparing multistratum models

Prof Brian D Ripley (
Fri, 4 Dec 1998 07:27:04 +0000 (GMT)

On Fri, 4 Dec 1998, Ray Webster wrote:

> I'm using varcomp in S-Plus 4.0 for Windows to compare unbalanced
> multistratum models. In particular, I have some variable measured at
> different times at various sites, but not all sites have a measurement at
> each time.
> I want to test the significance of various independent variables, and I'm
> doing this by comparing the model likelihoods obtained using
> method="reml". My problem is that sometimes the more complex model
> has a lower (more negative) likelihood. An example is:

It is a _restricted_ likelihood, based on fewer effective observations.
Use method="ml" to compare models. I think you will find lme() more
informative for such models, but you will need est.method="ML".

> mod1_varcomp(log(rep)~site,data,method="reml",tol=1e-10)
> mod3_varcomp(log(rep)~plant+site,data,method="reml",tol=1e-10)
> where log(rep) is the dependent var, site is a random factor, and plant is
> a quantitative independent variable (data is the name of the data frame).
> I can get the likelihoods from mod1$info[7] and mod3$info[7]:
> > mod1$info[7]
> objective
> -8.652036
> > mod3$info[7]
> objective
> -8.867806
> This makes no sense to me at all - how can the likelihood for the more
> complex model be less? Is there an error in the function, or am I being
> stupid here?

You need to know the likelihood of what: in this, lme (by default) and
arima.mle (which isn't using MLEs) the effective dataset changes with the
model. REML uses the marginal likelihood of the residuals from a
least-squares fit of the fixed effects (grand mean +/- plant in your

Brian D. Ripley,        
Professor of Applied Statistics,
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

----------------------------------------------------------------------- This message was distributed by To unsubscribe send e-mail to with the BODY of the message: unsubscribe s-news