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).
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.
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!
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)), family=nb(),data=dfe) 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")