A nice approach for checking the fit of your assumed model to the data, accounting for features, such as, over-dispersion, non-normality, zero-inflation is the simulated scaled residuals provided by the DHARMa package. If your assumed model is correct, these residuals should have a uniform distribution. You can find more details on the procedure they are defined and used in the vignette of the package.
As an example, using the simulated example above, I compare below the fit of the Gamma model to the fit of the wrong normal model:
set.seed(0)
N <- 250
x <- runif(N, -1, 1)
a <- 0.5
b <- 1.2
y_true <- exp(a + b * x)
shape <- 2
y <- rgamma(N, scale = y_true / shape, shape = shape)
gamma_model <- glm(y ~ x, family = Gamma(link = "log"))
normal_model <- glm(y ~ x, family = gaussian())
library("DHARMa")
check_gamma_model <- simulateResiduals(fittedModel = gamma_model, n = 500)
plot(check_gamma_model)
check_normal_model <- simulateResiduals(fittedModel = normal_model, n = 500)
plot(check_normal_model)
