# Is there a difference between a GLM with gaussian family and a standard linear model?

I am trying to model some data and I'm going around in circles trying to figure out what family I should be using when writing the generalised linear model.

Here is some example data:

df <- data.frame(x = c(12.349731, 11.171157, 5.690359, 12.291016, 13.510602, 12.563224, 9.630021, 9.919312, 9.962405, 7.148739,
8.964152, 11.024807, 9.114435, 9.683293, 9.211215, 7.799241),
y = c(-24.620238, -11.496512, -36.912983, -2.884713, 3.369998, -13.716537, -8.734567, -8.707669, -8.595272, -32.649879,
-30.428423, -20.543413, -16.415703, -10.638254, -28.197134, -33.543908))


Ploting x against y shows an apparently linear relationship, however the variance increases with the mean and the error distribution is decidedly not normal. This suggests that I need to use a generalised linear model.

However, how do I choose which family? I have seen this question asked here and I understand that it is a bit of a grey area.

I know I can rule out a number of them as this isn't binomial or poisson data, but I also do not think that inverse.gaussian or Gamma is correct because the data can clearly be negative.

So this leaves gaussian, but is this not just the same as a standard linear model with all the assumptuions that come along with that?

One way I was thinking of was to change the link function from identity to something else, but I am not sure how to choose a different link function aside from running through them all and seeing which model fits 'best'.

I would appreciate any insights on where to go from here.

• I don't see much evidence that the variance increases with the mean, nor that the error distribution is non-normal. You might be reading too much into chance variation. The principal message of the diagnostic plots is to highlight point 1 as being especially interesting: it's pretty far from all the other points. However, falling in the middle of the explanatory values, it has no leverage and so you don't need to worry about it unless your objective is to obtain an accurate estimate of the error variance. – whuber Jul 16 at 13:10
• Thank you for your reply. So the patterns I'm seeing aren't outside of what you might expect from random variation? Is that because the sameple size is relatively small? – tom91 Jul 16 at 13:52

I don't see much evidence that the variance increases with the mean, nor that the error distribution is non-normal. You might be reading too much into chance variation. The principal message of the diagnostic plots is to highlight point 1 as being especially interesting: it's pretty far from all the other points. However, falling in the middle of the explanatory values, it has no leverage and so you don't need to worry about it unless your objective is to obtain an accurate estimate of the error variance.

A good way to check is to simulate similar datasets and study the variations among the diagnostic plots. That's a bit of an effort but IMHO really pays off.

Among the many ways to create a simulation, parametric and nonparametric bootstrapping are appealing. Both of them add random terms to the fitted values to create new datasets. The parametric method generates iid Normal residuals of the size expected from the residual standard deviation of the fit. The nonparametric method samples (with replacement) from the original residuals. The parametric method gives you a sense of what results look like when variations are generated from a Normal distribution. The nonparametric method gives you a sense of what results look like when the variations come from the distribution of residuals you observed: this can give you a sense of whether non-normality of those residuals may be affecting your interpretation.

Here, for instance, are Normal QQ plots for 23 bootstrapped datasets using the parametric method. For comparison, the QQ plot for the original data appears at the upper left in red.

Let's examine the QQ plots, for instance.

The deviation from the diagonal line at the lower left in the original data may have caught your attention, as it should. The question is whether to assume it's something you should treat as an inherent property of the process or population you are studying, or whether to treat it as "random variation." Comparing this deviation to deviations apparent in other plots (as in the top (first) row, fourth column from the left, or in several along the bottom row) suggests random variation would be the better judgment.

You may perform similar simulation-based analyses to see whether you agree with my initial conclusions.

This R program generates comparable tableaux of diagnostic plots (and scatterplots, corresponding to a plot.type of 0) for any x,y data frame, using either parametric or nonparametric bootstrapping, and thus may be generally useful.

fit <- lm(y ~ x, df)
df$$r <- residuals(fit) df$$y.hat <- predict(fit)
df$$y.boot <- df$$y
sigma <- summary(fit)$sigma boot.type <- c("Parametric", "Nonparametric")[2] plot.type <- 2 # Usual plots are 1 (Res vs. fit), 2 (Q-Q), 3 (Scale-Loc), 5 (Res vs Leverage) n.iter <- 24 ncol <- floor(sqrt(n.iter * 2/(-1+sqrt(5)))) nrow <- ceiling(n.iter/ncol) par(mfcol=c(nrow, ncol)) set.seed(17) for (i in 1:n.iter) { if (i != 1) { if(boot.type == "Parametric") { df$$y.boot <- df$$y.hat + rnorm(nrow(df), 0, sigma) } else { df$$y.boot <- df$$y.hat + sample(df$r, nrow(df), replace=TRUE)
}
}
fit.boot <- lm(y.boot ~ x, df)
if(plot.type==0) {
with(df, plot(x,y.boot, main=paste(boot.type, "Sample", i)))
abline(fit.boot)
} else {
plot(fit.boot, which=plot.type, col=(i==1)+1)
}
}
par(mfcol=c(1,1))

• Thankyou for this indepth and well thought out response. You're right it looks like I had been reading too much into the diagnostic plots in the original model. I haven't done much in the way of bootstrapping before but I will look into it. Thanks again! – tom91 Jul 16 at 14:41