2
$\begingroup$

In R there is this function pcr that runs a principal components regression. How can I use the information yielded from pcr in order to predict my response's value given a testset without using the function predict.pcr. I am very curious to find out. For example, below i have the first 6 components of the pcr_model.

                          1 comps       2 comps       3 comps       4 comps       5 comps      6 comps
X1                   1.402973e-05  1.638770e-05  3.025826e-05  3.843174e-05  3.893891e-05 1.058217e-04
X2                   3.585756e-06  6.721506e-06  1.424065e-05  3.192438e-05  5.830454e-05 4.251218e-04
X3                  -1.611619e-05 -1.697127e-05 -2.024606e-05 -2.847303e-05 -2.091408e-05 3.076982e-04
X4                  -1.747680e-06 -1.595572e-06 -2.301613e-07  9.525593e-05  2.357615e-04 2.249315e-04
X5                  -6.762417e-06 -8.728910e-06 -1.207538e-05  7.892716e-05  2.120044e-04 1.456909e-04
X6                  -1.069051e-06 -4.037689e-07 -4.417756e-07  1.255295e-04  2.152321e-05 4.000915e-05

So instead of going about predict(pcr_model,testset,ncomp=3) how would one go via a manual route? I was considering the approach as such: mean(1st comp * testset,2nd comp*testset,3rd comp*testset). Is this approach correct? If not what is the correct one?

EDIT: there is a nice tutorial online that uses predict.pcr function if you would like to see a quick example https://www.r-bloggers.com/performing-principal-components-regression-pcr-in-r/

$\endgroup$
  • $\begingroup$ I had to do this recently for something else. Where I would start is the GitHub for the pcr function, which is github.com/cran/pls/blob/…, and do some digging. $\endgroup$ – Clarinetist Oct 24 '17 at 13:19
  • $\begingroup$ i have looked at the pcr in Github. Either i did something wrong, for looking in the function i didn't get much help on how to go with the prediction way that I look for. $\endgroup$ – Hercules Apergis Oct 24 '17 at 14:45
  • 1
    $\begingroup$ You should be able to extract the following things from your pcr_model: mean of the training set, PCs (that you showed in the Q), and regression coefficients. For ncomp=3 and 6 original variables, mean mu is 6 dimensional vector, PCs are 6x3 matrix V, and regression coeffs is a 3 dimensional vector beta. Let's say your test data is stored in matrix Xtest. Then the prediction equals to: (Xtest-mu)*V*beta (this is pseudocode, you need to write down proper R code using matrix multiplication). It is as simple as that!! $\endgroup$ – amoeba Oct 24 '17 at 15:18
  • $\begingroup$ Mhmm, makes sense, the beta you say are the regression coeffs derived from running an lm with the three components? $\endgroup$ – Hercules Apergis Oct 24 '17 at 16:06
  • $\begingroup$ Yes. But I assume you can get it directly from the pcr_model object, you don't need to run lm yourself. Principal component regression (PCR) is by definition PCA followed by regression. So pcr() is doing nothing else than running first prcomp and then lm. $\endgroup$ – amoeba Oct 24 '17 at 21:12
2
$\begingroup$

Here's an example of doing it using entirely base R. It's not exactly manual -- I use the prcomp function. But it walks you through the main steps.

head(iris)
#outcome
y <- iris$Sepal.Length
#data
X <- model.matrix(Sepal.Length~.-1, data = iris)
#training and test sets
tr <- as.logical(1:nrow(iris) %% 2)
te <- tr == FALSE
#take principal components
pca <- prcomp(X[tr,], center = TRUE, scale = TRUE)
summary(pca)
# take first two
Xrot_tr <- pca$x[,1:2]
# do principal components regression
m <- lm(y[tr] ~ ., as.data.frame(Xrot_tr))
summary(m)
# predict for test set.  first rotate the test set according to the identified axes
Xrot_te <- predict(pca, newdata = X[te,])[,1:2]
# NB:  you could do the same thing with matrix multiplication.  First you need to scale the test set the same way as you scaled the training set
Xte <- sweep(sweep(X[te,], MARGIN = 2, FUN = "-", STATS = pca$center), MARGIN = 2, FUN = "/", STATS = pca$scale)
# then you multiply the scaled data by the rotation matrix
hard_way <- (Xte %*% pca$rotation)[,1:2]
all.equal(Xrot_te, hard_way)

# then get the predictions
pr <- predict(m, newdata = as.data.frame(Xrot_te))
# calc MSE
mean((y[te] - pr)^2)

# alternatively, you could predict without predict.lm
pr_manual <- cbind(1, hard_way) %*% as.matrix(m$coef)
head(pr_manual)
all.equal(pr, pr_manual)
# yeah whatever:
mean((pr - pr_manual)^2)
$\endgroup$
  • $\begingroup$ But still in order to go into predictions you use the predict function! But I want to go finding the predictions without the predict function! $\endgroup$ – Hercules Apergis Oct 24 '17 at 13:35
  • $\begingroup$ Oh that's super easy. Give me a sec $\endgroup$ – generic_user Oct 24 '17 at 13:36
  • $\begingroup$ Updated. In so doing I found a weird bug with linear models, and how they behave differently when fed data frames and formulas versus matrices. Ultimately however, I would have been lost if I didn't know the math behind everything $\endgroup$ – generic_user Oct 24 '17 at 14:01
  • $\begingroup$ +1 but this answer would improve if you included the same thing in math and in text. To make the prediction one needs to apply the centering and scaling from the training set, then the PCA rotation from the training set, then the linear combination from regression. Done. $\endgroup$ – amoeba Oct 24 '17 at 14:05
  • $\begingroup$ So how does the math go given that you have 3 components? As in the question; in pcr i get the components (with the variable coefficients). So when i have principal component 1,2 and 3, what is the math for deriving the prediction? (cause in your example you follow the prcomp function, while I am curious for the way pcr function which both do differently. $\endgroup$ – Hercules Apergis Oct 24 '17 at 14:05

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.