1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Sorry you wording/model is a bit confusing. Surely you want to say: $Y_i = b_{i,1} \psi_1 + ... + b_{i,3} \psi_3$, $\Psi$ being your matrix of PCs right? Check Massy (1965) Principal Components Regression in Exploratory Statistical Research and Jolliffe (1982) A Note on the Use of Principal Components in Regression; while slightly wordy they are easy to follow and will give you a good understanding about this kind of PCR-related questions. $\endgroup$ – usεr11852 Nov 13 '13 at 17:50
  • $\begingroup$ Do you mean to use a (negative) lag in the regressors? $\endgroup$ – Michael M Nov 13 '13 at 20:57
2
$\begingroup$

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))

enter image description here

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))

enter image description here

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]))
})

enter image description here

(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 )

enter image description here

$\endgroup$
  • $\begingroup$ Please note that in your concrete example the PCA rotation was also accompanied by reflection (which is because the determinant of the rotation matrix, the eigenvectors, happened to be -1. Taken that, your 3rd ("mapping") picture is messy due to the reflection. $\endgroup$ – ttnphns Nov 14 '13 at 4:55
  • $\begingroup$ Thanks for pointing that out. I actually noticed that after I'd posted the example but couldn't think of a good way to articulate what was going on without getting complicated. $\endgroup$ – David Marx Nov 14 '13 at 5:04
  • $\begingroup$ If you change sign of any one (or any odd number of) eigenvector (column of the rotation matrix), you get rid of the reflection. If you wish you might do it and repost the new pictures. But, to me, your current ones are good enough. $\endgroup$ – ttnphns Nov 14 '13 at 5:12
  • $\begingroup$ Thanks for the explanation. Insted of using the raw data it is possible to use the covariance matrix as input of prcomb, isn't it? $\endgroup$ – Francesco Nov 14 '13 at 11:41
  • $\begingroup$ @Francesco PCA operates on the covariance matrix already. The rotation matrix is just the eigenvectors of the covariance matrix of the data you provided. $\endgroup$ – David Marx Nov 15 '13 at 18:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.