I have seven years of continuous insect population data, along with temperature and humidity parameters. I’d like to use this data to predict future populations in a given year using a generalized additive model (GAM) with package mgcv.

Note: This data is count data and contains many zeroes, so I'm using a negative binomial distribution.

My model looks like this:

model <- gam((population~s(temperature)+s(humidity)+s(date)+
         +s(trap, bs="re")+s(site, bs="re")+s(year, bs="re")), family=nb(),data=df)

I used delta AIC values to determine this as the best model. In the above, trap, site, and year are all random effects. Year is listed as a random effect because populations within each year are variable (largely due to temperature within that year) and are being treated as replication (i.e. independent). I'd like to predict the population based on the smoothing factors of temperature, humidity, and date. Since I want to predict the population in any given year, I've compiled my date data into one year, and then treat the year and a random effect (see graph below).

enter image description here

As I mentioned, I would like to predict populations across, say, 2018. (And ideally, feed this model current temperature/humidity data to-date)

I've tried mapping the model predictions on my existing data, and get something like this.

 fpred<- predict(model, se=TRUE, type="response")
plot(df$date, df$population, cex=1.1, pch=16,
     main="Negative binomial GAM",xlab="Date",ylab="Average popuation") 
I<- order(df$date)
lines(df$date[I], fpred$fit[I], lwd=2) 

In the image below, Black dots represent the insect population throughout the season, and the red lines indicate the predicted population using my GAM.

Black dots represent the insect population throughout the season, and the red lines indicate the predicted population using my GAM

Obviously this isn't accurate, and I think I'm missing something when is comes to GAMs. Any ideas or guidance on how to construct and predict with GAMs would be greatly appreciated!

Here's what happens when I do day of year, instead of date (I thought date itself may be causing some errors, but sadly this wasn't the fix for me). enter image description here

And on a side note, this is what happens when I simplify my model to just look at population ~ dateofyear(doy).

modelb <- gam((population~s(doy)), 
fpred<- predict(modelb, se=TRUE, type="response")
plot(dfe$doy, dfe$population, cex=1.1, pch=16,
     main="Negative binomial GAM",xlab="Day of year",ylab="Average population") 
I<- order(dfe$doy)
lines(dfe$doy[I], fpred$fit[I], lwd=2, col="red")

enter image description here

My day of year spline from the full model and gam.check results look like this: enter image description here enter image description here

Here's the histogram of my population variable.: enter image description here

  • 3
    $\begingroup$ Why are you fixing the theta parameter of the NegBin at 1? Use family = nb() instead and let gam() estimate theta. I doubt you want year to be random; any none-observed year will just get the same "population"-level prediction modified by whatever you chose as your scenarios for future covariates. Are years essentially independent? I'm not quite following why size = 11.8. You shouldn't be producing a confidence interval using the response scale - you are using the formulation that assumes things are additive and they can't be on the response scale as counts are bounded at 0. $\endgroup$ – Gavin Simpson Apr 5 '18 at 19:49
  • $\begingroup$ What do you mean by "compile" you data into one year? Are there multiple 'years' of data in the plot you show? $\endgroup$ – Gavin Simpson Apr 5 '18 at 19:50
  • $\begingroup$ Thanks for your response! I've updated the model to be family =nb() and I still get the same result. I've taken out the code specifying size = 11.8 (this wasn't really relevant to producing the graph anyways). There are multiple "years" of data in the plot I show - another graph has been added to illustrate that. $\endgroup$ – Heather Apr 5 '18 at 20:10
  • 1
    $\begingroup$ What's date? Is that plot essentially showing day of year? Are those real zeroes out in the left end of the x-axis? $\endgroup$ – Gavin Simpson Apr 5 '18 at 20:35
  • 1
    $\begingroup$ I think it would help if you could show a plot of the distribution of your outcome variable (just use the command plot(table(y)), where y is your outcome variable), as well as a plot of your model residuals versus your fitted values. As a side remark, in many applications, temperature and humidity tend to interact with each other. $\endgroup$ – Isabella Ghement Apr 11 '18 at 15:39

Your Answer

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

Browse other questions tagged or ask your own question.