5
$\begingroup$

I am trying to model behaviour data using the mgcv package. I am using a beta distribution because the data ranges between $[0,1]$. Here is the GAM output:

Family: Beta regression(4.134) 
Link function: logit 

Formula:
FA ~ s(contrastGroupForTable, bs = "cs", by = experimentalGroupForTable, 
    k = 2) + s(animalGroupForTable, bs = "re") + experimentalGroupForTable + 
    rigGroupForTable + sexGroupForTable

Parametric coefficients:
                              Estimate Std. Error z value Pr(>|z|)  
(Intercept)                    -0.6834     0.3655  -1.870   0.0615 .
experimentalGroupForTableeYFP   0.5792     0.3659   1.583   0.1135  
rigGroupForTableRigB(Right)     0.8689     0.3659   2.375   0.0176 *
sexGroupForTableMale            0.4624     0.3659   1.264   0.2063  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                                                           edf Ref.df Chi.sq p-value    
s(contrastGroupForTable):experimentalGroupForTableCasp3  1.991      2   1888  <2e-16 ***
s(contrastGroupForTable):experimentalGroupForTableeYFP   1.993      2   2000  <2e-16 ***
s(animalGroupForTable)                                  19.669     20   1401  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.765   Deviance explained =   81%
-REML = -1694.5  Scale est. = 1         n = 1681

I believe that I am fitting two splines, one for each experimentalGroupForTable, smoothing across contrastGroupForTable. I also have a random effect s(animalGroupForTable, bs = "re") and a couple of fixed effects.

The question I wanted to answer was: does the two GAM models fitted for each experimentalGroupForTable significantly differ from each other for predicting the value of FA? My interpretation now is since the $p$-value for experimentalGroupForTableeYFP is $> 0.05$, at this confidence interval, that my null hypothesis would be rejected. Rather, it looks like rigGroupForTable as a fixed effect has a significant effect to the model fit of the GAM.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
5
$\begingroup$

Without knowing a ton about your data here, here is what I see from your output:

  • You have rightfully included the grouping factor as both a parametric term and a by-factor smooth. This means you have included mean group differences which can already be directly determined by the intercept. We can see in the parametric coefficients that "Casp 3" is the reference as the intercept and "YFP" is the contrast against that intercept, with a conditional mean increase of about $ 0.5792$, which means that the mean for this group in the response is higher. How high that is in practical terms only you can know given there is a lack of context provided by the question posted.
  • While here we see that $p = 0.1135$, rating this as statistically significant or not isn't in and of itself important...there is about an $11$% chance that you would observe this estimate or a more extreme estimate if the null hypothesis is true. The null is almost never really true, so whether you think this distinction is important is up to you. Consider neoFisherian reporting if this is for a journal that cares about such things.
  • Note that $p$-values for smooth terms in GAMs are more approximate than typical $p$-values are concerned (technical discussion about this can be found here). With more "noise" in your data, I imagine $p$-values become less distinct in determining probability. The $\chi^2$ values are quite high, so I would imagine that this is not the case, but some visualization would help determine whether or not your data is fitting correctly.
  • The "rigGroupForTable" coefficient is indeed statistically significant if your alpha cutoff is $a = .05$.
  • Just for reference, this is not "two models" but simply one model with both parametric and smooth terms.

If you want more detailed answers, it may be advisable to provide further context in your question.

Improve this answer
$\endgroup$
3
  • $\begingroup$ I understand your point about that this is not "two models" but simply one model with a by variable smooths - thanks for this point. Since I am here, how do you interpret deviance explained on the output? Is the deviance explained equivalent to the $R^2$ for regressions? $\endgroup$
    – Betelgeuse
    Commented Jan 5 at 17:45
  • 1
    $\begingroup$ It is somewhat analogous to adjusted $R^2$ for a non-Gaussian response such as yours. $\endgroup$
    – Shawn Hemelstrand
    Commented Jan 5 at 18:02
  • $\begingroup$ Nice answer Shawn (+1) $\endgroup$
    – underflow
    Commented Jan 7 at 13:25
5
$\begingroup$

The question I wanted to answer was: does the two GAM models fitted for each experimentalGroupForTable significantly differ from each other for predicting the value of FA?

You can't actually tell that from this output (beyond the trivial difference that the groups have different mean values of the response, even though this difference is not statistically significant). The first two lines in the smooth terms section of the output are Wald-like tests against the null hypothesis of a flat constant function. As such, this output only tells you that each of the two estimated smooth functions differs from a constant, flat function. It doesn't provide any information whatsoever about the differences in the two smooths.

To answer this you'd need to use the model to estimate the difference between the two smooths (see ?difference_smooths for example in my {gratia} package for a way to do this), or reformulate the model such that it contains a difference smooth, by making experimentalGroupForTable an ordered factor (as one option) and fitting:

# make an ordered factor from your factor by variable
my_df <- my_df |>
  transform(expGrpOrd = ordered(experimentalGroupForTable))

# set the contrasts to be the usual treatment contrasts
contrasts(my_df$expGrpOrd) <_ "contr.treatment"

# fit the model, updated to use the ordered factor and a reference smooth
gam(FA ~ s(contrastGroupForTable, bs = "cs") +                     # !! <-- added
      s(contrastGroupForTable, bs = "cs", by = expGrpOrd, k = 2) + # !! <-- changed
      s(animalGroupForTable, bs = "re") +
      expGrpOrd +                                                  # !! <-- changed
      rigGroupForTable +
      sexGroupForTable,
    data = my_df, ....)

Now the first smooth, the one without the by argument will be for the reference group (level) of the factor, while the smooth generated by the s() term with the factor by will be a smooth difference from the reference smooth.

There are other ways to do this too; another is to fit the "average" smooth effect and then separate difference smooths one per level of the factor, which estimate how each group differs from the average. This is done using the sz basis, and importantly we don't want to include the factor by term as a parametric effect in this model:

gam(FA ~ s(contrastGroupForTable, bs = "cs") +            # !! <-- added
      s(contrastGroupForTable, experimentalGroupForTable,
        xt = list(bs = "cs"), k = 2) +                    # !! <-- changed
      s(animalGroupForTable, bs = "re") +
      rigGroupForTable +
      sexGroupForTable,
    data = my_df, ....)

In both cases I'm not sure what you hope to achieve by setting k so low?

My interpretation now is since the $p$-value for experimentalGroupForTableeYFP is >0.05, at this confidence interval, that my null hypothesis would be rejected.

This is back-to-front; you fail to reject the null hypothesis at the 95%

Improve this answer
$\endgroup$

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.