# Choosing GLM or GAM for my daily flower production data in R?

Thanks for stopping by to read this! I am a first time poster, hopefully I included the correct among of info.

I've struggled back and forth with this problem for awhile - how do I analyze my data- using GAM, GLM (repeated measures??) or something else? I don't know anyone who works with GAMs and I'm a new graduate student with only a very basic understanding of R/stats/coding. Any more in-depth explanations (ie pointing out what the terms mean in a code) would be very helpful for me, but I'll take whatever help I can get at this point!

This is an example of what my data looks like:

   DateNum   Tray  Treatment  Plant RIL Primary Total NonPrim
1     1    Early      1   128       0     0       0
2     1    Early      1   128       0     0       0
3     1    Early      1   128       1     1       0
4     1    Early      1   128       2     2       0
5     1    Early      1   128       3     5       2
6     1    Early      1   128       6     8       2


My research involves looking at the number of flowers produced by my plants each day of the plant's reproductive life. This information forms a flowering curve. So for above, plant 1 begins to flower on day 3. Then I'd come back on day 4, count flowers again, etc., until the plant has finished flowering a few weeks later. I have 840 plants total with 14 genetic lines (RIL) and 4 different treatments (Early, Peak, Late, Control). Primary just means the number of flowers on the first branch, while total includes all the branches, and NonPrim is Total-Primary. I'm more focused on Total for this analysis.

I want to see if the total flowering curves differ between the plants in different Treatments and RILs. I tried to use GAM but I can't seem to get it to work and I don't think it can with my data. My understanding of it is that if you have multiple continuous predictor variables (mm of rain, wind speed, soil moisture, etc) then you can add those terms to the gam model to see if they explain the variance better. So, for example, if I were to have measured height daily, or soil moisture daily, I could add this to the gam model. But I only measured daily flower production.

Based on a reddit suggestion, I tried this: gam1 <- gamm(Total~s(DateNum)+Treatment,random=list(Plant=~1),data=data1)

And the summary gave me this, which didn't make sense at all to me.

  > summary(gam1)
Length Class Mode
lme 18     lme   list
gam 31     gam   list


Then I noticed on the reddit post that they person said to do summary(gam1$lme). I haven't a clue what the$lme does to the output, but the output was very detailed. I don't know what the output means and if it shows that I have a difference in flowering between my treatments.

Linear mixed-effects model fit by maximum likelihood
Data: strip.offset(mf)
AIC      BIC    logLik
162019.7 162083.8 -81001.84

Random effects:
Formula: ~Xr - 1 | g
Structure: pdIdnot
Xr1      Xr2      Xr3      Xr4      Xr5      Xr6      Xr7
StdDev: 31.25945 31.25945 31.25945 31.25945 31.25945 31.25945 31.25945
Xr8
StdDev: 31.25945

Formula: ~1 | Plant %in% g
(Intercept) Residual
StdDev:    3.016979 8.962859

Fixed effects: y ~ X - 1
Value Std.Error    DF  t-value p-value
X(Intercept)    8.319116 0.2330588 21417 35.69534  0.0000
XTrtEarly       0.259773 0.3352030   836  0.77497  0.4386
XTrtLate       -1.028546 0.3342349   836 -3.07732  0.0022
XTrtPeak       -1.269514 0.3411782   836 -3.72097  0.0002
Xs(DateNum)Fx1 -1.803422 1.3364677 21417 -1.34939  0.1772
Correlation:
X(Int) XTrtEr XTrtLt XTrtPk
XTrtEarly      -0.695
XTrtLate       -0.697  0.485
XTrtPeak       -0.683  0.475  0.476
Xs(DateNum)Fx1  0.000  0.004 -0.002 -0.001

Standardized Within-Group Residuals:
Min          Q1         Med          Q3         Max
-2.81852865 -0.58269460 -0.08064353  0.39935691  7.48671038

Number of Observations: 22258
Number of Groups:
g Plant %in% g
1          840


I haven't tried the GLM yet, but I have for other types of data (ie does the date of peak flowering differ between RILs and Treatments (model1 = lme(DateNum~Treatment*RIL, data=data2, random=~1|Tray). Then I did anova1 <- aov(model1, data=data2) and it all made sense for me. I understand this and what the terms mean, and it seems a lot more intuitive to me and I can control for the block effect of each tray that held the plants in the greenhouse. I don't know if I can do this with GAMs.

Does anyone have any insight on this? I'm sorry if this is all over the place.

• re: "I haven't a clue what the $lme does", I haven't used the gamm function specifically, but in R models usually create objects of their model class. Objects are often basically lists with named elements. You can access named elements of an object in R with the $ operator, so gam1$lme gets the lme element of gam1. It appears that this is the element that actually has model output which summary() displays. Re: what model to choose and how to interpret it, that is an extremely broad question that might be better suited for cross validated. Mar 15, 2018 at 21:09 • Here are a couple of short tutorials on GAMM in R that might be helpful: this and this. Also a book on GAMM in R. Mar 15, 2018 at 21:15 • Thanks eipi10 for those tutorials - I haven't seen them before! I will have to follow along with them as I code. Thank you! Mar 22, 2018 at 20:09 ## 1 Answer Re $lme, this is the object returned by fitting the equivalent GAMM via nlme::lme(). Unless you are very familiar with the mixed model representation of splines, this output is going to be of limited use to you. In this representation, splines can be thought of as a random effects, with the range space of the spline basis (the wiggly basis functions) being columns in the random effects model matrix and a variance parameter plays the role of the smoothness parameter in a traditional GAM.

It is likely more useful to look at the equivalent GAM representation. The $gam component contains most of the components that a typical object fitted by mgcv::gam() does. summary(gam1$gam) will give you back Wald-like tests of the smoother in your model.

For such simple random effects, you could fit this as

gam1 <- gam(Total ~ s(DateNum) + Treatment + s(Plant, bs = 're'), data = data1)


using the opposite side of the splines-as-random-effects duality mentioned above.

You may want to have separate smooths per Plant, not just a different intercept:

gam2 <- gam(Total ~ Treatment + s(DateNum, Plant, bs = 'fs'), data = data1)

• Thanks for your response this was very helpful! I have 2 more questions: 1. What if I wanted to include tray as a random effect (block effect of location of plants in greenhouse) like I did in the GLM (model1 = lme(DateNum~Treatment*RIL, data=data2, random=~1|Tray)? Would it be something like gam3 <- gam(Total ~ Treatment + s(DateNum, Plant, bs = 'fs') + s(Tray), data= data1)? That doesn't seem right to me and I had trouble finding examples. 2. Do I need a term to tell gam that I am doing repeated measures (I measured flower production for each plant once everyday for 30 days)? Mar 22, 2018 at 20:07
• If you want a random effect of Tray you use s(Tray, bs = 're'), i.e. you want to special random effect basis (as per the first code block in my answer, where I did the same but for Plant). There are variations of this model if you want a population-level effect of DateNum and Plant-level random smooths: Total ~ Treatment + s(DateNum) + s(DateNum, Plant, bs = 'fs') + s(Tray, bs = 're') for example. Mar 22, 2018 at 22:15
• As for 2., not necessarily. Assuming DateNum is the time variable and assuming you use a large enough basis for these smooths (check with gam.check(), you may not need anything else. You will need to check residuals for autocorrelation, which is not so trivial if your data are irregularly spaced in time but pretty easy if regularly spaced. If the residuals are effectively uncorrelated then you don't need anything further. Mar 22, 2018 at 22:17
• Hi Gavin, it has been some time since I last posted due to working on other projects, but I have a few more questions about modelling since working more on it. When I included the random effect of tray term to my model (+s(Tray, bs='re'), my treatment effect becomes insignificant, and when I remove the tray term it becomes significant again (from p=0.2 to p=0.03). I don't so much care about labelling things as significant or not, more so understanding why this happens. Jun 13, 2018 at 18:56