I am using gam from the mgcv package to model the effects of climate on the size of small lakes. I have data from 70 lakes across 30 years. As change in lake area is non-linear over time, I have modeled Year as a smooth term:

gam(Size ~ s(Year, Lake_ID, bs="fs") + Climate)

where Lake_ID is a factor corresponding to each lake and Climateis the annual temperature and precipitation measurements. However, my PACF plot looks like:

enter image description here

Is my interpretation that this shows significant autocorrelation up to lag 6 correct? I have tried to fit the same model using gamm with an autocorrelation term.

gamm(Size ~ s(Year, Lake_ID, bs="fs") + Climate, correlation = corAR1(form=~Year|Lake_ID)

However, this gave the following error: Error in matrix(0, size.cg[i], size.cg[i]) : object 'size.cg' not found.

This error can be reproduced on a toy dataset below, which seems to only occur using the bs="fs" smooth:

dat <- gamSim(4,n=200,scale=2)
gamm1 <- gamm(y~s(f1,fac,bs="fs"), correlation = corAR1(form=~f1|fac),data=dat)

Is there some way to use the factor-smooth with corAR1 in gamm? Any ideas how to address the autocorrelation in the lake dataset using either gam or gamm?


  • $\begingroup$ I think this is a bug, but it occurred to me you could try the t2() specification for the random splines, which I just updated my answer with. $\endgroup$ – Gavin Simpson May 20 at 17:34

This should work, in principle and it seems like there might be some bugs in mgcv::extract.lme.cov2; some variables that need to be set (including size.cg, but others are needed) are not set but later are later referenced.

A key branch seems to be here:

if (n.levels >= start.level || n.corlevels >= start.level) {
    if (n.levels >= start.level) 
        Cgrps <- nlme::getGroups(b, level = start.level)
    else Cgrps <- grps
    Cind <- sort(as.numeric(Cgrps), index.return = TRUE)$ix
    rCind <- 1:n
    rCind[Cind] <- 1:n
    Clevel <- levels(Cgrps)
    n.cg <- length(Clevel)
    size.cg <- array(0, n.cg)
    for (i in 1:n.cg) size.cg[i] <- sum(Cgrps == Clevel[i])
else {
    n.cg <- 1
    Cind <- 1:n

This model/example is either not being captured by the statement in the if() or the else part needs to set more variables.

You should send an email to Simon Wood, at the address in ?bug.reports.mgcv and provide him with the requested info and the reproducible example.

As an alternative, you might try the equivalent t2() construct:

gamm1 <- gamm(y~ s(x1, fac, bs = c("cr", "re"), full = TRUE),
              correlation = corAR1(form = ~ x1 | fac),
              data = dat)

This doesn't work for the example (the model is over-specified) but might get you past the issue you encountered for your data.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.