Solved – Which method is correct? (generalized additive model, mgcv)

data visualizationgeneralized-additive-modelmgcvrtime series

If possible could you please let me know if I am on the right track? I don't know anyone who works with GAMs who I can ask and I would be so appreciative of any help I receive. I use mgcv in R.

My data is daily flower counts ('Total'), with different plant genetic lines ('RIL') and 4 different treatments ('Trt'). I basically want to compare the 4 treatments and see if there is a difference in flower production over time between them. Are the shapes of the flowering curves different?

An example of my gam:

gam1 <- gam(Total ~ Trt * RIL + s(DateNum, k = 9, bs = "fs") + s(Plant, bs = "re"), data=ce1230)

I understand that using anova.gam(gam1) I can see if there is a significant difference between treatments and rils. However when I plot gam1, it shows one curve that represents everything. I was hoping to see a plot with separate curves for each of my 4 treatments. Does this mean that I need to make gamtrt1, gamtrt2, gamtrt3, gamtrt4 and somehow compare the separate gams to each other? or use the original anova(gam1) and just plot each treatment separately? Or am I looking at this the entirely wrong way?

Best Answer

I would first try an ANOVA and look at that F-test. If it rejects that, there is evidence that there is variation between treatments. After that, since your outcome if count data, I would do a poisson regression. Unlike GAMs where it can be hard to test the important of effects, Poisson regression gives you your usually p-values which makes it simple to see which things are important to causing the differences in your data.

https://stats.idre.ucla.edu/r/dae/poisson-regression/

If there is significnet difference between your mean and variance, (poisson says that they should be the same), modify the poisson regression with a Poisson-Gamma model.

Related Solutions

Solved – Comparing GAM smooths in mgcv R package using AIC

It's not flawed, per se, just pointless. I can think of (and have experienced) situations where the specific type of basis (of the two you mention) can result in markedly different fits. However, where this has happened to me, it has usually been solved by increasing k for one of the bases or because the wrong model was fitted. In this instances trivial differences in the basis used were magnified by the real problem (needing larger basis dimension, fitting the right model), not because of any fundamental difference in the performance of the individual basis.

In most situations you are going to see trivial differences in the fits of models fitted with different basis (among standard bases) and these are going to result in trivial differences in AIC. In most cases AIC is going to tell you, therefore, that the model fits are equivalent.

If you are planning on using SCAM models, I might suggest you use P splines in the GAMs as the splines in the scam package are all based on P splines.

GAMs Plotting – Plotting Generalized Additive Models on Response Scale with Multiple Smooth and Linear Terms

If the model contains z then the effect of x estimated by the model is that given z is in the model. Hence the fitted response is the additive sum of the two effects, and we can't talk generally about the estimated values of the response for a range of values of x without also stating the value of z.

For Gaussian models, you can just add on the intercept in plot.gam() to shift around the smooth curve on the y-axis. See argument shift to plot.gam(). This assumes that as per the example x and z are unrelated in the model, and furthermore some value of z (in this case I think 0 as it is a linear term not subject to identifiability constraints).

A more general solution is just to predict from the model at a grid of values over x while holding z constant at some representative value, say its mean or median.

Here's a full example of doing this by hand:

library("mgcv")
library("ggplot2")

set.seed(1)
df <- gamSim()
m <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "REML")

new_data <- with(df, expand.grid(x2 = seq(min(x2), max(x2), length = 200),
                                 x0 = median(x0),
                                 x1 = median(x1),
                                 x3 = median(x3)))

ilink <- family(m)$linkinv
pred <- predict(m, new_data, type = "link", se.fit = TRUE)
pred <- cbind(pred, new_data)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

ggplot(pred, aes(x = x2, y = fitted)) +
  geom_ribbon(aes(ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line()

producing

That script should be fine for any of the standard family options in mgcv, but you'll have to take careful note of what predict() returns for some of the fancier families in mgcv.

Best Answer

Related Solutions

Solved – Comparing GAM smooths in mgcv R package using AIC

GAMs Plotting – Plotting Generalized Additive Models on Response Scale with Multiple Smooth and Linear Terms

Related Question