R help - Sep 2006 - Problem with Variance Components (and general glmm confusion)

Dear list, 

I am having some problems with extracting Variance Components from a
random-effects model:

I am running a simple random-effects model using lme: 

model<-lme(y~1,random=~1|groupA/groupB)

which returns the output for the StdDev of the Random effects, and model AIC etc
as expected.

Until yesterday I was using R v. 2.0, and had no problem in calling the variance
components of the above model using VarCorr(model), together with their 95%
confidence intervals using intervals() - although for some response variables a
call to intervals() returns the error: Cannot get confidence intervals on
var-cov components: Non-positive definite approximate variance-covariance.

I have now installed R v. 2.3.1 and am now experiencing odd behaviour with
VarCorr(lme.object), with an error message typically being returned:

Error in VarCorr(model) : no direct or inherited method for function
'VarCorr' for this call

Is this known to happen? For instance could it be due to the subsequent loading
of new packages? (lme4 for instance?).

To get around this problem I have tried running the same model using lmer: 

model2<-lmer(y~1 + (1|groupA) + (1|groupB))

Should this not produce the same model? The variance components are very similar
but not identical, making me think that I am doing something wrong. I am also
correct in thinking that intervals() does not work with lmer? I get: Error in
intervals(model2) : no applicable method for "intervals"

GLMM

I have a general application question - please excuse my ignorance, I am
relatively new to this and trying to find a way through the maze.  In short I
need to compile generalized linear mixed models both for (a) Poisson data and
(b) binonial data incorporating a two nested random factors, and I need to be
able to extract AIC values as I am taking an information-theoretic approach to
model selection.  Prior to sending an email to the list I have spent quite a few
days reading the background on a number of functions, all of which offer
potential for this; glmmML, glmmPQL, lmer, and glmmADMB.  I can understand that
glmmPQL is unsuitable because there is no way of knowing the maximised
likelihood, but is there much difference between the remaining three options? I
have seen simulation comparisons published on this list between glmmADMB and
glmmPQL and lmer, but it seems these are before the latest release of lmer, and
also they do not evaluate glmmML.  To a newcomer this myriad of options is
bewildering, can anyone offer advice as to the most robust approach?

Many thanks for your time and patience, 

Toby Gardner 

School of Environmental Sciences
University of East Anglia
Norwich, NR4 7TJ
United Kingdom
Email: t.gardner@uea.ac.uk
Website: www.uea.ac.uk/~e387495

	[[alternative HTML version deleted]]

On 9/4/06, Toby Gardner <t.gardner at uea.ac.uk>
wrote:
> Dear list,
>
> I am having some problems with extracting Variance Components from a
random-effects model:
>
> I am running a simple random-effects model using lme:
>
> model<-lme(y~1,random=~1|groupA/groupB)
>
> which returns the output for the StdDev of the Random effects, and model
AIC etc as expected.
>
> Until yesterday I was using R v. 2.0, and had no problem in calling the
variance components of the above model using VarCorr(model), together with their
95% confidence intervals using intervals() - although for some response
variables a call to intervals() returns the error: Cannot get confidence
intervals on var-cov components: Non-positive definite approximate
variance-covariance.
>
> I have now installed R v. 2.3.1 and am now experiencing odd behaviour with
VarCorr(lme.object), with an error message typically being returned:
>
> Error in VarCorr(model) : no direct or inherited method for function
'VarCorr' for this call
>
> Is this known to happen? For instance could it be due to the subsequent
loading of new packages? (lme4 for instance?).
Yes.  Avoid loading lme4 and nlme simultaneously.
>
> To get around this problem I have tried running the same model using lmer:
>
> model2<-lmer(y~1 + (1|groupA) + (1|groupB))
In recent versions of lme4 you can use the specification

model2 <- lmer(y ~ 1 + (1|groupA/groupB))

Your version may be correct or not.  It depends on what the distinct
levels of groupB correspond to.  The version with the / is more
reliable.
>
> Should this not produce the same model? The variance components are very
similar but not identical, making me think that I am doing something wrong. I am
also correct in thinking that intervals() does not work with lmer? I get: Error
in intervals(model2) : no applicable method for "intervals"
That is correct.  Currently there is no intervals method for an lmer
model.  You can use mcmcsamp to get a Markov chain Monte Carlo sample
to which you can apply HPDinterval from the "coda" package.  However,
these are stochastic intervals so it is best to try on a couple of
chains to check on the reproducibility or the intervals.
>
> GLMM
>
> I have a general application question - please excuse my ignorance, I am
relatively new to this and trying to find a way through the maze.  In short I
need to compile generalized linear mixed models both for (a) Poisson data and
(b) binonial data incorporating a two nested random factors, and I need to be
able to extract AIC values as I am taking an information-theoretic approach to
model selection.  Prior to sending an email to the list I have spent quite a few
days reading the background on a number of functions, all of which offer
potential for this; glmmML, glmmPQL, lmer, and glmmADMB.  I can understand that
glmmPQL is unsuitable because there is no way of knowing the maximised
likelihood, but is there much difference between the remaining three options? I
have seen simulation comparisons published on this list between glmmADMB and
glmmPQL and lmer, but it seems these are before the latest release of lmer, and
also they do not evaluate glmmML.  To a newcomer this myriad !
>  of options is bewildering, can anyone offer advice as to the most robust
approach?
Goran can correct me if I am wrong but I don't believe that glmmML can
be used with multiple levels of random effects.

I'm not sure what the status of glmmADMB is these days.  There was
some controversy regarding the license applied to some of that code a
while back.  I don't know if it has been resolved to everyone's
satisfaction.

When using lmer I would suggest using method = "Laplace" and perhaps
control = list(usePQL = FALSE, msVerbose = 1) as I mentioned in
another reply to the list a few minutes ago.

Let us know how it works out.
>
> Many thanks for your time and patience,
>
> Toby Gardner
>
> School of Environmental Sciences
> University of East Anglia
> Norwich, NR4 7TJ
> United Kingdom
> Email: t.gardner at uea.ac.uk
> Website: www.uea.ac.uk/~e387495
>
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>