Using GAM Gaussian family to Examine Significance in Data Pattern Across Two Categorical Variables I came across this post in using GAM with the default Gaussian family to examine the significant difference of the patterns across two categorical variables
Statistical differences between two hourly patterns
I am wondering why is Gaussian used here, instead of Gamma family since MET (metabolic rate) seems to not include negative values ?
Is it because it's due to the purpose of the task, which is not to predict, but rather examining the significance of the pattern across 2 categorical variables?
In other words, is the default Gaussian family fine to use when comes to examining the significance of the pattern of the data across 2 categorical variables when we do not worry about the predictions at all? We are just concern about the significance in the PATTERN of the data. And hence, the family used is not a concern as long as the R-square is reasonably good?
I am new to GAM and would appreciate some advice. Thank you.

The two most important questions you had here are below:

Why is Gaussian used here, instead of Gamma family since MET (metabolic rate) seems to not include negative values ?

Unfortunately we can't read the minds of other posters here, but the logic that gamma should be used purely because of positive values is not that useful an assertion in my opinion. Consider these two distributions I just simulated from gamma and Gaussian families, which both only have positive values:

Only the shape of the distribution matters and we can only infer that from the data.

The family used is not a concern as long as the R-square is reasonably good?

A great example of where this is problematic is a logistic GAM, which is constrained to only have values between 0 and 1. Running a Gaussian regression with assumed normally distributed error variance normally has a number of problems and should be avoided in a linear model. This changes to a degree with a GAM model, but still has it's own issues. See below a GAM model I fit to the same data in gamair, one with a logistic model and the other with a Gaussian model:
library(mgcv)
library(gamair)
data("wesdr")

fit.log <- gam(ret
         ~s(dur,k=5),
         data=wesdr,
         family=binomial,
         method="REML")

fit.gauss <- gam(ret
                 ~s(dur,k=5),
                 data=wesdr,
                 family=gaussian,
                 method="REML")

If you plot them as is:
par(mfrow=c(1,2))
plot(fit.log,
     main = "Logistic GAM")
plot(fit.gauss,
     main = "Gaussian GAM")

You'll notice they are remarkably similar:

However, the y-axis for the logistic regression is plotted by its logit link, whereas the Gaussian model to the right is plotted literally by decimal values between 0 and 1. Already we can see that the Gaussian model isn't that interpretable, as we want to know the outcome for what should be a binary response. This isn't made that useful by having logistic data in logit form either, so we can make one subtle change to the plotting to get a meaningful result...by plotting instead by logistic probability using trans=plogis:
par(mfrow=c(1,2))
plot(fit.log,
     trans = plogis,
     main = "Logistic GAM")
plot(fit.gauss,
     main = "Gaussian GAM")

You can see now that the logistic model now has much more meaningful interpretation, as the y-axis has been replaced with predicted probability values:

Now we can make an approximation that the outcome has a certain probability based off the dur variable. For example, we know that a duration of about 10 years of diabetes yields around a 60% chance of getting retinopathy. We can't really get anything from the plot to the right because they are not data that has a any real interpretable value. So in short, specifying your distribution is always important for regression, though the consequences of that decision vary based off what you're doing.
Edit
As pointed out by Gavin Simpson in the comments, you would need to include the intercept into the predicted probability plot. You can do this with the following code:
plot(fit.log,
     trans=plogis,
     shift=coef(fit.log)[1],
     main = "Logistic GAM Shifted Intercept")


