[R] Simplification of Generalised Linear mixed effects models using glmmPQL

[Andrew Robinson]

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/