# Heteroscedasticity problem violating assumptions for lm and glm

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()

• Welcome to Cross Validated! family=Gaussian means that the conditional distribution is Gaussian.
– Dave
Jun 27, 2023 at 21:11
• It's complicated. A linear model does not have to be estimated using least squares, and this question of mine asks about linear models with binomial likelihoods. Likewise, using a Gaussian likelihood does not require a linear model, and this answer of mine at least kind of alludes to this. However, in code, family=Gaussian in the glm function gives the same results as lm, so an assumption of a linear model with $iid$ Gaussian errors (same assumptions as lm has for standard errors).
– Dave
Jun 27, 2023 at 21:43
• It may be helpful if you provided scatterplots of your data or the residual plots you refer to. Jun 28, 2023 at 0:09
• I’d like to hear from others more experienced than me, but I don’t think the QQ plot suggests a worrisome violation of normality, also considering that LM is robust to mild violation to normality when the sample size is relatively large.
– kk68
Jun 28, 2023 at 19:38
• Those plots all look fine to me, better than most I have seen in my work to be honest. Jun 28, 2023 at 19:40

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.

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 :).

• And, dont use rank tests, unless you have a very good argument. They don't help you to interpretate the result, and they compare the order of the data, not the actual values of the data. It is often said that they compare medians - this is not true! you can generate two samples with same median but with clearly significant difference in rank test. Jun 28, 2023 at 9:49
• Thanks! "simply just ignore this" is also what I've heard from other people that have more knowledge within stats. Especially when transformations does not work (I've tried both log and square root transformation) and that clearly messes up the normal distribution but solves the heteroscedasticity.. I've been thinking about using white.adjust as well because of the heteroscedasicity but I still have the problem with the normal distribution. I've heard people say that if you have a lot of observations you can basically assume normal distribution even thoughshapiro wilk test says otherwise. Jun 28, 2023 at 19:16
• I also tried mixed model with random effect but to no success.. My "limited" knowledge of stats takes me to the non-parametric tests. You say I should avoid them. Although Ive been looking into Aligned Ranks Transformation ANOVA which is a non-parametric test where you can test interaction effects. Hear me out, I know you said avoid the non-parametric tests. But could you use them to justify fitting a linear model? that is, if i get similar results in the non-parametric i go for the lm/glm because its more common in the field, but i also used non-parametric to avoid type 1 error? Jun 28, 2023 at 19:19
• @Kirsten There are some who would argue to use rank tests unless there is a very good reason not to.
– Dave
Jun 28, 2023 at 20:11
• I'm not sure employing splines would really be that necessary for this model. The residuals would likely look a lot more bizarre if there was some severe departures from linearity. Jun 28, 2023 at 23:07

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.