# Can I violate assumptions of normality for categorical linear regression models?

I'm using packages included in this R/rStudio tutorial to set up some linear regression models comparing a continuous dependent variable (eccentricity) to three categorical variables (year, bird tissue type, and their interaction). The endgame is to rank these by AIC. Everything (including the models) runs fine in R. However, the tutorial indicates that for categorical linear regression, I should check assumptions of normality before running the linear models (explanation here), and every category of my data but one (not shown) fails the Shapiro-Wilk normality test. I'm reading elsewhere (e.g. here) that with a sample size higher than thirty (each group has a minimum of 100 points) central limit theorem should obviate the need for a normality test, and that normality tests become oversensitive with high sample size. That being said, some of my data is obviously skewed, so it's more a question of whether I can run the same analysis with skewed data than a question of whether the underlying distribution is normal. Can I proceed with linear models for these groups, or should I be using a different statistical test that doesn't rely on assumptions of normality being met?

Here are the density histograms, qq-plots, and Shapiro-Wilks p-values for some of my groups:

• Are your data bounded by 0 and 1? Can eccentricity be less than 0 or greater than 1? Feb 13 at 0:34
• Eccentricity should be bound by 0 and 1. Feb 13 at 0:37

I think we need to back up.

Questions of normality are -- in my opinion -- the least important when it comes to linear regression. That being said, your model and data are showing other signs that something is wrong.

Your QQ plots are indicating that your model can not accomodate the boundedness of eccentricity. See how in the top right QQplot the residuals seem to top out? This is telling me linear regression is not appropriate.

How is eccentricity measured? It might be preferable to use a logistic or beta regression rather than linear regression.

EDIT:

OK, we have determined that the eccentricity is bounded between 0 and 1 and that it isn't the result of some discrete process (e.g. successes divided by trials).

I think a beta regression might be more appropriate. Because I don't have access to your data, I will simulate something similar, shown below.

If this data is in a data frame d, then we can call the beta regression as follows

# you might need to install this
library(betareg)

fit <- betareg(y ~ grp, data=d)


To test your hypothesis, you can use a likelihood ratio test. Unfortunately, I don't think R implements this for the betareg class, but we can program it fairly easily.

library(betareg)

fit <- betareg(y ~ grp, data=d)
# Compare against a null model
fit0 <- betareg(y ~ 1, data=d)

loglikstat <- as.numeric(2*(logLik(fit) - logLik(fit0)))

p.value <- pchisq(loglikstat, df=3, lower.tail = FALSE)
p.value
>[1] 0 # The p value is very small


Here is the full code to replicate my analysis

library(tidyverse)
library(betareg)

n <- 2000
mu <- c(0.9, 0.9, 0.5, 0.5)
k <- c(10, 10, 10, 10)

a <- mu*k
b <- (1-mu)*k
grp <- sample(1:4, size = n, replace=TRUE)

y <- rbeta(n, a[grp], b[grp])

d <- tibble(grp, y)

d %>%
ggplot(aes(y)) +
geom_histogram() +
facet_wrap(~ grp, ncol=1)

fit <- betareg(y ~ grp, data=d)
# Compare against a null model
fit0 <- betareg(y ~ 1, data=d)

loglikstat <- as.numeric(2*(logLik(fit) - logLik(fit0)))

p.value <- pchisq(loglikstat, df=3, lower.tail = F)
p.value

• My raw data is an xy scatterplot of carbon (x) and nitrogen (y) values that are continuous, can be positive or negative, and are in theory normally distributed. I used the siberR package to estimate bayesian ellipses for every combination of the above categorical variable using 95% of my raw data, and the eccentricity included here is the eccentricity of those ellipses (calculated by dividing the distance from the centroid to the foci by the length of the semi-major axis e=c/a). I am trying to see if the differences in eccentricity for each category are driven by one variable in particular. Feb 13 at 0:53
• @ElizaBeso000 From what you describe, I think a beta regression might be more appropriate. See my Edit Feb 13 at 1:11
• Thank you! This is so helpful! I just want to check - when specifying the df for pchisq, are you using the df from the alternative hypothesis model? My models vary between 3 and 7 dfs so I want to make sure I'm modifying the pchisq df appropriately for each comparison. Feb 13 at 2:28
• The degrees of freedom (df) should be the difference between the number of parameters in each model. I think the package l test has a function called lrtest that calculates the degrees of freedom automatically for me so you don’t make any errors. Feb 13 at 2:33
• So coming here with a follow up question - I am now trying to use a similar approach on the areas of the ellipses, which are obviously all positive, with a gamma regression based on your feedback above (i.e. glm(y ~x, data = data, family = Gamma (link = "log"). The above lrtest() code produces an output even though the summary for the gamma regressions and gamma null model don't show a log-likelihood value like the beta models do. Can I trust the output of lrtest() for a gamma regression or is there another approach I should use to produce a comparable output? Feb 16 at 6:11