# Stop using bivariate correlations for variable selection

20 Mar 2014# Stop using bivariate correlations for variable selection

Something I've never understood is the widespread calculation and reporting of univariate and bivariate statistics in applied work, especially when it comes to model selection. Bivariate statistics are, at best, useless for multi-variate model selection and, at worst, harmful. Since nearly all questions of interest are highly dimensional, the use of bivariate statistics is worrisome.

What typically happens is a researcher sits down with their statistical software of choice and they compute a correlation between their response variable and their collection of possible predictors. From here, they toss out potential predictors that either have low correlation or for which the correlation is not significant.

The concern here is that it is possible for the correlation between the marginal distributions of the response and a predictor to be almost zero or non-significant and for that predictor to be an important element in the data generating pathway.

Consider the following simulation that illustrates the difference between the
marginal relationship between the response and a predictor and the true multi-
variable relationship. The function `genY`

generates the response using user
provided values for \( \beta \) and adds in a little white noise.

```
genY <- function(betas, X, stdPercent = 0.25) {
modelMatrix <- cbind(1, X)
yhat <- modelMatrix %*% betas
y <- yhat + rnorm(nrow(yhat), 0, stdPercent * mean(yhat))
y
}
```

The function `cor`

returns the correlation between two values and works with
a vector and a matrix argument but `cor.test`

requires two vectors. For this
I wrote a simple function `sigCor`

that returns `TRUE`

if a correlation between
the vector `x`

and `y`

is significant at a user provided threshold and `FALSE`

otherwise. This function is designed to work with `apply`

to filter our
predictor set.

```
sigCor <- function(x, y, threshold = 0.2) {
p <- cor.test(x, y)$p.value
p <= threshold
}
```

The main function then generates a random covariance matrix and draws from a
multivariate normal centered at 0 with that covariance matrix. The values of
the matrix `X`

are then checked for correlation with the response `y`

. If the
absolute value of the correlation exceeds some threshold or if the p value for
the significance of the correlation is less than some threshold, the variable is
kept. Otherwise it is dropped. The function returns the number of retained
variables for each selection method.

```
simulateSelection <- function(nx, n, corLim = 0.2, pLim = 0.2) {
require(MASS)
require(clusterGeneration)
betas <- rep(1, nx + 1)
mu <- rep(0, nx)
sigma <- genPositiveDefMat(nx, covMethod = "unifcorrmat")$Sigma
X <- mvrnorm(n, mu, sigma)
y <- genY(betas, X)
xcor <- cor(y, X)
limitedX <- X[, abs(xcor) >= corLim]
pX <- X[, apply(X, 2, sigCor, y, pLim)]
nCoef <- c(ncol(X) + 1, ncol(limitedX) + 1, ncol(pX) + 1)
}
```

If we keep the defaults of `corLim = 0.2`

and `pLim = 0.2`

(keep any correlation
where the absolute value exceeds 0.2 and any correlation with a p value less
than 0.2) and use `replicate`

to run the function 1000 times with 10 predictors
and \( n = 1000 \), we can calculate the probability that the different method
selects the correct model. (Note that we set all of values of `betas`

to 1, so
all 11 coefficient estimates should be non-zero and all of the predictors are
correlated with the outcome).

```
results <- replicate(1000, (simulateSelection(10, 1000)))
correctModelFraction <- apply(results == 11, 1, mean)
correctModelFraction
```

```
## [1] 1.000 0.026 0.573
```

Using a filter based on the size of the correlation only yields the correct model 2.6% of the time. Using a relatively high p value as a threshold does better returning the correct model 57.3% of the time.

We can explore the effect of varying the value of `pLim`

using a few simple
modifications to the above code.

```
simulateSelectionP <- function(nx, n, pLim) {
require(MASS)
require(clusterGeneration)
betas <- rep(1, nx + 1)
mu <- rep(0, nx)
sigma <- genPositiveDefMat(nx, covMethod = "unifcorrmat")$Sigma
X <- mvrnorm(n, mu, sigma)
y <- genY(betas, X)
pX <- X[, apply(X, 2, sigCor, y, pLim)]
correct <- (ncol(pX)) == nx
}
varyP <- function(nx, n, p) {
fractionCorrect <- numeric(100 * length(p))
j <- 1
for (i in 1:length(p)) {
set.seed(42)
fractionCorrect[(j):(j + 99)] <- replicate(100, simulateSelectionP(nx,
n, p[i]))
j <- j + 100
}
fractionCorrect
}
```

Using these new functions, we can calculate the fraction of the time that a
threshold of `p`

yields the correct model. If we vary `p`

from \( 0.01--0.50 \) and
plot the resulting curve,

```
prCorrect <- varyP(10, 1000, seq(0.01, 0.5, length.out = 20))
results <- data.frame(prCorrect = prCorrect, limit = rep(seq(0.01, 0.5, length.out = 20),
each = 100))
meanByLimit <- aggregate(results$prCorrect, list(results$limit), mean)
sdByLimit <- aggregate(results$prCorrect, list(results$limit), sd)
meanByLimit$sd <- sdByLimit$x
meanByLimit$upper <- meanByLimit$x + qnorm(0.975) * meanByLimit$sd/sqrt(100)
meanByLimit$lower <- meanByLimit$x - qnorm(0.975) * meanByLimit$sd/sqrt(100)
ggplot(meanByLimit, aes(x = Group.1, y = x, ymin = lower, ymax = upper)) + geom_line() +
geom_ribbon(alpha = 0.25) + theme_bw() + scale_x_continuous("Threshold for p value") +
scale_y_continuous("Fraction of models in simulation that are correct")
```

we see that, as expected, the probability of getting the correct model goes up as we increase the threshold value. We have about a 75% shot at getting the right model using this method when we set the limit on the p value to be 0.50 compared to under 30% at a p value of 0.01.

However, this raises an interesting question: if we increase the threshold for the p value, what is the value of using the bivariate correlation at all?

```
randomRate <- function(n, pLim) {
x <- rnorm(n)
y <- rnorm(n)
sigCor(x, y, pLim)
}
p <- seq(0.01, 0.5, length.out = 20)
rate <- numeric(length(p))
for (i in 1:length(p)) {
rate[i] <- mean(replicate(1000, randomRate(1000, p[i])))
}
rbind(p, rate)
```

```
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## p 0.010 0.03579 0.06158 0.08737 0.1132 0.1389 0.1647 0.1905 0.2163
## rate 0.009 0.03900 0.06400 0.10100 0.1110 0.1290 0.1660 0.1870 0.2210
## [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]
## p 0.2421 0.2679 0.2937 0.3195 0.3453 0.3711 0.3968 0.4226 0.4484 0.4742
## rate 0.2250 0.2920 0.3180 0.3180 0.3210 0.4010 0.3720 0.4560 0.4130 0.4660
## [,20]
## p 0.5
## rate 0.5
```

As the threshold increases, the value of the bivariate correlation as a filter goes down and goes down fast. At the limit required for there to be a 50% chance of the model containing all relevant predictors, the rate of improperly included variables would go up to nearly 20%.

The real problem comes into play here — the bivariate comparsions selects for the wrong variables by over-emphasizing the relationship between the marginal distributions. Because the selected model is incomplete and important variables are omitted, the resulting parameter estimates are biased and inaccurate.

The bivariate comparsion is a terrible way to select relevant variables for a highly dimensional model as the function of interest is relating the all of the predictors to the outcome. It can neither rule in nor rule out a predictor from a model. A better approach is to generate the model using domain knowledge and a model of the expected data generating process. If the process is to remain strictly data driven, methods like cross-validation or AIC/BIC provide better measures of model quality and predictor importance than the bivariate correlation coefficient.