Is there normality in my data? Which statistical test should I use? I runned two GLMs using the same dependent and independent variables, but modelling each analysis according a different type of distribution. Then, I compared its AIC values to find what distribution best fits the data:
> glm_1 <- glm(data$Y ~ data$X, family = "gaussian")

> glm_2 <- glm(data$Y ~ data$X, family = "Gamma")

> AIC(glm_1, glm_2)
       df      AIC
glm_1  3       38.52803
glm_2  3       26.09031

As you can see, the Gamma distribution better explains my data distribution. However, when performing a Shapiro-Wilk test, I got an indication that there is normality in Y:
> shapiro.test(data$Y)

    Shapiro-Wilk normality test

data: data$Y
W = 0.88314, p-value = 0.1417

When viewing the distribution of Y on a histogram and on a Q-Q plot, I couldn't convince myself that the distribution is normal:
hist(data$Y, main = " ", xlab = "X")


qqnorm(data$Y)
qqline(data$Y)


According to the results above, would it be incorrect if I performed a simple linear regression of Y ~ X, or would it be more appropriate to make a GLM considering the Gamma distribution?

*

*Remember that failing to reject the null hypothesis is NOT accepting the null hypothesis. The Shapiro-Wilk test is a test of the null hypothesis of normality. The power of a test is in part a function of sample size. Therefore, if you have few observations, you're unlikely to be able to reject the null (although there are some caveats to this I won't go into) as the test is not well-powered to detect small departures from normality. The main takeaway being the test is not very useful when you have very little data.


*However, the test is ALSO not very useful when you have lots of data! Imagine you have infinite data but you distribution is only a tiny bit different from what you'd expect of a normal distribution. The test will give you a p-value pretty much at 0 indicating non-normality. But of course, that's not very useful -- your test told you it's not normal, but what you really care about is just how not normal is it. The test won't get you there. This is a general feature of "testing" assumptions using point nulls.


*You seem to think the data being normal is the Gauss-Markov assumption. That is decidedly not the case! See the first two paragraphs: https://en.wikipedia.org/wiki/Gauss%E2%80%93Markov_theorem


*Even if it was, the assumption is not about the data -- it's about the error term. Note that the error term is NOT the model residual. The error term is something you do not observe. The residual is specific to your model.


*The Gamma distribution might better fit your data, but model fit statistics can be not-so-great. Ideally you'd be able to test on out-of-sample prediction with some more data not used to train your model.


*More generally, the Gamma distribution is bounded from below by 0: i.e. it only has support on $(0,\inf)$. If you data can EVER be negative, then it would be quite silly to use that likelihood since you would be using a model that cannot possibly produce the population data, even if it fits your sample well enough.
