Using PC scores or cluster analsis in predictions

Question

I have very big data and low number of observations. So I decided to use PCA to reduce dimension of the data. The following is R example (just an dummy example - for workout):

xmat <- matrix(sample(-1:1, 100000, replace = TRUE), ncol = 1000)
colnames(xmat) <- paste ("V", 1:1000, sep ="")
rownames(xmat) <- paste("S", 1:100, sep = "")

In this example dataset I have 1000 variables and 100 observations / subjects.

I am doing PCA. Lets say.

out <- princomp(xmat)
Error in princomp.default(xmat) : 
  'princomp' can only be used with more units than variables

Q1: is there a way to reduce dimensionality with p > n ? I would like to use all variables information as opposed to representative ones. Without having proper solution I went anyway to use cluster analysis of variables to categorize the variables and pick the randomly from the clusters.

To create a list of representative variables I tried to cluster the variables.

# cluster variables 
d <- dist(t(xmat), method = "euclidean") # distance matrix
fit <- hclust(d, method="ward")
plot(fit)
groups = cutree(fit,40)
groupd <- data.frame(var = names(groups), group = groups)

What I am thinking is randomly pick one variable from each group above and use this in PCA. Assume that I have the following y variable.

set.seed(1234)
yvar.d <- data.frame (subject = c(paste("S", 1:100, sep = "")), yvar = rnorm (100, 50,10))

Here is my question:

1. What could be statistical challenge of using cluster analysis ?

2. Can we use PCA scores in predictions of y. How ? Just multiple regression or we can introduce something such as variance explained by each components in the model ?

   Edits:

   Based on the discussions (see the comments below), I am using different function to do PC analysis.

"The calculation is done by a singular value decomposition of the (centered and possibly scaled) data matrix, not by using eigen on the covariance matrix. This is generally the preferred method for numerical accuracy. The print method for these objects prints the results in a nice format and the plot method produces a scree plot." - from function help.

     out1 <- prcomp(xmat)
      out1$x[1:3,1:3]
                      PC1        PC2       PC3
S1  2.940862 -2.7379835  6.527103
S2 -1.081124 -0.5294796 -0.276591
S3  2.375710  0.4505205 -4.236289

   out1$sdev
 screeplot(out1,npcs=30, type="lines",col=3) # 30 PCA plotted

 out1$rotation

I also come to see an example in SO how to use PCA in prediction. Here is my workout:

## take our training and test sets
YY <-  yvar.d$yvar 
prop <- 0.5
train = sample(1:length(YY), round(length(YY)*prop,0))


# data for testing model purpose 
testid = setdiff (1:length(YY), train)
YY1 <- YY
newXPCA <- data.frame(out1$x)
test.data <- data.frame (y = YY1[testid],newXPCA[testid,]) 
test.data[1:10,1:10]

train.data <- data.frame(y= YY1[train],newXPCA [train,])
train.data[1:10,1:10]

## fit the PCA
pc <- prcomp(train.data[, -1])
trainwPC <- data.frame (y = train.data$y, pc$x)

model1 <- lm(y ~ ., data = trainwPC)

#predict() method for class "prcomp"
test.p <- predict(pc, newdata = test.data)
pred <- predict(model1, newdata = data.frame(test.p), type = "response")
pred 
Warning message:
In predict.lm(model1, newdata = data.frame(test.p), type = "response") :
  prediction from a rank-deficient fit may be misleading

I just adopted this script from the SO link, I am not sure about accuracy of the script.

I still have technical questions remaining such as clarification to remaining question 2 above:

(1) If I want to split data into training and test set by sampling 50% of data (as show in the script). Should I do just multiple regression with y and the out1$x ? how many components to use ? is variance of each component play role in good model selection such as avoid over-fitting ? How ?

(2) Clustering (using x clusters) vs PCA analysis (with subset of x components vs all ) what would be statistically favorite for predictions in the situations where have p > n ? As I said to my mind the PCA analysis can use all information but I do not know if there is downside of such information such as over-fitting and "error consumption".

Worked example appreciated.

Answer 1

For (1), I would recommend 10-fold cross-validation. So, split the data (randomly) into 10 chunks:

cvGroup = sample(1:10,length(YY),replace=T)

Then, leave one of the cross-validation groups out and fit the model. Use the model built on the 9 groups to predict on the 10th. Repeat this process for each group (leaving it out of the model fitting process and using the final model to predict on this group). At the end, you'll now have predictions for every observation, and the model making the prediction did not use that observation. So, you can now use your favorite metric to measure model performance: mean squared error between predicted and actual, median absolute deviation, ROC (if Y is categorical), etc.

So, to address your question of "How many components should I use?", I would say "however many give you the best model". Start with 1, follow the procedure above, and get an MSE. Then, try two components. If it's not better than 1, stop and just use one. If it is, then try three. Continue until you've found the best. Alternatively, you could just try all of 1:n components and the plot the MSE as a function of the number of components. This should give you a feel of how many to use.

Cross-validation helps protect against over-fitting. The rationale is that if the model is truly over-fit to 9/10 of the data, then it won't perform well on the 1/10 that's left out. So, using this approach helps you to get a rough understanding of where over-fitting starts to become a problem.

For (2), I don't understand what you're asking. Are you wanting to cluster the x components to get a new variable? And then possibly build models for each clustered subset? You could try this, I'd again compare the MSE of this approach with the overall MSE and see what does best.

EDIT

Here's a example, using glm and MSE as the model and model performance metric, respectively:

set.seed(321)
xmat <- matrix(sample(-1:1, 100000, replace = TRUE), ncol = 1000)
colnames(xmat) <- paste ("V", 1:1000, sep ="")
rownames(xmat) <- paste("S", 1:100, sep = "")
yvar.d <- data.frame (subject = c(paste("S", 1:100, sep = "")), yvar = rnorm (100, 50,10))
new.vars = data.frame(prcomp(xmat)$x)
new.vars = cbind( Y=yvar.d[,2], new.vars )

cvGroup = sample(1:10, size=100, replace=T)
preds = rep(0, nrow(yvar.d))
RMSE = c()
for(varsUsed in 1:40){
  for(i in 1:10){
    mod = glm( Y ~ ., data=new.vars[cvGroup!=i,1:(1+varsUsed)] )
    preds[cvGroup==i] = predict(mod, newdata=new.vars[cvGroup==i,1:(1+varsUsed)])
  }
  RMSE = c(RMSE,mean( (preds-new.vars$Y)^2 ))
}
plot( RMSE )

Using set.seed() should give you the same results that I got. So, when you plot RMSE at the end, you should see an increasing curve: the more components you use in this model, the worse your model's RMSE is. This makes sense, because in this case there's no relationship between Y and the independent variables. You'd want to choose the model the number of components that gives you the lowest RMSE (or highest AUC, or lowest MAD, etc.)