# Could one say that one dataset is distributed more normally than another?

I am trying to fit a simple linear model: experiment ~ calculated_1. From the basic model, I get residuals_1. And I know that another calculated data set can theoretically improve prediction of the experimental data. So I create the other model: experiment ~ I(calculated_1 + calculated_2). From this I get residuals_2.

residuals_1 <- c(-0.6512352, -0.7809935, 0.4313386, -0.1061949, -0.5986882, -0.8560606,
0.4278563, 0.6283887, -1.0142389, -0.8876992, 0.1680664, 0.2761103,
0.6310700, 0.2358954, -0.1694672, -0.1466784, 0.2525157, 1.6015804,
-0.1060606, 1.2452776, -0.8153085, 0.2490334, -0.0145075)

residuals_2 <- c(-0.684552455, -1.020307239, 0.868944787, 0.191034959, -0.318848920,
-0.961003963, -0.687738963, 1.460222114, -1.092411033, -0.480501852,
-0.043288749, -0.231946550, 0.504105347, 0.753821856, -0.036888490,
0.226492631, 0.197110863, 0.350608105, -0.156642627, 2.078234441,
-1.508965041 ,0.004402459, 0.588118321)
qqnorm(residuals_1)


qqnorm(residuals_2)


As one can notice the variance of residuals_2 is slightly greater then the variance of residual_1. But is it correct to say that adding calculated_2 term make the residuals to be more likely normally distributed? And provide a better description of the experimental data?

• Why do you think more normal residuals improve the predictive strength of a model? Aug 3, 2015 at 20:35

It is possible for one set of data to more closely approximate a given theoretical distribution than another. One way to think about how close one distribution is to another is to use the Kullback-Leibler divergence. This is a measure of how much information is lost by using one distribution instead of the other. More specifically, we can determine the information loss associated with using the densities from a true normal fitted to your data instead of a more flexible kernel density fitted to your data.

library(flexmix)
sim = function(x, ref){     # find similar values
lr   = length(ref)
indx = vector(length=lr)
for(i in 1:lr){
indx[i] = which(abs(x-ref[i])==min(abs(x-ref[i])))
}
return(indx)
}
get.kl = function(res){     # get densities & K-L divergence
d.k = density(res)
d.k = d.k$y[sim(d.k$x, res)]
d.n = dnorm(res, mean=mean(res), sd=sd(res))
return(KLdiv(cbind(d.k, d.n))[1,2])
}

get.kl(residuals_1)  # [1] 0.01309022
get.kl(residuals_2)  # [1] 0.005502055


More than twice as much information would be lost when approximating residuals_1 by a normal than approximating residuals_2.

It is also possible for a misspecified model to yield non-normal residuals:

set.seed(3)
x = runif(30, min=0, max=10)
g = rep(0:1, each=15)
y = 17 - 0.3*x + 2*g + rnorm(30)

m1 = lm(y~x)
m2 = lm(y~x+g)
get.kl(resid(m1))  # [1] 0.04971849
get.kl(resid(m2))  # [1] 0.008378342


Having said those things, using how closely your residuals approximate a normal distribution to select between models with different covariates would be rather odd. It is much more typical for people to select between models using a hypothesis test (if the models are nested), the AIC (if they are not nested), out of sample predictive error via cross-validation (if prediction is the ultimate goal of the model), etc. People do (and should) evaluate the residuals of their models, but they might use that to determine if a transformation is necessary, or if they should bootstrap instead of rely on the standard errors, for example.

• Great answer. And nice bit at the end. Could you expand a bit more on this point, as I am not sure what it is saying. "or if they should bootstrap instead of rely on the standard errors". Bootstrap to do what? Aug 4, 2015 at 5:28
• Thank you very much for so complete and helpful answer! Predictive perfomance evaluation is not the main goal o this study but I shall do it. THanks for the hint. Am I right that according to K-L divergence information lost in case of residuals_2 is lower in comparison with residuals_1? Aug 4, 2015 at 10:46
• @MHH, the standard errors that come with your output assume the sampling distribution is distributed as a normal or t. If the sampling distribution isn't, then your inferences would be invalid. You can bootstrap to get an empirical estimate of the sampling distribution, which you can use for inference instead. Aug 4, 2015 at 11:44
• @dkarlov, yes, according to K-L divergence information lost in case of residuals_2 is lower in comparison with residuals_1. Aug 4, 2015 at 11:45