# Multivariate regression or PCA to reduce response variables?

I hope the title is self-explanatory, but essentially I want to know which method is better: does it make sense to use a PCA to reduce a number of response Y variables and then conduct a univariate multiple regression OR conduct a multivariate multiple regression using all response variables as they are?

My variables are all continuous. I will be conducted the analysis in R.

## migrated from stackoverflow.comJul 24 '14 at 20:59

This question came from our site for professional and enthusiast programmers.

If I understand your question correctly, you have multiple response variables (y, not x), and you want to know if it makes more sense to use multivariate multiple regression on the responses directly, or to extract the first principal component of the y's first, and then regress on that.

As others have said, you are likely to get a much more rigorous answer on https://stats.stackexchange.com/, and I encourage you to try that, but hopefully this will help a little bit.

IMO the answer is: use mlm. There are two main reasons. First, you are likely to loose information about structure if you run PCA first, and second, the goodness-of-fit tests available for an mlm model are more robust than those available for a simple lm model (which is what you would end up with in the PCA approach).

To illustrate the first problem (in R), consider the following highly contrived dataset:

set.seed(1)     # for reproducible example
df <- data.frame(x1=sample(1:100,100),x2=sample(1:100,100),x3=sample(1:100,100))
df$y1 <- with(df,1+0.5*x1+3*x3+rnorm(100,sd=50)) df$y2 <- with(df,1+0.5*x2+3*x3+rnorm(100,sd=50))


In this dataset, we have two responses (y1 and y2), and three predictors (x1, x2, and x3). Both y1 and y2 depend strongly on x3 (by design); y1 also depends on x1 but not x2, and y2 depends on x2 but not x1. Because both y1 and y2 depend on x3, they are highly correlated.

If we run mlm on this dataset, we see the dependencies clearly.

# multivariate multiple regression
fit.mlm <- lm(cbind(y1,y2) ~ ., df)
summary(fit.mlm)
# Response y1 :
# ...
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)   1.3142    16.1056   0.082    0.935
# x1            0.5874     0.1849   3.177    0.002 **
# x2           -0.1303     0.1847  -0.705    0.482
# x3            3.0592     0.1857  16.470   <2e-16 ***
# ---
#
# Response y2 :
# ...
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)  17.3293    14.8825   1.164   0.2471
# x1           -0.2371     0.1709  -1.387   0.1685
# x2            0.4432     0.1707   2.597   0.0109 *
# x3            3.0696     0.1716  17.884   <2e-16 ***
# ---


So from this we can see clearly that x1 is important in explaining the variation in y1, and x2 is important in explaining the variation in y2. But look what happens when we run PCA first:

# principal components analysis
pc <- prcomp(df[c("y1","y2")],retx=T)
summary(pc)   # first pc explains 86% of variation in y
# Importance of components:
#                             PC1     PC2
# Standard deviation     136.9136 54.7439
# Proportion of Variance   0.8622  0.1378
# Cumulative Proportion    0.8622  1.0000
df$pc <- pc$x[,"PC1"]
fit.pca <- lm(pc~x1+x2+x3,df)
summary(fit.pca)
# ...
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)
# (Intercept) -242.6575    15.4029 -15.754   <2e-16 ***
# x1             0.2590     0.1768   1.464    0.146
# x2             0.2134     0.1766   1.208    0.230
# x3             4.3327     0.1776  24.391   <2e-16 ***
# ---


This result suggests that neither x1 nor x2 are important, which is false by design! In other words, we have lost the structural information about y1 and y2.

The second reason, regarding robustness of the goodness-of-fit tests is complex (there is a pretty good treatment here, while the classic paper on this is here). To motivate the treatment, consider two extreme cases:

In the first case, suppose the y`s are completely uncorrelated (orthogonal). In that case the principal components will just be the original y's, and PCA will contribute nothing.

In the second case, suppose the y's are perfectly correlated (colinear). In that case, the y's in a sense represent replicate measurements for each combination of the x's. The goodness of fit tests for mlm take into account these addition degrees of freedom, whereas an approach that uses PCA first essentially discards them. So in this case as well, you are better off with mlm.