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

Question

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)[1]
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.

Answer 1

@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.

Answer 2

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?

Answer 3

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.