I’m trying to fit some gamms using gamm4 in Rstudio with the goal of doing some model comparison, but I’m getting stuck with wrapping my head around GAM(M)s and structuring these models.

To give you a bit of context I aim to answer 1. how imperviousness has affected the species occupancy during my study period (ten years) and 2. which method of measuring imperviousness (between ndvi, ndbi and urban land cover percentage [ULC]) is best to reflect the impact of urbanisation on this species? I’ve done this with glms but I’d like more flexibility than these models provide.

For my data:
Binary response variable: occurrence (species occurrence; 0=244, 1=233)
Predictors: measurement of imperviousness (ndvi, ndbi, ULC) and year (0:10)
Random effects: siteID (n=67) nested within catchments (n=5)

I plan on tackling the second question with GAM(M) cross-validation (either leave-one-out or k-fold).
For some reason my results are indicating that the relationships are linear, which I know is not true from my glms and from visualising the occurrence probability across the years.

Here’s the code I’ve been using

# packages
library(mgcv)
library(gamm4)

# Load and format species data
mydata<-read.table("00_Data/mydata.txt",header=T)
mydata$siteID<-as.factor(mydata$siteID)
mydata$catchment<-as.factor(mydata$catchment)
mydata$year<-mydata$year-min(mydata$year)
head(mydata);dim(mydata)

# Fit gamm with year as a smooth term
> gammmod1 <- gamm4(occurrence~s(year, k=10), random=~(1|siteID/catchment), data=mydata, family=binomial)
> summary(gammmod1$gam)
Family: binomial 
Link function: logit 

Formula:
occurrence ~ s(year, k = 10)

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.1858     0.2305  -0.806     0.42

Approximate significance of smooth terms:
        edf Ref.df Chi.sq p-value
s(year)   1      1  1.096   0.295

R-sq.(adj) =  -0.00317   
glmer.ML = 444.03  Scale est. = 1         n = 477


# Fit gamm without smooth term
> gammmodb <- gamm4(occurrence~year, random=~(1|siteID/catchment), data=mydata, family=binomial)
> summary(gammmodb$gam)
Formula:
occurrence ~ year

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.69740    0.52779  -1.321    0.186
year         0.08059    0.07698   1.047    0.295


R-sq.(adj) =  -0.00317   
glmer.ML = 444.03  Scale est. = 1         n = 477

# Fit gamm with ndvi as smooth term
gammmod2 <- gamm4(occurrence~s(year,k=15)+s(ndvi,k=15), random=~(1|siteID/catchment),data=mydata, family=binomial)
summary(gammmod2$gam)
Error in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots) : 
  A term has fewer unique covariate combinations than specified maximum degrees of freedom

I’m still learning about GAM(M)s so any advice or paper recommendations would be greatly appreciated.