Regression through PCA

Question

I am working with principal component analysis (PCA) to have a regression model. Let's say we have 3 variables that we get every month, and so far we gathered 10 months of data. We have built the PCA on the raw data matrix.

p = prcomp(na.omit(rdatat), center=TRUE,scale=TRUE)

loadings=p$rotation[]

pca=cbind(p$x[,1], p$x[,2], p$x[,3])

Now we can put p$x[,1] as input of the regression and the result will be in the form of $Y=$ intercept$+ K\,$Component1, where $K$ is a coefficient from regression.

How can I build Component1 the month after in order to predict $Y$? I have the raw values of the 3 variables I use in PCA, but how can I use them?

Answer 1

prcomp has a predict S3 method you can use to apply the same transformations to new data quickly. Pass in the data for the new month and the prcomp object like so:

new.pca = predict(p, newdata=x.new)

But, the fact that you are asking this suggests that you are missing something fundamental about what PCA is doing, because you can also do this with the "loadings" rotation matrix. So I'm just going to give you a full PCA example so you can better understand what's going.

# sample data
set.seed(123)
n=500
mu=c(1,2)
sigma <- matrix(c(4,3,3,3), ncol=2)
x = rmvnorm(n, mu, sigma)

# This is what the "center" and "scale" prcomp options are doing for us.
sample.means = apply(x,2,mean) 
sample.sd    = apply(x,2,sd) 
x.scaled = (x-sample.means)/sample.sd

p = prcomp(x.scaled, center=F,scale=F)

x.scaled %*% p$rotation == p$x
# TRUE TRUE TRUE TRUE....

# new data
x.new = rmvnorm(10, mu, sigma)
x.new.scaled = (x.new-sample.means)/sample.sd

# The predict method does the transformation for us, which in this
# example is just multiplying by the rotation matrix. If we had asked
# prcomp to center and scale the data, the predict method would take
# care of this for us as well.
x.new.pca = predict(p, x.new.scaled) 

x.new.scaled %*% p$rotation == x.new.pca
# TRUE TRUE TRUE TRUE....

So what's so magical about this loadings matrix? Well, all PCA is doing is rotating our data into a new basis. That's why the components are always orthogonal to each other: they're the basis of the space after the transformation. I chose a two-dimensional example because it's much easier to see what's happening in two dimensions. Here's our data before we applied PCA (but after we scaled and centered):

plot(x.scaled, main='original, scaled data',xlim=c(-5,4.5), ylim=c(-3,3))

Here's the data after we applied PCA. As you can see, we just rotated it so that the direction along which we have the most variance (the first principle component) becomes our new x-axis.

plot(p$x, col='red', xlim=c(-5,4.5), ylim=c(-3,3))

To make this rotation really explicit, here's how the original location of our data maps to the new location after it's transformed by pca:

plot(x.scaled, main='Mapping of points from before and after rotation', xlim=c(-5,4.5), ylim=c(-3,3))
points(p$x, col='red')
sapply(1:n,function(i){
  lines(c(x.scaled[i,1], p$x[i,1]), c(x.scaled[i,2], p$x[i,2]))
})

(It's subtle, but you might have noticed that the transformed data is actually flipped along the x-axis. This is due to a negative PC, which isn't actually a problem)
And since we've gone this far, he's a sample regression on dimension-reduced data.

# let's do the regression you described to complete the example.
# sample data
set.seed(123)
n=500
mu=c(1,2)
sigma <- matrix(c(4,3,3,3), ncol=2)
x = rmvnorm(n, mu, sigma)      
b = c(2,3,4)
y = cbind(rep(1, n),x)%*%b
y.obs = y + rnorm(n, 0,5)    

# Dimensionality reduction
p = prcomp(x, scaled=T, centered=T)    
x.new = rmvnorm(10, mu, sigma)
x.new.pca = predict(p,x.new)  

# The predict method was giving me problems here so I just
# got the predicted y-values the old fashioned way.
y.new = cbind(rep(1,10),x.new.pca[,1]) %*% model.pca$coeff

# Plot
plot(p$x[,1], y.obs)
model.pca=lm(y.obs~p$x[,1])    
abline(model.pca, col="red", lty=2)
points(x.new.pca[,1], y.new, col="blue", pch=16 )