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 foundPresumably 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" =====
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
Rob, Thanks for sending the example...> > 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-- this relates to the fact that `gamm' assumes that the correlation structure is nested within the random effect grouping structure, so in fact the `|garden' in the correlation formula is redundant (see ?gamm on correlation). Clearly this should be dealt with more gracefully... I'll try to do so for the next release. In the meantime just remove the redundant grouping factor. best, Simon> > 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