# In R, find function F(x) to transform values in a vector to a normal distribution?

I have a PDF (Probability Density Function) generated from a vector of 1,000,000 empirical values. This empirical PDF is heavily skewed to the right.

In this form, I can't make accurate predictions using a linear regression.

To fix this, is there some method to find the function F(x) to transform (i.e. "squash") the values in the vector into a standard normal distribution, so I can feed said transformed vector into a linear regression?

Of course, this would also involve finding the inverse of F(x) that transforms (i.e. "de-squashes") any predictions back into the original empirical PDF.

What I have tried

So far, I have managed to generate the density function from the empirical data: Here is the R code:

par(mfrow=c(2,1))

install.packages("bootstrap")
library(bootstrap)
data(stamp)
nobs <- dim(stamp)
hist(stamp$Thickness,col="grey",breaks=100,freq=F) dens <- density(stamp$Thickness)
lines(dens,col="blue",lwd=3)

plot(density(stamp$Thickness),col="black",lwd=3, main="Simulation to choose density plot") for(i in 1:10) { newThick <- rnorm(nobs,mean=stamp$Thickness,sd=dens$bw*1.5) lines(density(newThick,bw=dens$bw),col="grey",lwd=3)
}

# If I wanted to do a linear regression to predict stamp thickness,
# what is the function F(x) to "squash" (i.e. transform) the "stamp"
# vector into a normal distribution, and the corresponding inverse
# function Finv(x) to "desquash" (i.e. untransform) any predictions back
# into the original prediction?


Update 1

@Andre Silva sugggested that:

What need to have normal distribution are the residuals (predicted versus observed) derived from your (multiple) linear regression model.

According to post on Multiple Linear Regression:

After fitting the regression line, it is important to investigate the residuals to determine whether or not they appear to fit the assumption of a normal distribution. A normal quantile plot of the standardized residuals y - is shown to the left. Despite two large values which may be outliers in the data, the residuals do not seem to deviate from a random sample from a normal distribution in any systematic manner. Update 2

See Left skewed vs. symmetric distribution observed for R code that illustrates that the only relevant concern is if the residuals are normally distributed.

• Look at a plot of the natural log of the data. I bet that looks a lot closer to normal! – wcampbell Apr 7 '13 at 15:44
• @wcampbell Thanks for suggesting that I use natural log, this is a good first step. The example that I gave was just an example, the actual problem has a PDF that is so skewed that a log transform won't work. – Contango Apr 7 '13 at 16:50

@AndreSilva is right that regression does not require the data to be normal. The assumption of linear regression is that the residuals are normal. It may help to read this thread: what-if-residuals-are-normally-distributed-but-y-is-not, to clarify this issue.

However, this point does not go far enough. First, the normality of the residuals serves to ensure that you can trust the standard p-values that software outputs. Even then, you can trust your p-values without normal residuals, if your N is sufficiently large. With $N>1,000,000$, there is likely to be little reason to worry about the validity of your p-values.

At any rate, if you want to make a predictive model, whether or not the residuals are normal is irrelevant. OLS regression methods are unbiased whether the residuals are normal or not, without regard for N. Thus, if you want to make point predictions (i.e., $\hat y_i$, the predicted mean of the conditional response distribution where $X=x_i$), you will be fine. If you want to make interval predictions, you can do that as well, you just shouldn't use the normal distribution (which would be the default) to do so. Instead, you can use the estimated density of your residuals to make prediction intervals.

• Thanks, @StéphaneLaurent. I still wonder if we're miscommunicating. The fact that the model is $Y\sim\mathcal N(X\beta,\sigma^2I)$ (& assuming the residuals are normal) is insufficient for $Y$ to be normal (of course, it could be though), it would depend on the distribution of $X$. To make sure we're on the same page, would you mind taking a look at my answer at the thread I link above? There you can see that $Y$ is very clearly non-normal (tri-modal, in fact), but the residuals are normal & everything is fine. That's what I'm getting at. Is that what you disagree with? – gung - Reinstate Monica Apr 8 '13 at 12:50
• @gung Each $y_i$ has a normal distribution, but with a mean $\mu_i$ depending on $i$. So the data are normal, but they are not identically distributed. You are claiming that $\mu + \epsilon$ has not a normal distribution when $\mu$ is a fixed number and $\epsilon \sim {\cal N}(0, \sigma^2)$ ! – Stéphane Laurent Apr 8 '13 at 14:36
• @StéphaneLaurent, I wondered if this might be a semantic issue. I now gather that by "the data are normal", you mean '$y_i$ is normal conditional on $x_i$', which is almost identical to what I mean by "residuals" (the difference being that the residuals have had the $X\beta$ subtracted off). I'm fine w/ that, of course. The problem for the OP, if you look at the top histogram (& for many people, in my experience), is that they interpret "the data are normal", as meaning 'the marginal distribution of $Y$ is normal'. Thus the OP is concerned about that top histogram, .../... – gung - Reinstate Monica Apr 8 '13 at 15:27
• .../... despite the fact that the residuals are normal enough, which you can see from the qq-plot of the residuals at the bottom (& hence $y_i$ would be normally distributed conditional on $x_i$). I think what's important here is to get people clear on the fact that the marginal distribution of $Y$ is irrelevant. That's why I use the term "residuals", instead of 'the data', which I think is ambiguous, & have found to lead to misinterpretations (an example of which is above). – gung - Reinstate Monica Apr 8 '13 at 15:30
• @gung I have elaborated my point in an answer. The $x_i$ are not random in a linear model, or they are random but the model is conditional to the $x_i$, hence the $y_i$ are normal. – Stéphane Laurent Apr 8 '13 at 16:12

Your data do not need to be "Normal". What need to have normal distribution are the residuals (predicted versus observed) derived from your (multiple) linear regression model. Did you test this presupposition?

• Thanks! I think this may be a clue to solve this. I'm trying to find the code to do a multiple linear regression model on the components returned by a call to density() in R. – Contango Apr 7 '13 at 16:49

A classical linear model (such as simple linear regression) for a sample $y=(y_1, \ldots, y_n)$ has form $y_i = \mu_i + \epsilon_i$ where the "error terms" $\epsilon_1, \ldots, \epsilon_n \sim_{i.i.d} {\cal N}(0, \sigma^2)$ and the means $(\mu_1, \ldots, \mu_n)$ are unknown numbers satisfying some linear constraints.

Thus each $y_i$ is assumed to be generated from a normal distribution: $y_i \sim {\cal N}(\mu_i, \sigma^2)$ but drawing a histogram or an estimated density for the $y_i$ does not allow to check normality because they are not identically distributed (the distribution of $y_i$ depends on $i$ through $\mu_i$). In other words the $y_i$ are assumed to be generated from a normal distribution but not from a common distribution. If you have an "i.i.d" sample from an unknown distribution then you can estimate this distribution with a histogram or an estimated density, but if the sample is not "i.i.d" the histogram or the estimated density is useless.

However the error terms $\epsilon_i$ are identically distributed. One never knows the realizations of the $\epsilon_i$ because the $\mu_i$ are unknown, but the residuals $\hat\epsilon_i = y_i - \hat\mu_i$ approximate the realizations of the error terms and one can assess departure from the normality on the sample of residuals.

In some cases, for instance a one-way ANOVA, you can separately check normality by drawing a histogram of the $y_i$ in each group of individuals defined by the factor because the model assumes the $y_i$ are normal and i.i.d in each group. If the group sizes are small it is better to check normality with the residuals.

• Thank you, this makes it a lot clearer. For those of us who are just beginning to learn stats, IID = Independent and Identically Distributed. – Contango Apr 9 '13 at 21:05