T.C. Cameron
2007-Feb-20 13:13 UTC
[R] Simplification of Generalised Linear mixed effects models using glmmPQL
Dear R users I have built several glmm models using glmmPQL in the
following structure:
m1<-glmmPQL(dev~env*har*treat+dens, random = ~1|pop/rep, family Gamma)
(full script below, data attached)
I have tried all the methods I can find to obtain some sort of model fit
score or to compare between models using following the deletion of terms
(i.e. AIC, logLik, anova.lme(m1,m2)), but I cannot get any of them to
work.
Yet I see on several R help pages that others have with similar models?
I have tried the functions in lme4 as well and lmer or lmer2 will not
accept my random terms of "rep" (replicate) nested within
"pop"
population.
I have read the appropriate sections of the available books and R help
pages but I am at a loss of where to move from here
data<-read.table("D:\\bgytcc\\MITES\\Data\\Analysis\\test.txt",header=T)
attach(data)
names(data)
m1<-glmmPQL(dev~env*har*treat+dens, random = ~1|pop/rep, family = Gamma)
summary(m1)
anova.lme(m1)
m2<-update(m1,~.-env:har:treat)
anova.lme(m1,m2)###this does not work
AIC(m1)##this does not work
logLik(m1)##this does not work?
##################this does not work
class(m1) <- "lme"
class(m2) <- "lme"
anova.lme(m1,m2)
#################################
m3<-lmer(dev~env*har*treat+dens + (1|pop/rep), family = Gamma)
## this generates an error
Error in lmerFactorList(formula, mf, fltype) :
number of levels in grouping factor(s) 'rep:pop', 'pop'
is too
large
In addition: Warning messages:
1: numerical expression has 1851 elements: only the first used in:
rep:pop
2: numerical expression has 1851 elements: only the first used in:
rep:pop
m4<-lmer(dev~env*har*treat + dens + (1|rep) +(1|pop), family = Gamma,
method = "Laplace")
## this works but it does not give me an anova output with p values
anova(m4)
Analysis of Variance Table
Df Sum Sq Mean Sq
env 1 17858 17858
har 2 879 439
treat 2 2613 1306
dens 1 1016476 1016476
env:har 2 870 435
env:treat 2 1188 594
har:treat 4 313 78
env:har:treat 4 1188 297
........................................................................
............
Dr Tom C Cameron
Genetics, Ecology and Evolution
IICB, University of Leeds
Leeds, UK
Office: +44 (0)113 343 2837
Lab: +44 (0)113 343 2854
Fax: +44 (0)113 343 2835
Email: t.c.cameron at leeds.ac.uk
Webpage: click here
<http://www.fbs.leeds.ac.uk/staff/profile.php?tag=Cameron_TC>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: test.txt
Url:
https://stat.ethz.ch/pipermail/r-help/attachments/20070220/82646d15/attachment.txt
Andrew Robinson
2007-Feb-20 23:43 UTC
[R] Simplification of Generalised Linear mixed effects models using glmmPQL
Hello Tom,
the problem is because R has assumed that pop and rep are integers,
not factor levels. Try:
test <- read.table("test.txt",header=T)
sapply(test, class)
test$pop <- factor(test$pop)
test$rep <- factor(test$rep)
then try fitting the models. Also, there has been substantial
discussion about the production of p-values for mixed-effects models
in R; it's now a FAQ:
http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-are-p_002dvalues-not-displayed-when-using-lmer_0028_0029_003f
The package writer (Professor Bates) recommends the use of mcmcsamp to
obtain inferential information about your model - see
?mcmcsamp
in the lme4 package for examples of its use.
Cheers,
Andrew
ps it's a bad idea to call things "data" :)
On Tue, Feb 20, 2007 at 01:13:25PM +0000, T.C. Cameron
wrote:> Dear R users I have built several glmm models using glmmPQL in the
> following structure:
>
> m1<-glmmPQL(dev~env*har*treat+dens, random = ~1|pop/rep, family >
Gamma)
>
> (full script below, data attached)
>
> I have tried all the methods I can find to obtain some sort of model fit
> score or to compare between models using following the deletion of terms
> (i.e. AIC, logLik, anova.lme(m1,m2)), but I cannot get any of them to
> work.
> Yet I see on several R help pages that others have with similar models?
>
> I have tried the functions in lme4 as well and lmer or lmer2 will not
> accept my random terms of "rep" (replicate) nested within
"pop"
> population.
>
> I have read the appropriate sections of the available books and R help
> pages but I am at a loss of where to move from here
>
>
>
>
data<-read.table("D:\\bgytcc\\MITES\\Data\\Analysis\\test.txt",header=T)
> attach(data)
> names(data)
>
>
>
> m1<-glmmPQL(dev~env*har*treat+dens, random = ~1|pop/rep, family = Gamma)
> summary(m1)
> anova.lme(m1)
> m2<-update(m1,~.-env:har:treat)
> anova.lme(m1,m2)###this does not work
> AIC(m1)##this does not work
> logLik(m1)##this does not work?
>
>
>
> ##################this does not work
> class(m1) <- "lme"
> class(m2) <- "lme"
> anova.lme(m1,m2)
> #################################
>
> m3<-lmer(dev~env*har*treat+dens + (1|pop/rep), family = Gamma)
>
> ## this generates an error
> Error in lmerFactorList(formula, mf, fltype) :
> number of levels in grouping factor(s) 'rep:pop',
'pop' is too
> large
> In addition: Warning messages:
> 1: numerical expression has 1851 elements: only the first used in:
> rep:pop
> 2: numerical expression has 1851 elements: only the first used in:
> rep:pop
>
>
> m4<-lmer(dev~env*har*treat + dens + (1|rep) +(1|pop), family = Gamma,
> method = "Laplace")
> ## this works but it does not give me an anova output with p values
> anova(m4)
> Analysis of Variance Table
> Df Sum Sq Mean Sq
> env 1 17858 17858
> har 2 879 439
> treat 2 2613 1306
> dens 1 1016476 1016476
> env:har 2 870 435
> env:treat 2 1188 594
> har:treat 4 313 78
> env:har:treat 4 1188 297
>
>
>
> ........................................................................
> ............
> Dr Tom C Cameron
> Genetics, Ecology and Evolution
> IICB, University of Leeds
> Leeds, UK
> Office: +44 (0)113 343 2837
> Lab: +44 (0)113 343 2854
> Fax: +44 (0)113 343 2835
>
>
> Email: t.c.cameron at leeds.ac.uk
> Webpage: click here
> <http://www.fbs.leeds.ac.uk/staff/profile.php?tag=Cameron_TC>
>
>
--
Andrew Robinson
Department of Mathematics and Statistics Tel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/