# Non-normally distributed residuals for multivariate linear regression.Still a valid model?

I try to know if the independent variables are affecting the outcome of the dependent variable, but while the Shapiro-Wilk test shows residuals non-normally distributed, the autocorrelation of errors and heteroscedasticity are non-significant.

Sample size of 500. Dependent variable: mean = 4.8; SD=86.54; Range= 0.0 - 2000. 3 independent variables, all of them continuous.

I start with a general linear regression (glm) with the IV log-transformed:

> glm_data <- glm (log(DV+1) ~ IV1 + IV2 + IV3, data=data.raw)
> summary (glm_data)

Call:
glm(formula = log(DV + 1) ~ IV1 + IV2 + IV3, data = data.raw)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.2791  -0.1750  -0.1088  -0.0368   7.4115

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.841757   0.571616  -1.473   0.1415
IV1          0.011198   0.010277   1.090   0.2764
IV2          0.014509   0.005877   2.469   0.0139 *
IV3          0.073306   0.044355   1.653   0.0990 .
---

(Dispersion parameter for gaussian family taken to be 0.3914142)

Null deviance: 212.28  on 536  degrees of freedom
Residual deviance: 208.62  on 533  degrees of freedom
(20 observations deleted due to missingness)
AIC: 1026.2

Number of Fisher Scoring iterations: 2


I check the distribution of the residuals

> shapiro.test (glm_data$residuals) Shapiro-Wilk normality test data: glm_data$residuals
W = 0.2666, p-value < 2.2e-16


Transformed DV or not, the value is the same.

The density of the residuals is just one peak very close to the Y axis.

I check the autocorrelation of errors using Durbin-Watson test (durbinWatsonTest {car})

> durbinWatsonTest (glm_data)

lag Autocorrelation D-W Statistic p-value
1       0.1091337      1.780399   0.066
Alternative hypothesis: rho != 0


The p-value is always over 0.05, but not by much. It's over 0.6 if I don't log-transform the DV.

I check the heteroscedasticity using Breusch-Pagan test (bptest{lmtest})

> bptest (glm_data)

studentized Breusch-Pagan test

data:  glm_data
BP = 3.8858, df = 3, p-value = 0.2741


And correlation in my independent variables with vif{car}:

> vif (glm_data)

IV1      IV2      IV3
6.078988 1.607718 5.236179


I feel they could be lower, but I'll go with them.

I don't know how I can keep going, but I bootstrapping (Boot{car}) and calculate the confident intervals for the glm mode:

> boot_glm_data <- Boot (glm_data, R=2000)
> summary (boot.glm.data)
R  original  bootBias   bootSE   bootMed
(Intercept) 1000 -39.27793 -1.581097 36.70436 -36.08043
IV1         1000   0.84846  0.036708  0.81795   0.79462
IV2         1000   1.33761  0.049331  1.21655   1.31777
IV3         1000   0.93736  0.035413  1.18028   0.87289

> confint (boot.glm.data)
Bootstrap quantiles, type =  bca

2.5 %   97.5 %
(Intercept) -214.52868331 1.858335
IV1           -0.05086198 4.932204
IV2            0.06636597 7.352573
IV3           -1.01565092 3.731557
Warning message:
In norm.inter(t, adj.alpha) : extreme order statistics used as endpoints


I don't know what else I can do. Should I ignore the Shapiro-Wilk normality test? Can I use a different regression or test?

• Normality, autocorrelation and heteroskedasticity need not all go hand in hand; hence, your title does not indicate a paradox or a problem. A different title might be more informative. As I understand, you are trying to build a good model so that the model residuals would satisfy the model assumptions. You might want to reshape or refocus the question slightly to get what you really want. By the way, is normality an important assumption for this class of models? For example, in linear regression normality is generally not important, especially in large samples due to the central limit theorem. – Richard Hardy Jun 4 '15 at 18:26
• Thanks @RichardHardy. Maybe the question was not the most appropriate. I know that those three characteristics do not need to go together, but I thought the three should be satisfied to validate the regression model. I'll think about a better question, but yes, I'm trying to build a model with the 3 IV as predictors of the 1 DV. I heard and read about "large samples", but I don't know when a sample could be considered "large" enough to not care about normality of the residuals. – Saccharo Jun 4 '15 at 18:33
• Sample size of 500 observations seems large enough if you only have 3 regressors. For 3 regressors, even 100 (if not less) would be OK, I think. – Richard Hardy Jun 4 '15 at 18:53
• Normality is generally only important for hypothesis tests and intervals - and even then only if you're using normal-theory procedures - and even then, the assumption is often not so critical at large sample sizes (the major exception being prediction intervals for individual values, for which normality can matter even with large samples). With 500 points, you're very likely to reject normality even when the distribution is close to normal ... and yet with so many points, it may not actually matter terribly much. ... (ctd) – Glen_b Jun 4 '15 at 19:12
• (ctd) ... It's the extent of (and kind of) non-normality that matters and its degree of impact on your inference (i.e. "effect size"), rather than p-value in the normality test (which may be mostly telling you "your sample is large enough to detect even mild non-normality").. – Glen_b Jun 4 '15 at 19:13