[R] Help with gamm errors

[Simon Wood]

Rob, 

It might be worth upgrading at least mgcv and MASS to the current versions 
(latest mgcv is 1.3-28, just gone to CRAN). Way back I vaguely remember that 
there was an issue with glmmPQL (called by gamm) not picking up correlation 
structure variables from the dataframe, but I can't now find the details and 
this should be resolved now. The second issue I've not seen before: I could 
imagine it occurring if you have any groups of size zero. Could you let me 
know if either of theseproblems  are not resolved by upgrading, please?, and 
I'll dig further.

best,
Simon

On Tuesday 09 October 2007 18:28, Rob Robinson wrote:
> Dear All
>  Hopefully someone out there can point out what I am missing! I have a
> (large, several hundred) dataset of gardens in which over two years the
> presence/absence of a particular bird species is noted each week. I have
> good reason to believe there is a difference between the two years in the
> weekly proportion of gardens and would like to assess this, before going on
> to look in more detail at particular weeks. There is nothing special about
> the particular gardens used (though they are not strictly a random sample)
> and the weekly 'counts' are obviously autocorrelated (with weeks closer
> together being more similar though the correlation may differ between
> gardens). Thus (I suspect) a gamm statement such as:
>
>
> m1=gamm(present~s(week,bs="cc")+s(week,bs="cc",by=y1),random=list(garden=~1
>) ,
>   correlation=corAR1(form=~week|garden),family=binomial,data=count.data2)
>
> is required (y1 is year and the others are self explanatory, all variables
> are in count.data2). This yields the following output
>
> > Maximum number of PQL iterations:  20
> > iteration 1
> > Error in eval(expr, envir, enclos) : object "week" not found
>
> Presumably something is not getting passed internally to glmmPQL (see
> results of traceback())
>
> [large data vector scrolls off screen]
> 6: eval(expr, envir, enclos)
> 5: eval(mcall)
> 4: glmmPQL(y ~ X - 1, random = rand, data = strip.offset(mf), family =
> family,
>        correlation = correlation, control = control, weights = weights,
>        niter = niterPQL, verbose = verbosePQL)
> 3: eval(expr, envir, enclos)
> 2: eval(parse(text = paste("ret$lme<-glmmPQL(", deparse(fixed.formula),
>
> ",random=rand,data=strip.offset(mf),family=family,correlation=correlation,c
>o ntrol=control,",
>        "weights=weights,niter=niterPQL,verbose=verbosePQL)", sep = "")))
> 1: gamm(present ~ s(week, bs = "cc") + s(week, bs = "cc", by = y1),
>        random = list(garden = ~1), correlation = corAR1(form = ~week |
>            garden), family = binomial, data = count.data2)
>
> Question 1: I have followed the example in Simon Wood's excellent GAM book
> in specifying the form of the correlation term, but have I done this right
> or do I need extra information to get it to recognise the week variable?
>
> Leaving that aside, I altered to the form of correlation to
>
>
> m1=gamm(present~s(week,bs="cc")+s(week,bs="cc",by=y1),random=list(garden=~1
>) ,
>   correlation=corAR1(form=~1|garden),family=binomial,data=count.data2)
>
> (ie removing the week term). This model proceeds - to a point...
>
>  Maximum number of PQL iterations:  20
> iteration 1
> iteration 2
> iteration 3
> iteration 4
> iteration 5
> iteration 6
> Error: attempt to select less than one element
>
> Traceback() suggests that the model fits, but that lme can't calculate
> something?
>
> 2: extract.lme.cov2(ret$lme, mf, n.sr + 1)
> 1: gamm(present ~ s(week, bs = "cc") + s(week, bs = "cc", by = y1),
>        random = list(garden = ~1), correlation = corAR1(form = ~1 |
>            garden), family = binomial, data = count.data2)
>
> Question 2: Any hints on what might be causing this? Am I fitting the wrong
> (or too complicated a model)?
>
> Btw if it is relevant I am using mgcv_1.3-19 and R 2.3.1.
>
> Many thanks for any hints on where I should start digging
> Cheers
> rob
>
> *** Want to know about Britain's birds? Try  www.bto.org/birdfacts ***
>
> Dr Rob Robinson, Senior Population Biologist
> British Trust for Ornithology, The Nunnery, Thetford, Norfolk, IP24 2PU
> Ph: +44 (0)1842 750050         E: rob.robinson at bto.org
> Fx: +44 (0)1842 750030         W: http://www.bto.org
> eSafe scanned this email for viruses, vandals and malicious content (!)
>
> ==== "How can anyone be enlightened, when truth is so poorly lit" =====
>
> ______________________________________________
> R-help at r-project.org 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.

-- 
> Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK
> +44 1225 386603  www.maths.bath.ac.uk/~sw283