XonlyPCA.Rmd

---
title: "Principal Components Regression, Pt.1: The Standard Method"
author: "Nina Zumel, Win-Vector LLC"
date: "May 7, 2016"
output:
  md_document:
    variant: markdown_github
---

# Principal Components Regression, Pt.1: The Standard Method

This article is by [Dr. Nina Zumel](http://www.win-vector.com/site/staff/nina-zumel/) of [Win-Vector LLC](http://www.win-vector.com/) and is hosted at: [http://www.win-vector.com/blog/2016/05/pcr_part1_xonly](http://www.win-vector.com/blog/2016/05/pcr_part1_xonly).


In this note, we discuss principal components regression and some of the issues with it:

* The need for scaling.
* The need for pruning.
* The lack of "_y_-awareness" of the standard dimensionality reduction step.

The purpose of this article is to set the stage for presenting dimensionality reduction techniques more appropriate for predictive modeling, such as _y_-aware principal components analysis, variable pruning, L2-regularized regression, supervised PCR, or partial least squares.  We do this by working detailed examples and building the relevant graphs.  In our follow-up article we describe and demonstrate the idea of _y_-aware scaling.

Note we will try to say "principal components" (plural) throughout, following Everitt's _The Cambridge Dictionary of Statistics_, though this is not the only common spelling  (e.g. Wikipedia: [Principal component regression](https://en.wikipedia.org/wiki/Principal_component_regression)).  We will work all of our examples in [R](https://cran.r-project.org/).
 
## Principal Components Regression

In principal components regression (PCR), we use principal components analysis (PCA) to decompose the independent (_x_) variables into an orthogonal basis (the principal components), and select a subset of those components as the variables to predict _y_. PCR and PCA are useful techniques for dimensionality reduction when modeling, and are especially useful when the independent variables are highly colinear.

Generally, one selects the principal components with the highest variance -- that is, the components with the largest singular values -- because the subspace defined by these principal components captures most of the variation in the data, and thus represents a smaller space that we believe captures most of the qualities of the data. Note, however, that standard PCA is an "_x_-only" decomposition, and as Jolliffe (1982) shows through examples from the literature, sometimes lower-variance components can be critical for predicting _y_, and conversely, high variance components are sometimes not important.

> Mosteller and Tukey (1977, pp. 397-398) argue similarly that the components with small variance are unlikely to be important in regression, apparently on the basis that nature is "tricky, but not downright mean". We shall see in the examples below that without too much effort we can find examples where nature is "downright mean".
-- Jolliffe (1982)

The remainder of this note presents principal components analysis in the context of PCR and predictive modeling in general. We will show some of the issues in using an _x_-only technique like PCA for dimensionality reduction. In a follow-up note, we'll discuss some _y_-aware approaches that address these issues.

```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.width=7, fig.height=7,
                      echo=TRUE, warning=FALSE, message=FALSE)

library('ggplot2')
library('tidyr')
library('WVPlots') # devtools::install_github('WinVector/WVPlots',build_vignettes=TRUE)

barbell_plot = function(frame, xvar, ymin, ymax, colorvar=NULL) {
  if(is.null(colorvar)) {
    gplot = ggplot(frame, aes_string(x=xvar))
  } else {
    gplot = ggplot(frame, aes_string(x=xvar, color=colorvar))
  }
  
  gplot + geom_point(aes_string(y=ymin)) + 
    geom_point(aes_string(y=ymax)) +
    geom_linerange(aes_string(ymin=ymin, ymax=ymax)) +
    ylab("value")
}

dotplot_identity = function(frame, xvar, yvar, colorvar=NULL) {
  if(is.null(colorvar)) {
    gplot = ggplot(frame, aes_string(x=xvar, y=yvar, ymax=yvar))
  } else {
    gplot = ggplot(frame, 
                   aes_string(x=xvar, y=yvar, ymax=yvar, 
                              color=colorvar))
  }
  gplot + geom_point() + geom_linerange(aes(ymin=0))
}

extractProjection <- function(ndim,princ) {
  # pull off the rotation.  
  proj <- princ$rotation[,1:ndim] 
  # sign was arbitrary, so flip in convenient form
  for(i in seq_len(ndim)) {
    si <- sign(mean(proj[,i]))
    if(si!=0) {
      proj[,i] <- proj[,i]*si
    }
  }
  proj
}

rsq <- function(x,y) {
  1 - sum((y-x)^2)/sum((y-mean(y))^2)
}
```

First, let's build our example. In this sort of teaching we insist on toy or synthetic problems so we actually
_know_ the right answer, and can therefore tell which procedures are better at modeling the truth.

In this data set, there are two (unobservable) processes: one that produces the output `yA` and one that produces the output `yB`.  We only observe the mixture of the two: `y =  yA + yB + eps`, where `eps` is a noise term.  Think of `y` as measuring some notion of success and the `x` variables as noisy estimates of two different factors that can each drive success.  We'll set things up so that the first five variables (`r paste('x',formatC(1:5,width=2,flag=0),sep='.')`) have all the signal. The odd numbered variables correspond to one process (`yB`) and the even numbered variables correspond to the other (`yA`).

Then, to simulate the difficulties of real world modeling, we'll add lots of pure noise
variables (`noise*`). The noise variables are unrelated to our _y_ of interest -- but are related to other "y-style" processes that we are not interested in. As is common with good statistical counterexamples, the example looks like something that should not happen or that can be easily avoided.  Our point is that the data analyst is usually working with data just like this. 

Data tends to come from databases that must support many different tasks, so it is exactly the case that there may be columns or variables that are correlated to unknown and unwanted additional processes. The reason PCA can't filter out these noise variables is that without use of _y_, standard PCA has no way of knowing _what_ portion of the variation in each variable is important to the problem at hand and should be preserved.  This _can_ be fixed through domain knowledge (knowing which variables to use), variable pruning and _y_-aware scaling. Our next article will discuss these procedures; in this article we will orient ourselves with a demonstration of both what a good analysis and what a bad analysis looks like.

All the variables are also deliberately mis-scaled to model some of the difficulties of working with under-curated real world data.

```{r mkdata}
# build example where even and odd variables are bringing in noisy images
# of two different signals.
mkData <- function(n) {
  for(group in 1:10) {
    # y is the sum of two effects yA and yB
    yA <- rnorm(n)
    yB <- rnorm(n)
    if(group==1) {
      d <- data.frame(y=yA+yB+rnorm(n))
      code <- 'x'
    } else {
      code <- paste0('noise',group-1)
    }
    yS <- list(yA,yB)
    # these variables are correlated with y in group 1,
    # but only to each other (and not y) in other groups
    for(i in 1:5) {
      vi <- yS[[1+(i%%2)]] + rnorm(nrow(d))
      d[[paste(code,formatC(i,width=2,flag=0),sep='.')]] <- ncol(d)*vi
    }
  }
  d
}
```

Notice the copy of _y_ in the data frame has additional "unexplainable variance"
so only about 66% of the variation in _y_ is predictable.

Let's start with our train and test data.

```{r makedata}
# make data
set.seed(23525)
dTrain <- mkData(1000)
dTest <- mkData(1000)
```

Let's look at our outcome _y_ and a few of our variables.

```{r firstlook}
summary(dTrain[, c("y", "x.01", "x.02",
                   "noise1.01", "noise1.02")])
```

Usually we recommend doing some significance pruning on variables before moving on -- see [here](http://www.win-vector.com/blog/2014/02/bad-bayes-an-example-of-why-you-need-hold-out-testing/) for possible consequences of not pruning an over-abundance of variables, and [here](http://www.win-vector.com/blog/2015/08/how-do-you-know-if-your-data-has-signal/) for a discussion of one way to prune, based on significance. For this example, however, we will deliberately attempt dimensionality reduction without pruning (to demonstrate the problem). Part of what we are trying to show is to _not_ assume PCA performs these steps for you.

## Ideal situation

First, let's look at the ideal situation. If we had sufficient domain knowledge (or had performed significance pruning) to remove the noise, we would have no pure noise variables.  In our example we know which variables carry signal and
therefore can limit down to them before doing the PCA as follows.

```{r ideal}
goodVars <-  colnames(dTrain)[grep('^x.',colnames(dTrain))]
dTrainIdeal <- dTrain[,c('y',goodVars)]
dTestIdeal <-  dTrain[,c('y',goodVars)]
```

Let's perform the analysis and look at the magnitude of the singular values.

```{r idealsv}
# do the PCA
dmTrainIdeal <- as.matrix(dTrainIdeal[,goodVars])
princIdeal <- prcomp(dmTrainIdeal,center = TRUE,scale. = TRUE)

# extract the principal components
rot5Ideal <- extractProjection(5,princIdeal)

# prepare the data to plot the variable loadings
rotfIdeal = as.data.frame(rot5Ideal)
rotfIdeal$varName = rownames(rotfIdeal)
rotflongIdeal = gather(rotfIdeal, "PC", "loading",
                       dplyr::starts_with("PC"))
rotflongIdeal$vartype = ifelse(grepl("noise", 
                                     rotflongIdeal$varName),
                               "noise", "signal")

# plot the singular values
dotplot_identity(frame = data.frame(pc=1:length(princIdeal$sdev), 
                            magnitude=princIdeal$sdev), 
                 xvar="pc",yvar="magnitude") +
  ggtitle("Ideal case: Magnitudes of singular values")
```

The magnitudes of the singular values tell us that the first two principal components carry most of the signal. We can also look at the variable loadings of the principal components. The plot of the variable loadings is a graphical representation of the coordinates of the principal components. Each coordinate corresponds to the contribution of one of the original variables to that principal component.

```{r idealsvld}
dotplot_identity(rotflongIdeal, "varName", "loading", "vartype") + 
  facet_wrap(~PC,nrow=1) + coord_flip() + 
  ggtitle("x scaled variable loadings, first 5 principal components") + 
  scale_color_manual(values = c("noise" = "#d95f02", "signal" = "#1b9e77"))
```

We see that we recover the even/odd loadings of the original signal variables. `PC1` has the odd variables, and `PC2` has the even variables. The next three principal components complete the basis for the five original variables. 

Since most of the signal is in the first two principal components, we can look at the projection of the data into that plane, using color to code *y*.

```{r idealproj}
# signs are arbitrary on PCA, so instead of calling predict we pull out
# (and alter) the projection by hand
projectedTrainIdeal <-
  as.data.frame(dmTrainIdeal %*% extractProjection(2,princIdeal),
                                 stringsAsFactors = FALSE)
projectedTrainIdeal$y <- dTrain$y
ScatterHistN(projectedTrainIdeal,'PC1','PC2','y',
               "Ideal Data projected to first two principal components")
```

Notice that the value of _y_ increases both as we move up and as we move right. We have recovered two orthogonal features that each correlate with an increase in y (in general the signs of the principal components -- that is, which direction is "positive" -- are arbitrary, so without precautions the above graph can appear flipped).  Recall that we constructed the data so that the odd variables (represented by `PC1`) correspond to process _yB_ and the even variables (represented by `PC2`) correspond to process _yA_. We have recovered both of these relations in the figure.

This is why you rely on domain knowledge, or barring that, at least prune your variables.  For this example variable pruning would have gotten us to the above ideal case.  In our next article we will show how to perform the significance pruning.

## *X*-only PCA

To demonstrate the problem of *x*-only PCA on unpruned data in a predictive modeling situation, let's analyze the same data without limiting ourselves to the known good variables.  We are pretending (as is often the case) we don't have the domain knowledge indicating which variables are useful _and_ we have neglected to significance prune the variables before PCA. In our experience, this is a common mistake in using PCR, or, more generally, with using PCA in predictive modeling situations.

This example will demonstrate how you lose
modeling power when you don't apply the methods in a manner appropriate to your
problem.  Note that the appropriate method for your data may not match the doctrine of another field, as they may
have different data issues.

### The wrong way: PCA without any scaling 

We deliberately mis-scaled the original data when we generated it.  Mis-scaled data is a common problem in data science situations, but perhaps less common in carefully curated scientific situations. In a messy data situation like the one we are emulating, the best
practice is to re-scale the *x* variables; however, we'll first naively apply PCA to the data as it is. This is to demonstrate the sensitivity of PCA to the units of the data.

```{r noscale }
vars <- setdiff(colnames(dTrain),'y')

duTrain <- as.matrix(dTrain[,vars])
prinU <- prcomp(duTrain,center = TRUE,scale. = FALSE) 

dotplot_identity(frame = data.frame(pc=1:length(prinU$sdev), 
                            magnitude=prinU$sdev), 
                 xvar="pc",yvar="magnitude") +
  ggtitle("Unscaled case: Magnitudes of singular values")

```

There is no obvious knee in the magnitudes of the singular values, so we are at a loss as to how many variables we should use. In addition, when we look at the variable loading of the first five principal components, we will see another problem:

```{r noscaleloading}
rot5U <- extractProjection(5,prinU)
rot5U = as.data.frame(rot5U)
rot5U$varName = rownames(rot5U)
rot5U = gather(rot5U, "PC", "loading",
                       dplyr::starts_with("PC"))
rot5U$vartype = ifelse(grepl("noise", 
                                     rot5U$varName),
                               "noise", "signal")

dotplot_identity(rot5U, "varName", "loading", "vartype") + 
  facet_wrap(~PC,nrow=1) + coord_flip() + 
  ggtitle("unscaled variable loadings, first 5 principal components") + 
  scale_color_manual(values = c("noise" = "#d95f02", "signal" = "#1b9e77"))
```

The noise variables completely dominate the loading of the first several principal components. Because of the way we deliberately mis-scaled the data, the noise variables are of much larger magnitude than the signal variables, and so the true signal is masked when we decompose the data.

Since the magnitudes of the singular values don't really give us a clue as to how many components to use in our model, let's try using all of them. This actually makes no sense, because using all the principal components is equivalent to using all the variables, thus defeating the whole purpose of doing PCA in the first place. But let's do it anyway (as many unwittingly do).

```{r noscalemodel}
# get all the principal components
# not really a projection as we took all components!
projectedTrain <- as.data.frame(predict(prinU,duTrain),
                                 stringsAsFactors = FALSE)
vars = colnames(projectedTrain)
projectedTrain$y <- dTrain$y

varexpr = paste(vars, collapse="+")
fmla = paste("y ~", varexpr)

model <- lm(fmla,data=projectedTrain)
summary(model)
estimate <- predict(model,newdata=projectedTrain)
trainrsq <- rsq(estimate,projectedTrain$y)
```

Note that most of the variables that achieve significance are the very last ones! We will leave it to the reader to confirm that using even as many as the first 25 principal components -- half the variables -- explains little of the variation in *y*. If we wanted to use PCR to reduce the dimensionality of the problem, we have failed. This is an example of what Jolliffe would have called a "downright mean" modeling problem, which we caused by mis-scaling the data. Note the r-squared of `r format(trainrsq, digits=4)` for comparison, later.

So now let's do what we should have done in the first place: scale the data.

### A better way: Preparing the training data with _x_-only scaling

Standard practice is to center the data at mean zero and scale it to unit standard deviation, which is easy with the `scale` command.

```{r xonlyexample}
dTrainNTreatedUnscaled <- dTrain
dTestNTreatedUnscaled <- dTest

# scale the data
dTrainNTreatedXscaled <- 
  as.data.frame(scale(dTrainNTreatedUnscaled[,colnames(dTrainNTreatedUnscaled)!='y'],
                      center=TRUE,scale=TRUE),stringsAsFactors = FALSE)
dTrainNTreatedXscaled$y <- dTrainNTreatedUnscaled$y
dTestNTreatedXscaled <- 
  as.data.frame(scale(dTestNTreatedUnscaled[,colnames(dTestNTreatedUnscaled)!='y'],
                      center=TRUE,scale=TRUE),stringsAsFactors = FALSE)
dTestNTreatedXscaled$y <- dTestNTreatedUnscaled$y

# get the variable ranges
ranges = vapply(dTrainNTreatedXscaled, FUN=function(col) c(min(col), max(col)), numeric(2))
rownames(ranges) = c("vmin", "vmax") 
rframe = as.data.frame(t(ranges))  # make ymin/ymax the columns
rframe$varName = rownames(rframe)
varnames = setdiff(rownames(rframe), "y")
rframe = rframe[varnames,]
rframe$vartype = ifelse(grepl("noise", rframe$varName),
                        "noise", "signal")

summary(dTrainNTreatedXscaled[, c("y", "x.01", "x.02", 
                                  "noise1.01", "noise1.02")])
barbell_plot(rframe, "varName", "vmin", "vmax", "vartype") +
  coord_flip() + ggtitle("x scaled variables: ranges") + 
  scale_color_manual(values = c("noise" = "#d95f02", "signal" = "#1b9e77"))
```

Note that the signal and noise variables now have commensurate ranges.

### The principal components analysis

```{r xscaledPCA}
vars = setdiff(colnames(dTrainNTreatedXscaled), "y")

dmTrain <- as.matrix(dTrainNTreatedXscaled[,vars])
dmTest <- as.matrix(dTestNTreatedXscaled[,vars])
princ <- prcomp(dmTrain,center = TRUE,scale. = TRUE) 
dotplot_identity(frame = data.frame(pc=1:length(princ$sdev), 
                            magnitude=princ$sdev), 
                 xvar="pc",yvar="magnitude") +
  ggtitle("x scaled variables: Magnitudes of singular values")
sum(princ$sdev^2)
```

Now the magnitudes of the singular values suggest that we can try to model the data with only the first twenty principal components. But first, let's look at the variable loadings of the first five principal components.

```{r xscaledload}
rot5 <- extractProjection(5,princ)
rotf = as.data.frame(rot5)
rotf$varName = rownames(rotf)
rotflong = gather(rotf, "PC", "loading", dplyr::starts_with("PC"))
rotflong$vartype = ifelse(grepl("noise", rotflong$varName), 
                          "noise", "signal")

dotplot_identity(rotflong, "varName", "loading", "vartype") + 
  facet_wrap(~PC,nrow=1) + coord_flip() + 
  ggtitle("x scaled variable loadings, first 5 principal components") + 
  scale_color_manual(values = c("noise" = "#d95f02", "signal" = "#1b9e77"))
```

The signal variables now have larger loadings than they did in the unscaled case, but the noise variables still dominate the projection, in aggregate swamping out the contributions from the signal variables. The two processes that produced *y* have diffused amongst the principal components, rather than mostly concentrating in the first two, as they did in the ideal case. This is because we constructed the noise variables to have variation and some correlations with each other -- but not be correlated with _y_. PCA doesn't know that we are interested only in variable correlations that are due to _y_, so it must decompose the data to capture as much variation, and as many variable correlations, as possible.

In other words, PCA must represent all processes present in the data, regardless of whether we are trying to predict those particular processes or not. Without the knowledge of the _y_ that we are trying to predict, PCA is forced to prepare for _any_ possible future prediction task.

#### Modeling 

Let's build a model using only the first twenty principal components, as our above analysis suggests we should. 

```{r quant2}
# get all the principal components
# not really a projection as we took all components!
projectedTrain <- as.data.frame(predict(princ,dmTrain),
                                 stringsAsFactors = FALSE)
projectedTrain$y <- dTrainNTreatedXscaled$y

ncomp = 20
# here we will only model with the first ncomp principal components
varexpr = paste(paste("PC", 1:ncomp, sep=''), collapse='+')
fmla = paste("y ~", varexpr)

model <- lm(fmla,data=projectedTrain)
summary(model)
projectedTrain$estimate <- predict(model,newdata=projectedTrain)
ScatterHist(projectedTrain,'estimate','y','Recovered 20 variable model versus truth (train)',
            smoothmethod='identity',annot_size=3)
trainrsq <- rsq(projectedTrain$estimate,projectedTrain$y)
```

This model explains `r format(100*trainrsq, digits=4)`% of the variation in the training set.  We do about as well on test.

```{r test}
projectedTest <- as.data.frame(predict(princ,dmTest),
                                 stringsAsFactors = FALSE)
projectedTest$y <- dTestNTreatedXscaled$y
projectedTest$estimate <- predict(model,newdata=projectedTest)
testrsq <- rsq(projectedTest$estimate,projectedTest$y)
testrsq
```

This is pretty good; recall that we had about 33% unexplainable variance in the data, so we would not expect any modeling algorithm to get better than an r-squared of about 0.67.

We can confirm that this performance is as good as simply regressing on all the variables without the PCA, so we have at least not lost information via our dimensionality reduction. 

```{r noprep}
# fit a model to the original data
vars <- setdiff(colnames(dTrain),'y')
formulaB <- paste('y',paste(vars,collapse=' + '),sep=' ~ ')
modelB <- lm(formulaB,data=dTrain)
dTrainestimate <- predict(modelB,newdata=dTrain)
rsq(dTrainestimate,dTrain$y)
dTestestimate <- predict(modelB,newdata=dTest)
rsq(dTestestimate,dTest$y)
```

We will show in our next article how to get a similar test r-squared from this data using a model with only two variables.

### Are we done?

Scaling the variables improves the performance of PCR on this data relative to not scaling, but we haven't completely solved the problem (though some analysts are fooled into thinking thusly). We have not explicitly recovered the two processes that drive *y*, and recovering such structure in the data is one of the purposes of PCA -- if we did not care about the underlying structure of the problem, we could simply fit a model to the original data, or use other methods (like significance pruning) to reduce the problem dimensionality.

It is a misconception in some fields that the variables must be orthogonal before fitting a linear regression model. This is *not* true. A linear model fit to collinear variables can still predict well; the only downside is that the coefficients of the model are not necessarily as easily interpretable as they are when the variables are orthogonal (and ideally, centered and scaled, as well). If your data has so much collinearity that the design matrix is ill-conditioned, causing the model coefficients to be inappropriately large or unstable, then regularization (ridge, lasso, or elastic-net regression) is a good solution. More complex predictive modeling approaches, for example random forest or gradient boosting, also tend to be more immune to collinearity.

So if you are doing PCR, you presumably are interested in the underlying structure of the data, and in this case, we haven't found it. Projecting onto the first few principal components fails to show much of a relation between these components and _y_.

We can confirm the first two _x_-scaled principal components are not informative with the following graph.

```{r xscaledplot}
proj <- extractProjection(2,princ)
# apply projection
projectedTrain <- as.data.frame(dmTrain %*% proj,
                      stringsAsFactors = FALSE)
projectedTrain$y <- dTrainNTreatedXscaled$y
# plot data sorted by principal components
ScatterHistN(projectedTrain,'PC1','PC2','y',
               "x scaled Data projected to first two principal components")
```

We see that _y_ is not well ordered by `PC1` and `PC2` here, as it was in the ideal case, and as it will be with the *y*-aware PCA. 

In our next article we will show that we can explain almost 50% of the _y_ variance in this data using only two variables.  This is quite
good as even the "all variable" model only picks up about that much of the relation and _y_ by design has about 33% unexplainable variation. In addition to showing the standard methods (including variable pruning) we will introduce a technique we call "_y_-aware scaling."

Click [here](http://www.win-vector.com/blog/2016/05/pcr_part2_yaware) for part 2 (*y*-aware methods).

### References

Everitt, B. S. _The Cambridge Dictionary of Statistics_, 2nd edition, Cambridge University Press, 2005.

Jolliffe, Ian T. "A Note on the Use of Principal Components in Regression," _Journal of the Royal Statistical Society. Series C (Applied Statistics)_, Vol. 31, No. 3 (1982), pp. 300-303