[R] Help with gamm errors

[Rob Robinson]

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,co
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" =====