Heteroscedasticity problem violating assumptions for lm and glm

Question

I’ve been trying to fit some data for a manuscript to a model for the past two days and I keep running into problems with violation of assumptions for the lm and glm functions in R and it’s driving me crazy.

The dependent variable is continuous so I thought I could plot the data in a linear model in R (function lm()) with an interaction effect between the two independent variables. However, when doing diagnostics my data is non-normally distributed and heteroscedastic. One solution to this would be transformation, when I log transform the dependent variable I get equal variances (homoscedacity) but instead my residuals get even more non-normally distributed when looking at both diagnostic plots and running a Shapiro-Wilks test. That is, they look even worse than before. I think this is very odd, especially since all my data is > 1.

I then thought of fitting the data to a GLM instead. This should work since the data doesn’t have to be normally distributed (even with family=gaussian if I’m not mistaken). However, the problem with heteroscedasticity is still present and if I use a log transform that influences the interaction effect (this is what I see and I’ve heard this before but can’t recall the source) and I’m interested in the interaction effect. I’ve been looking a bit on non-parametric tests and Kruskal Wallis tests seems to be able to look at interactions with the ranks() in R but I don’t have any experience with this.

I’m really not a stats guru which makes it even harder to navigate how to proceed. I know it might be hard to give direct answers but I was mostly wondering whether:

1. Someone has experienced a similar situation where they had the statistics thought out and when the data was collected they realized that, due to assumptions being violated, they couldn’t stick to the planned procedure.

2. Someone has a solution for my problem, or valuable input on how to proceed. E.g. non-parametric ways in R with interaction effects or a solution to the violations of assumptions.

EDIT:

I added output from diagnostic figures with this code:

plot(model)

ggplot(aa, aes(x = interaction(A,B), y = .resid)) + geom_boxplot()

The car::leveneTest(.resid ~ A*B, data = dat) results show the following:

P = Levene's Test for Homogeneity of Variance (center = median)
       Df F value   Pr(>F)    
group   8   3.647 **0.000515** **

And using shapiro.test(modartificialglm$residuals), I get this:

non-normal distribution - p-value = **0.01105**

shapiro-Wilk normality test

data:  model$residuals
W = 0.98432, p-value = 0.01105

I should also mention that the model looks like this:

model<-glm(y~ A* B,  dat, family="gaussian")

Second to last model (blue circles) are plot(glsmod)

GLS model : glsmod <- gls(y~ A* B, dat, weights = varIdent(form = ~ 1 | A* B) all the other diagnostics are done on glm()

Answer 1

As the other answers say, the unweighted lm() model (equivalent to your glm() with a Gaussian family; in the form you used, glm() and lm() are based on the same assumptions) wasn't that bad. Even with the complication of potential heteroscedasticity, the Normal Q-Q plot wasn't very bad at all. You have a reasonably large number of observations, so it's possible to have "statistically significant" violations of normality and heteroscedasticity that aren't of practical importance. See, for example, the discussion of whether normality testing is essentially useless.

In a later version of the question you added some results from a generalized least squares (gls()) model, in this case a weighted least squares model. That's a well accepted way to deal with heteroscedasticity. You allowed for different variances among all combinations of A and B, with observations weighted inversely to their variance estimates (weights = varIdent(form = ~ 1 | A* B)). If you are worried about the heteroscedasticity in the unweighted models, that would be a good way to go.

The terminology can be very confusing when you're first starting to learn these methods. In particular, "generalized least squares" sounds a lot like "generalized linear model" even though they can have quite different applications. Further complicating matters, either of those could implement a "generalized additive model" to fit continuous predictors flexibly. Make sure you know what type of "generalized" model you're dealing with.

Answer 2

Based on your first comment about continuous outcomes, a Gaussian model seems the right choice. There are of course some situations where it may be better to dichotomize your outcome, say if you have blood test results that have a lower or upper detection boundary.

Leaving that for now, lm is basically a special case of glm. Yes there is a difference in their estimation method, so you may get slightly different results but for all practical applications it should not change your results. For both of them, you will have an underlying assumption about normally distributed errors (that is mistakes in measurements or natural variation) - thus the considerations about residuals.

Going back to the linear model, often written in R as y~x+c, where y is our outcome x (an explanatory variable) and c covariates. If the model does not fulfill the criterion, non-normally distributed residuals and heteroscedasticity, then something went wrong when we specified our model. It simply does not reflect reality.

A simple solution is to ignore this - This is obviously not a first choice, but we are allowed to say: In our opinion this is the best possible model, it does not fit perfect but we think it is acceptable for our purpose. In this argument it is nice to know, that the linear model is robust even under reasonable misspecification (unfortunately you have to find the reference yourself, I don't remember the paper).

Before we accept a wrong model - as you suggest yourself, transformation is a good start. It should be noted here, that this will make the model result more difficult to interpret, so if we are talking small misspecifications, I would often prefer the misspecification. You log-transformed your data but another option is to square it. Basically there is no limit for the choice.

Another question you should ask yourself is - did I specify the correct model? Did I leave out any covariates that may explain this? Or have I mis-specified the variance structure. I don't know how familiar you are with regression models, but there is a model type called mixed models. These allow you to specify the variance structure of you data, and are often used for repeated observations. I will not explain it here but googling "introduction to linear mixed model" will give you many hits with easy explanations.

Closer to transformation is transforming covariates, creating functions of them - fx splines for continuous variables such as age. I will recommend that you do not us advanced transformations for variables you aim to interpret. This will make your life too hard. If you know you have interaction terms, I would actually start by considering two or more separate models. This will simplify you presentation of the data - and even if it does not solve your problem it may help you to identify it.

Hope this helps you to get a little further. And remember that model fit is important, but a model where you can identify and explain problems may be worth more.

I would also recommend you look up bootstrapping. If you have a mis-specified model variation, this method will help you to create more "correct" confidence intervals. Its strength over the robust sandwich estimator is that bootstrap relies purely on sample distribution -> less assumptions :).

Answer 3

Comments

Thank you for providing the plots. I am in agreement with Kirsten and Bowman that your residual and QQ plots actually look a lot better than I anticipated. Therefore I won't add anything there.

I saw in additional comments that you were considering other models. Stick to lm. As others have already noted, coding glm with the Gaussian family does pretty much the same thing as lm. Additionally, running a mixed model is only advisable if you actually need one (McNeish et al., 2016). You haven't mentioned any reason for why you would even need that or what random effects are present in your data. You can certainly comment below if there are indeed these influences in your data. Though Kirsten mentioned the use of splines, I don't see any reason why you would need them added to your model, but that would of course be dependent on whether or not your data is non-linear. I would imagine your residuals would look way worse if that was an actual issue. However, if you feel you would like to learn about them for your own sanity, I have included a reference below that covers the basics.

Reference

• McNeish, D., Stapleton, L. M., & Silverman, R. D. (2016). On the unnecessary ubiquity of hierarchical linear modeling. Psychological Methods, 22(1), 114–140. https://doi.org/10.1037/met0000078
• Perperoglou, A., Sauerbrei, W., Abrahamowicz, M., & Schmid, M. (2019). A review of spline function procedures in R. BMC Medical Research Methodology, 19(1), 46. https://doi.org/10.1186/s12874-019-0666-3