Reducing the number of covariates - a couple of “simple” strategies

Retaining covariates with a “real” relationship

In the example we’ve been playing with so far, we know that only x1, x2, and x3 have a “real” relationship with y, because the process we used to generate the data has 0s for the coefficient relating x4-x9 to y. Let’s regenerate the data from our last example, with \(R^2 \approx 0.4\), restrict our analysis to only x1, x2, and x3 and see if we encounter the instability in parameter estimates and predictions we saw last time.

library(tidyverse)
library(reshape2)
library(ggplot2)
library(cowplot)
library(mvtnorm)
library(corrplot)

## intetcept
##
beta0 <- 1.0
## regression coefficients
##
beta <- c(1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
## pattern of correlation matrix, all non-zero entries are set to saem
## correlation, covariance matrix caldulated from individual variances and a 
## single association parameter governing the non-zero correlation coefficients
##
## Note: Not just any pattern will work here. The correlation matrix and
## covariance matrix generated from this pattern must be positive definite.
## If you change this pattern, you may get an error when you try to generate
## data with a non-zero association parameter.
##
Rho <- matrix(nrow = 9, ncol = , byrow = TRUE, 
              data = c(1,0,1,0,1,0,1,0,1,
                       0,1,0,1,0,1,0,1,0,
                       1,0,1,0,1,0,1,0,1,
                       0,1,0,1,0,1,0,1,0,
                       1,0,1,0,1,0,1,0,1,
                       0,1,0,1,0,1,0,1,0,
                       1,0,1,0,1,0,1,0,1,
                       0,1,0,1,0,1,0,1,0,
                       1,0,1,0,1,0,1,0,1
                       ))
## vector of standard deviations for covariates
##
sigma <- rep(1, 9)

## construct a covariance matrix from the pattern, standard deviations, and
## one parameter in [-1,1] that governs the magnitude of non-zero correlation
## coefficients
##
## Rho - the pattern of associations
## sigma - the vector of standard deviations
## rho - the association parameter
##
construct_Sigma <- function(Rho, sigma, rho) {
  ## get the correlation matris
  ##
  Rho <- Rho*rho
  for (i in 1:ncol(Rho)) {
    Rho[i,i] <- 1.0
  }
  ## notice the use of matrix multiplication
  ##
  Sigma <- diag(sigma) %*% Rho %*% diag(sigma)
  return(Sigma)
}

## set the random number seed manually so that every run of the code will 
## produce the same numbers
##
set.seed(1234)

n_samp <- 100
cov_str <- rmvnorm(n_samp,
                   mean = rep(0, nrow(Rho)),
                   sigma = construct_Sigma(Rho, sigma, 0.8))

resid <- rep(2.0, n_samp)

y_str <- rnorm(nrow(cov_str), mean = beta0 + cov_str %*% beta, sd = resid)
dat_1 <- data.frame(y_str, cov_str, rep("Strong", length(y_str)))

cov_str <- rmvnorm(n_samp,
                   mean = rep(0, nrow(Rho)),
                   sigma = construct_Sigma(Rho, sigma, 0.8))
y_str <- rnorm(nrow(cov_str), mean = beta0 + cov_str %*% beta, sd = resid)
dat_2 <- data.frame(y_str, cov_str, rep("Strong", length(y_str)))

column_names <- c("y", paste("x", seq(1, length(beta)), sep = ""), "Scenario")
colnames(dat_1) <- column_names
colnames(dat_2) <- column_names

Now that the data are generated, let’s run the analyses

lm_for_pred_1 <- lm(y ~ x1 + x2 + x3, data = dat_1)
lm_for_pred_2 <- lm(y ~ x1 + x2 + x3, data = dat_2)
summary(lm_for_pred_1)

Call:
lm(formula = y ~ x1 + x2 + x3, data = dat_1)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5057 -1.3596  0.0372  1.0451  4.9972 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8951     0.2076   4.311 3.93e-05 ***
x1            0.9918     0.3274   3.029 0.003150 ** 
x2           -0.5504     0.2232  -2.466 0.015419 *  
x3            1.0247     0.2942   3.483 0.000748 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.07 on 96 degrees of freedom
Multiple R-squared:  0.4006,    Adjusted R-squared:  0.3819 
F-statistic: 21.39 on 3 and 96 DF,  p-value: 1.079e-10

summary(lm_for_pred_2)

Call:
lm(formula = y ~ x1 + x2 + x3, data = dat_2)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5624 -1.3991 -0.0468  1.3177  4.1804 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.0206     0.2001   5.099 1.71e-06 ***
x1            0.6424     0.3398   1.891 0.061704 .  
x2           -1.0993     0.1900  -5.787 8.99e-08 ***
x3            1.1957     0.3199   3.737 0.000316 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.93 on 96 degrees of freedom
Multiple R-squared:  0.4848,    Adjusted R-squared:  0.4687 
F-statistic: 30.11 on 3 and 96 DF,  p-value: 8.283e-14

Well, the first thing we discover is that the parameter estimates are not as stable as we might hope. They’re all in the right direction - x1 and x3 have positive coefficients, x2 has a negative coefficient -, but the magnitude of x2 is a long way off in the first data set and the magnitude of x1 is a long way off in the second data set. If we compare predictions, here’s what we get:

new_data <- data.frame(x1 = 4.0, x2 = 4.0, x3 = 4.0,
                       x4 = 4.0, x5 = 4.0, x6 = 4.0,
                       x7 = 4.0, x8 = 4.0, x9 = 4.0)
pred_from_1 <- predict.lm(lm_for_pred_1, new_data)
pred_from_2 <- predict.lm(lm_for_pred_2, new_data)

cat("Prediction from data set 1: ", pred_from_1, "\n",
    "Prediction from data set 2: ", pred_from_2, "\n",
    "True answer:                ", beta0 + as.matrix(new_data) %*% beta, "\n", sep = "")

Prediction from data set 1: 6.7596
Prediction from data set 2: 3.975488
True answer:                5

That doesn’t look very good either. This isn’t good news, because we know that x1, x2, and x3 have an association with y and even though we excluded irrelevant covariates our parameter estimates aren’t very stable and our predictions can be pretty variable. Let’s try another approach.

Selecting (relatively) uncorrelated covariates

One thing we could observe from the data even if we didn’t know how it was constructed is how the covariates are associated with one another. Let’s visualize those associations using corrplot() on each of the data sets.

corrplot(cor(dat_1[,2:10]), order = "hclust", title = "Data set 1")

corrplot(cor(dat_2[,2:10]), order = "hclust", title = "Data set 2")

Not surprisingly, both data sets show pretty clearly that there are two sets of coefficients within which there is a high correlation and between which there’s very little correlation at all.¹ Let’s see what happens if we take one covariate from the first cluster and one from the second cluster, say x2 and x1. You should get similar results regardless of which covariate you pick from each cluster, but I’ll let you check that out on your own.

lm_for_pred_1 <- lm(y ~ x1 + x2, data = dat_1)
lm_for_pred_2 <- lm(y ~ x1 + x2, data = dat_2)
summary(lm_for_pred_1)

Call:
lm(formula = y ~ x1 + x2, data = dat_1)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5904 -1.2933 -0.0243  1.3791  5.3443 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8728     0.2191   3.984 0.000131 ***
x1            1.7472     0.2590   6.745 1.11e-09 ***
x2           -0.5170     0.2354  -2.196 0.030459 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.185 on 97 degrees of freedom
Multiple R-squared:  0.3249,    Adjusted R-squared:  0.311 
F-statistic: 23.34 on 2 and 97 DF,  p-value: 5.302e-09

summary(lm_for_pred_2)

Call:
lm(formula = y ~ x1 + x2, data = dat_2)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.1753 -1.5097  0.0831  1.5921  4.2268 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.9310     0.2116   4.400 2.78e-05 ***
x1            1.5923     0.2401   6.632 1.88e-09 ***
x2           -1.0565     0.2019  -5.232 9.66e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.055 on 97 degrees of freedom
Multiple R-squared:  0.4099,    Adjusted R-squared:  0.3977 
F-statistic: 33.68 on 2 and 97 DF,  p-value: 7.789e-12

new_data <- data.frame(x1 = 4.0, x2 = 4.0, x3 = 4.0,
                       x4 = 4.0, x5 = 4.0, x6 = 4.0,
                       x7 = 4.0, x8 = 4.0, x9 = 4.0)
pred_from_1 <- predict.lm(lm_for_pred_1, new_data)
pred_from_2 <- predict.lm(lm_for_pred_2, new_data)

cat("Prediction from data set 1: ", pred_from_1, "\n",
    "Prediction from data set 2: ", pred_from_2, "\n",
    "True answer:                ", beta0 + as.matrix(new_data) %*% beta, "\n", sep = "")

Prediction from data set 1: 5.793525
Prediction from data set 2: 3.07432
True answer:                5

That’s a bit better, but not a lot. Just for the record, a Bayesian version of the analysis gives similar results.

library(rstanarm)

Loading required package: Rcpp
rstanarm (Version 2.18.2, packaged: 2018-11-08 22:19:38 UTC)
- Do not expect the default priors to remain the same in future rstanarm versions.
Thus, R scripts should specify priors explicitly, even if they are just the defaults.
- For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores())
- Plotting theme set to bayesplot::theme_default().

options(mc.cores = parallel::detectCores())

stan_lm_for_pred_1 <- stan_glm(y ~ x1 + x2, data = dat_1, family = gaussian(), refresh = 0)
stan_lm_for_pred_2 <- stan_glm(y ~ x1 + x2, data = dat_2, family = gaussian(), refresh = 0)
summary(stan_lm_for_pred_1, digits = 3)

Model Info:

 function:     stan_glm
 family:       gaussian [identity]
 formula:      y ~ x1 + x2
 algorithm:    sampling
 priors:       see help('prior_summary')
 sample:       4000 (posterior sample size)
 observations: 100
 predictors:   3

Estimates:
                mean     sd       2.5%     25%      50%      75%      97.5% 
(Intercept)      0.875    0.218    0.452    0.727    0.875    1.021    1.311
x1               1.744    0.254    1.247    1.579    1.745    1.913    2.238
x2              -0.512    0.232   -0.956   -0.670   -0.513   -0.353   -0.067
sigma            2.200    0.162    1.908    2.088    2.189    2.301    2.543
mean_PPD         0.901    0.314    0.277    0.688    0.902    1.111    1.516
log-posterior -227.611    1.414 -231.265 -228.243 -227.265 -226.592 -225.875

Diagnostics:
              mcse  Rhat  n_eff
(Intercept)   0.003 0.999 4858 
x1            0.004 1.000 4197 
x2            0.003 1.000 4573 
sigma         0.002 1.000 4911 
mean_PPD      0.005 1.000 3894 
log-posterior 0.034 1.002 1706 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

summary(stan_lm_for_pred_2, digits = 3)

Model Info:

 function:     stan_glm
 family:       gaussian [identity]
 formula:      y ~ x1 + x2
 algorithm:    sampling
 priors:       see help('prior_summary')
 sample:       4000 (posterior sample size)
 observations: 100
 predictors:   3

Estimates:
                mean     sd       2.5%     25%      50%      75%      97.5% 
(Intercept)      0.933    0.220    0.498    0.787    0.934    1.085    1.355
x1               1.588    0.245    1.105    1.422    1.590    1.759    2.066
x2              -1.058    0.201   -1.460   -1.195   -1.057   -0.926   -0.655
sigma            2.070    0.152    1.798    1.966    2.061    2.163    2.401
mean_PPD         1.336    0.300    0.750    1.134    1.334    1.533    1.919
log-posterior -221.590    1.482 -225.291 -222.307 -221.224 -220.513 -219.769

Diagnostics:
              mcse  Rhat  n_eff
(Intercept)   0.003 1.000 5020 
x1            0.004 1.000 4604 
x2            0.003 1.000 4388 
sigma         0.002 1.000 4219 
mean_PPD      0.005 0.999 4417 
log-posterior 0.036 1.001 1696 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Conclusion

What have we learned so far? Here’s my quick summary:

• Multiple regression is useful. It may under the appropriate circumstances allow us to distinguish “real” from “spurious” associations.

• There are some significant challenges to using multiple regression. Specifically, if you include all of the covariates you’ve measured and they are highly correlated

  • Your estimates of association are likely to be unstable.²

  • As a result, predictions you make from your regression that extrapolate beyond the bounds of the data you fit are also likely to be unstable.

• Even if you know ahead of time which covariates are important,³ your estimates of the magnitude of any association are likely to be very imprecise.

• And if you restrict you covariates to those that are roughly independent, you aren’t much better off.⁴

So what are we to do? Well, there are some other alternatives to reducing the number of covariates that we’ll explore in notebooks yet to come.