# GAM with opposite outcomes with different families

I'm building GAMs and I have some doubts regarding the family to use. I'm fitting GAMs because I expect some non-linear relationships between the response variable and some covariates. I've checked a lot of different sources but some details are still not clear to me.

Specifically, I have a positive, non-zero response variable, with the following distribution:

At first, and because I didn't check this in detail, I ran a simple model with the default family (gaussian), and the outcome made sense. Then, when I looked at the residuals, I realized they were not that good. Some model details first:

 s_status_gaussian <- gam(inter_k ~ factor(s_status) + s(month, by=factor(s_status)) +
s(factor(ds_2$id), bs="re"), data=ds_2, method = "REML")  Here are the residuals: So I thought about changing the family. From what I've read, looking at the histogram above, the distribution that would fit best is the gamma distribution, or eventually the inverse gaussian. Also the residuals look better: s_status_gamma <- gam(inter_k ~ factor(s_status) + s(month, by=factor(s_status)) + s(factor(ds_2$id), bs="re"),
data=ds_2,
family = Gamma(link = "inverse"),
method = "REML")


Lastly, I also tried log(inter_kill) and using again Gaussian as family, but the residuals are not so good again:

So, from these results, the most correct family to choose would be the gamma distribution, am I right? I guess this question is rather basic, but I would like to make sure my approach is correct because, when using the gamma distribution, the residuals do look better but the results of the model are basically the opposite of what I previously obtained with the Gaussian family... So I would like to make sure I'm choosing the right approach and interpreting the results correctly.

Additionally, why would the results be exactly the opposite when using another family?

Thanks in advance for any help/suggestions!

First, I would like to point out that you don't have to use GAM if you "expect some non-linear relationships between the response variable and some covariates". You can use nonlinear basis functions, like polynomials, with e.g. lm or glm, too.
Second, note that you changed not only the error family to Gamma, but you also changed the link function to "inverse", i.e. you are actually not computing the regression for inter_k but for 1/inter_k. I am not sure what you mean by "the results be exactly the opposite", but that could be the reason.
Apart from that, your approach seems reasonable. Trying a couple of models and checking out the residuals is the way to go. I would recommend to always prepare some 2D plots of the data, using pairs(), to get a feeling for your data. For the comparison of different models, you might want to consider some basic model selection metrics, like e.g. AIC or adjusted R squared.