GAM: Plotting gam models of different treatments?

Question

I want to figure out the best way to plot gams, but I am getting confused about the best way to do this for my data. Please note that I am not a statistician, so 'stats/coding for dummies' answers are much appreciated.

Briefly, my experiment looked at how flower production changed over time in response to 4 different flower removal treatments that differed by the timing of removal (Early, Peak, Late, Control/None). My data looks like this, where Trt is treatment, RIL is genetic line and Total is total flower production (response).

> head(data2)
   DateNum  Date Tray TrtNum   Trt Treat Plant RIL Primary Total NonPrim Dist Dist_num_tray Dist_num_plant
1        1 43426    1      1 Early    E1     1 128       0     0       0    2            50           66.8
3        3 43428    1      1 Early    E1     1 128       0     0       0    2            50           66.8
5        5 43430    1      1 Early    E1     1 128       0     0       0    2            50           66.8
7        7 43432    1      1 Early    E1     1 128       0     0       0    2            50           66.8
9        9 43434    1      1 Early           1 128       0     0       0    2            50           66.8
11      11 43436    1      1 Early           1 128       0     0       0    2            50           66.8

During the times of flower removal (3 day period), flower count for that entire treatment is 0 for 3 days. I knew this would obviously cause there to be a treatment effect, but I am interested in the response after this removal. So I divided my data into different subsets to compare flowering after damage of one of the treatments to the control (no damage).

Here is an example of a simple gam comparing Control vs Early after damage occurred:

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

I needed to make this to see the results of anova(gam1) to determine if there is a treatment (Trt) effect. However the plot() function (and some other things I've played around with like plotGAM from vortex package) shows only one curve, I think because it is combining the control and the early data into one curve. What I really want to see is two separate gam curves, one for early and one for control.

Should I just use gam1 for testing if there is a treatment effect, but make separate gams for each treatment, then extract the values (predict() function?) and graph those separately? There must be an easier way, or no?

Thank you! Sorry about the length, I just wanted to make sure I was giving enough info

Answer 1

I think there are at least two answers here, and hopefully @GavinSimpson, or anyone more qualified than I in gam theory can provide a better answer.

1. See @GavinSimpson 's answer to a related question. Note how he is nesting using the by term within mgcv::gam.
2. Specify the interaction term using ti smooth. From mgcv documentation:

   ti smooths exclude the basis functions associated with the ‘main effects’ of the marginal smooths, plus interactions other than the highest order specified. These provide a stable an interpretable way of specifying models with main effects and interactions. For example if we are interested in linear predicto f1(x) + f2(z) + f3(x,z), we might use model formula y~s(x)+s(z)+ti(x,z) or y~ti(x)+ti(z)+ti(x,z).