2
$\begingroup$

I am modeling a three-way interaction with a random error in GAM using mgcv package. The data is from a plant growth experiment with 2 categorical predictors f1 (mois; moisture = 3 levels) and f2 (spec; plant species = 4 levels), a continuous x1 predictor (temp_cont; temperature = 5 levels of a warming treatment treated as continuous in the model), and block, f3 (4 levels).

Datapoints with geom_smooth(loess) line for temperature - response_var relationship at different factor levels: Raw data

The aim is to determine (1) whether the relationship of plant growth (and carbon fluxes) with temperature is nonlinear or not (we anticipate it to be non-linear) and (2) whether the shape (slope) of the temperature response curves differ between plant species.

I would first fit a three-way interaction model and then compare this model to a model without the 'species' term using the models' AICs.

I struggle with fitting the model which also accounts for the block effect, and also how to check if this model fits the data well, and after reading many online posts, my latest try is:

gam(response_var ~ spec*moist+
                s(temp_cont, by = spec, k = 5, bs = 'tp')+
                s(temp_cont, by = moist, k = 5, bs = 'tp')+
                s(temp_cont, by = spmo, k = 5, bs = 'tp')+
                s(block, bs = 're'),
                data = data_var,
                method = 'REML')

spmo is the interaction term between spec and moist (data_var %>% mutate(spmo = as.factor(interaction(spec, moist, drop = TRUE))))

Is this model correctly and sufficiently taking into account the block effect? When I plot the smooths for the different levels of combinations of the interacting categorical factors, some of the plots look strange and edf numbers does not reflect the non-linearity of those relationships plotting the raw data. I wonder whether the three-way interaction model requires the two-way interactions (i.e., s(temp_cont, by = spec, k = 5, bs = 'tp') and s(temp_cont, by = moist, k = 5, bs = 'tp')).

The original model as per the script above: the original model as per the script above

Updated model without the two-way interactions: updated model without two-way interactions

I have read about options to account for different levels of wiggliness for the levels of a categorical factor included in the model and am not certain that I understood it correctly. What will be the way to determine whether different levels of wiggliness for a categorical factor in te model are needed?

Then, based on the edf parameter of the model, I can determine the linearity of relationships between temperature and response variable at different levels of a categorical factor (whether it is the three-way interaction or the reduced model).

I would also compare such a model with the same linear model (lme4 model) using AIC to further show that GAM would be preferred over a linear model, if this is the case. Will such a comparison be valid? I have read that the AICs are differently calculated for the GAM vs. the lme4 model.

The other issue is that the overall model adjusted.r2 is 0.9 (for the three-way interaction gam model shown above), which suggests overfitting, but we are not concerned about this as we are not trying to predict the response for any other plant species and we want to fit the model to this data and asses whether the relationship is nonlinear or not (non-linearity is the novelty of our research and we aim to show that it exists).

The link to the original data: https://drive.google.com/file/d/1HiDbZULJeDSOw8zbmseT2iqR513JQLni/view?usp=drivesdk

I would greatly appreciate constructive comments and recommendations for my model improvement. Many thanks.

--------------------------additional questions after the initial response--------------------------

The location-scale Tweedie with the power parameter allowed to vary with the 'moist' factor is a more appropriate choice than without this allowance, as you have shown in your initial response. I'm curious about the decision-making process in determining which factor should govern the variation of the power parameter ('moist' or 'spec' in this case). To explore this, I compared model, designated as m5 in your response, with another model where the power parameter varied with the 'spec' factor (and smoothly with temp_cont). After comparing both models using their AIC values, your selected model exhibited a lower AIC than the other model (518.9 vs. 522.6). Is this method consistent with your approach for selecting the factor influencing the power parameter of the Tweedie distribution? Alternatively, is there any indication on the selection of the factor when screening the plots produced by plotting the model?

My overall aim is to determine whether the response variable exhibits a nonlinear relationship with temperature and whether additional predictor variables ('spec' and/or 'moist') affect the shape of this relationship (or slope of this relationship, in the case when the linear relationship is preferred). Initially, I fitted the three-way interaction model with all the other 2-way interaction terms, denoted as the m5 model to address this. By systematically evaluating AIC values for models with individual interaction terms removed, I determined that the three-way interaction model (m5 model) provided the best fit, retaining all the other terms. Am I then right to consider the m5 model as the more appropriate model?

You also showed the GAM model with the Gamma family which has the shape parameter allowed to vary with moist and smoothly in temp_cont(model m7) has the lowest AIC (although within 2 AIC units of the AIC of m5). Which of the models do you prefer and why?

I recently learned online that an effective degree of freedom (edf) equal to 1 or lower indicates a linear relationship, while edf greater than 1 suggests nonlinearity. The model parameters of model m5 include 'edf', 'edf1', and 'edf2', with 'edf2' being identical to 'edf' values. Upon extracting the 'edf' values from the model, none exceeded 1. Does this simply imply that all relationships within the model are linear and not wiggly?

Additionally, I fitted a linear model using gam function, which will be the linear equivalent of the m5 model and compared AIC values between this model and the m5 model.

m5_linear <- gam(list(resp.var ~ spec * moist * temp_cont+
                   s(block, bs = 're'),
                   ~ 1,
                   ~ moist + s(temp_cont, k = 5)),
            data = data_var,
            family = twlss(),
            method = 'REML')

This linear model exhibited a higher AIC compared to the GAM model (545 vs. 518, respectively). Have I made this comparison correctly? What does the discrepancy between the fact that AIC would prefer the GAM model wile all 'edf' values are equal to 1 or lower signify?

I'm uncertain about which parameters to consider when assessing GAM model fit. For the m5 model, I've examined the normal Q-Q plot, which displays deviations from a straight line at both ends.

enter image description here

In my interpretation, these deviations might still be acceptable. Would you agree with this assessment?

Which other model checks should be done?

The 'k.check' function gives this output which does not look like standard output from the k.check function with smooth interaction terms:

                        k'       edf  k-index p-value
s(temp_cont)             4  3.520319 1.139948  0.9775
s(temp_cont,spec)       15  7.000248       NA      NA
s(temp_cont,moist)      10  4.059748       NA      NA
s(temp_cont,spec,moist) 30 13.342120       NA      NA
s(block)                 4  2.567347       NA      NA
s.2(temp_cont)           4  2.206703 1.139948  0.9875
$\endgroup$

1 Answer 1

4
+100
$\begingroup$

A problem with the first formulation is that there is nothing in the way the factor by smooths are set up that makes the smooths generated by the different terms orthogonal to one another and that can have implications for model identification.

A better way might be to fit your model as

library("readr")
library("mgcv")
library("dplyr")
library("gratia")

data_var <- readr::read_csv("~/Downloads/data_var.csv") |>
  mutate(block = factor(block),
  spec = factor(spec),
  moist = factor(moist))

m <- gam(response_var ~ s(temp_cont, k = 5) + 
    s(temp_cont, spec, k = 5, bs = 'sz')+
    s(temp_cont, moist, k = 5, bs = 'sz')+
    s(temp_cont, spec, moist, k = 5, bs = 'sz')+
    s(block, bs = 're'),
  data = data_var,
  family = Gamma(link = "log"),
  method = 'REML')

which does two things:

  1. we create a series of smooth-factor interactions that are orthogonal to lower order effects (i.e. the first order interaction is orthogonal to the main effect smooth s(temp_cont), and the second order interaction is orthogonal to the first order and main effects smooths), and
  2. we fit the model with a more appropriate random component via family = Gamma(link = "log")

This model produces

enter image description here

and

> summary(m)

Family: Gamma 
Link function: log 

Formula:
response_var ~ s(temp_cont, k = 5) + s(temp_cont, spec, k = 5, 
    bs = "sz") + s(temp_cont, moist, k = 5, bs = "sz") + s(temp_cont, 
    spec, moist, k = 5, bs = "sz") + s(block, bs = "re")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.02427    0.03801   26.95   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                           edf Ref.df       F  p-value    
s(temp_cont)             3.263  3.708 180.002  < 2e-16 ***
s(temp_cont,spec)        7.049  7.701  18.191  < 2e-16 ***
s(temp_cont,moist)       5.428  6.128 112.672  < 2e-16 ***
s(temp_cont,spec,moist) 13.052 13.633   2.988 0.000432 ***
s(block)                 2.356  3.000   3.697 0.003353 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.896   Deviance explained = 87.8%
-REML =  326.5  Scale est. = 0.073748  n = 238

I used Gamma(link = "log") because your concentration response can't be negative and therefore has a mean-variance relationship and the Gamma is a reasonable starting point to model such data.

What will be the way to determine whether different levels of wiggliness for a categorical factor in te model are needed?

You already did that; by definition, the by smooths and the sz smooths I use allow for different wigglinesses as each smooth gets it's own smoothness parameter. If you want them to share a single smoothness parameter than you need to use the id = 1 for example:

m2 <- gam(response_var ~ s(temp_cont, k = 5) + 
    s(temp_cont, spec, k = 5, bs = 'sz', id = 1)+
    s(temp_cont, moist, k = 5, bs = 'sz', id = 2)+
    s(temp_cont, spec, moist, k = 5, bs = 'sz', id = 3)+
    s(block, bs = 're'),
  data = data_var,
  family = Gamma(link = "log"),
  method = 'REML')

which results in far few smoothing parameters being estimated in m2 because I set all the smooths for a give s() call the same id

r$> m$sp                                                                                               
             s(temp_cont)        s(temp_cont,spec)1        s(temp_cont,spec)2 
             8.807430e-02              1.695839e+00              2.810039e+04 
       s(temp_cont,spec)3        s(temp_cont,spec)4       s(temp_cont,moist)1 
             1.803149e+00              7.822437e-05              2.702676e-05 
      s(temp_cont,moist)2       s(temp_cont,moist)3  s(temp_cont,spec,moist)1 
             6.038187e-01              4.021637e+00              3.395644e+04 
 s(temp_cont,spec,moist)2  s(temp_cont,spec,moist)3  s(temp_cont,spec,moist)4 
             2.008389e-01              2.922855e-01              3.173709e-04 
 s(temp_cont,spec,moist)5  s(temp_cont,spec,moist)6  s(temp_cont,spec,moist)7 
             5.505667e-01              3.936865e+00              6.172809e-05 
 s(temp_cont,spec,moist)8  s(temp_cont,spec,moist)9 s(temp_cont,spec,moist)10 
             4.159845e-05              2.151653e+04              1.870382e+01 
s(temp_cont,spec,moist)11 s(temp_cont,spec,moist)12                  s(block) 
             2.412331e+04              2.915978e+04              1.625096e+01 

r$> m2$sp                                                                                              
           s(temp_cont)       s(temp_cont,spec)      s(temp_cont,moist) s(temp_cont,spec,moist) 
           9.183549e-02            3.698062e+00            6.679556e-01            1.990463e+03 
               s(block) 
           1.672849e+01

But doing so results in a worse fit

> AIC(m, m2)
         df      AIC
m  36.03123 577.4602
m2 33.93994 580.5186

I would also compare such a model with the same linear model (lme4 model) using AIC to further show that GAM would be preferred over a linear model, if this is the case. Will such a comparison be valid? I have read that the AICs are differently calculated for the GAM vs. the lme4 model.

You should ideally fit the linear model using {mgcv} and gam() and then compare the two gam() models. You don't say what linear mixed model you want to fit, so I can't help with the gam() equivalent, but it is usually possible to get the same result as lmer() and glmer() with gam() except for the use of correlated random effects.

The other issue is that the overall model adjusted.r2 is 0.9 (for the three-way interaction gam model shown above), which suggests overfitting

Your model is not over fitted, nor is my version of it.

That's not to say there aren't problems - it seems the implied mean variance relationship of the Gamma is not appropriate for these data and a Tweedie random component does a better job (the variance scales as $\mu^{1.2}$ in the models below while the Gamma assumes it scale as $\mu^{2}$, for example.

# Fit as a Tweedie instead of a Gamma GAM
m3 <- gam(response_var ~ s(temp_cont, k = 5) + 
    s(temp_cont, spec, k = 5, bs = 'sz')+
    s(temp_cont, moist, k = 5, bs = 'sz')+
    s(temp_cont, spec, moist, k = 5, bs = 'sz')+
    s(block, bs = 're'),
  data = data_var,
  family = tw(link = "log"),
  method = 'REML')

# Fit a location-scale Tweedie, where we allow the power parameter
# to vary with `moist`
m4 <- gam(list(response_var ~ s(temp_cont, k = 5) + 
    s(temp_cont, spec, k = 5, bs = 'sz')+
    s(temp_cont, moist, k = 5, bs = 'sz')+
    s(temp_cont, spec, moist, k = 5, bs = 'sz')+
    s(block, bs = 're'),
    ~ 1,
    ~ moist),
  data = data_var,
  family = twlss(),
  method = 'REML')

# Fit a location-scale Tweedie, where we allow the power parameter
# to vary with `moist` and smoothly in `temp_cont`
m5 <- gam(list(response_var ~ s(temp_cont, k = 5) + 
    s(temp_cont, spec, k = 5, bs = 'sz')+
    s(temp_cont, moist, k = 5, bs = 'sz')+
    s(temp_cont, spec, moist, k = 5, bs = 'sz')+
    s(block, bs = 're'),
    ~ 1,
    ~ moist + s(temp_cont, k = 5)),
  data = data_var,
  family = twlss(),
  method = 'REML')

Each of the Tweedie models fits better than the Gamma

> AIC(m, m2, m3, m4, m5)
         df      AIC
m  36.03123 577.4602
m2 33.93994 580.5186
m3 39.47912 535.0788
m4 36.38935 522.2200
m5 37.75694 518.8516

although I didn't try similar location-scale-shape models for the Gamma, which could be done via the gammals() family:

# Fit a location-scale-shape Gamma, where we allow the shape parameter
# to vary with `moist`
m6 <- gam(list(response_var ~ s(temp_cont, k = 5) + 
    s(temp_cont, spec, k = 5, bs = 'sz')+
    s(temp_cont, moist, k = 5, bs = 'sz')+
    s(temp_cont, spec, moist, k = 5, bs = 'sz')+
    s(block, bs = 're'),
    ~ moist),
  data = data_var,
  family = gammals(),
  method = 'REML')

# Fit a location-scale-shape Gamma, where we allow the shape parameter
# to vary with `moist` and smoothly in `temp_cont`
m7 <- gam(list(response_var ~ s(temp_cont, k = 5) + 
    s(temp_cont, spec, k = 5, bs = 'sz')+
    s(temp_cont, moist, k = 5, bs = 'sz')+
    s(temp_cont, spec, moist, k = 5, bs = 'sz')+
    s(block, bs = 're'),
    ~ moist + s(temp_cont, k = 5)),
  data = data_var,
  family = gammals(),
  method = 'REML')

and it does seem that allowing the shape parameter of the Gamma to vary with moist and smoothly in temp_cont that we get the better fit and model diagnostics:

> AIC(m, m2, m3, m4, m5, m6, m7)
         df      AIC
m  36.03123 577.4602
m2 33.93994 580.5186
m3 39.47912 535.0788
m4 36.38935 522.2200
m5 37.75694 518.8516
m6 39.80994 524.5466
m7 42.57224 517.0190

I'm not suggesting you go that far, but I think you are worrying about the wrong things with the concern about whether you are over fitting, whether to use varying or the same smoothing parameters etc, when there are bigger problems such as your fitting a model that is objectively wrong for the type of response you are modelling - assuming it is conditionally Gaussian when you can't have negative concentrations is a bigger issue than worrying about how many smoothing parameters you end up using.


A few comments on the additional questions raised by the OP in their edited question.

I'm curious about the decision-making process in determining which factor should govern the variation of the power parameter ('moist' or 'spec' in this case)

Beyond noting that this was a 5 minute quick look at your data and not a proper analysis (that last bit is up to you to do, or I am available at standard consultancy rates ;-)) the main things to note here are

  1. building a distributional model is quite difficult and there often aren't good rules for how to go about building them.
  2. there is much less information in the data about these other distributional parameters than there is about the mean. Hence, all else equal you will typically want to keep the linear predictors for the other distributional terms simpler than the linear predictor for the location/mean parameter.
  3. Matteo Fasiolo's mgcViz package has some tools to help with diagnosing which linear predictors to enchance through predictors.

More specifically, I think I tried moist and spec in both the linear predictors for the $\sigma$ and $p$ parameters, but once moist was included spec didn't seem that useful. A model with moist * spec seems to be a slight improvement over the one I mentioned (m5)

m8 <- gam(list(response_var ~ s(temp_cont, k = 5) +  
    s(temp_cont, spec, k = 5, bs = 'sz')+ 
    s(temp_cont, moist, k = 5, bs = 'sz')+ 
    s(temp_cont, spec, moist, k = 5, bs = 'sz')+ 
    s(block, bs = 're'), 
    ~ 1, 
    ~ moist * spec + s(temp_cont, k = 5)), 
  data = data_var, 
  family = twlss(), 
  method = 'REML')

but it seems this is mainly due to a single species. So, perhaps it is best to keep spec in the model as that reflects the design of the experiment. It also seems to lead to a more precise estimate of the smooth effect of temp_cont in the linear predictor for $p$.

And for completeness, here's the same model with the gammals() fit:

m9 <- gam(list(response_var ~ s(temp_cont, k = 5) +  
    s(temp_cont, spec, k = 5, bs = 'sz')+ 
    s(temp_cont, moist, k = 5, bs = 'sz')+ 
    s(temp_cont, spec, moist, k = 5, bs = 'sz')+ 
    s(block, bs = 're'), 
    ~ 1, 
    ~ moist * spec + s(temp_cont, k = 5)), 
  data = data_var, 
  family = gammals(), 
  method = 'REML')

and there's isn't much to choose between the two, with the twlss() version (m8) performing slightly better

> AIC(m, m3, m4, m5, m6, m7, m8, m9)
         df      AIC
m  36.03123 577.4602
m3 39.47912 535.0788
m4 36.38935 522.2200
m5 37.75694 518.8516
m6 39.80994 524.5466
m7 42.57224 517.0190
m8 49.24810 514.6135
m9 50.75916 516.0505

IIRC adding terms to the linear predictor for $\sigma$ started to lead to computational issues so I just focused on $p$.

I recently learned online that an effective degree of freedom (edf) equal to 1 or lower indicates a linear relationship, while edf greater than 1 suggests nonlinearity.

That's not always strictly true with {mgcv; using select = TRUE we can often get EDF < 1 while have non-linear estimated functions.

I'm not sure what you were extracting but the EDF of the smooths are all > 1 in m5:

> edf(m5)          
# A tibble: 6 × 2
  .smooth                  .edf
  <chr>                   <dbl>
1 s(temp_cont)             3.52
2 s(temp_cont,spec)        7.02
3 s(temp_cont,moist)       4.02
4 s(temp_cont,spec,moist) 13.4 
5 s(block)                 2.57
6 s.2(temp_cont)           2.21

If you're just looking at what is in m5$edf1 or m5$edf2 then what you are looking at are the contributions to the EDF of the model/smooth per model basis function. These should, IIRC, all be <= 1 as we've penalised the coefficients for these basis functions and they can't be greater than 1 as each basis function contributes max 1 degree of freedom to the model - it would be 1 in an unpenalized model. What you need to do is to sum the EDF by smooth, which is what the gratia::edf() function did to generate that output.

Additionally, I fitted a linear model using gam function, which will be the linear equivalent of the m5 model and compared AIC values between this model and the m5 model.

Note that this model isn't completely linear — you left the smooth in the $p$ linear predictor. But otherwise, yes, this would be a linear equivalent of the model.

What does the discrepancy between the fact that AIC would prefer the GAM model wile all 'edf' values are equal to 1 or lower signify?

There is no discrepancy; you don't seem to have computed the EDF for the smooths correctly. Look at the EDF in the output from summary(m5) to see that the EDFs are > 1.

In my interpretation, these deviations might still be acceptable. Would you agree with this assessment?

Not normally, no. But this QQ-plot is assuming that the residuals are Gaussian. Ideally we'd use reference bands to guide our choice of reference quantiles, but the methods available for that don't work with twlss() it seems.

We see somewhat similar deviations for the gammals() model (m9), especially in the upper tail.

enter image description here

However, this discrepancy is well with the range of discrepancies from the 1:1 line we see if we simulate residuals from the model, which by definition must have the correct distribution. As shown by the reference band. I'm not 100% sure why the more nuanced methods aren't working for twlss() - they don't work for all the families in mgcv, but the necessary random number generator is in mgcv.

The 'k.check' function gives this output which does not look like standard output from the k.check function with smooth interaction terms

This is typical. The test that k.check() does doesn't work for bs = "re" and also it seems for the bs = "sz" basis, I suspect this is because of the way this basis is implemented.

For further diagnostics, I would see the mgcViz package.

$\endgroup$
4
  • 1
    $\begingroup$ Thank you very much for your response. I added my further questions to the original post below the dashed line. $\endgroup$ Feb 25 at 22:58
  • 2
    $\begingroup$ holy cow that's a thorough answer. $\endgroup$
    – Ben Bolker
    Feb 25 at 23:02
  • $\begingroup$ @RadimSarlej I've added some responses to the new questions, but if further clarification is needed, consider doing that in new questions that can focus on a specific statistical issue at a time. $\endgroup$ Feb 28 at 9:22
  • $\begingroup$ @Gavin Simpson - I am still digesting it all. To make sure I understand correctly, the model with bs = 'sz' does not include the main effect of factors in the model formula as opposed to a model with bs = 'tp' specified in the smooths. I understand that the main effects are still calculated for such a model with bs = 'sz'. Am I right? $\endgroup$ Mar 7 at 12:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.