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:
- What could be statistical challenge of using cluster analysis ?
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 ?
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
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.
'princomp' can only be used with more units than variablesPCA can be used with p>n. Maybe you were using the wrong function or using it wrongly? $\endgroup$
princomp()performs PCA using a Eigendecomposition.
princomp()will fail if the number of variables is larger than the number of observations. Use the SVD based
prcomp()instead; this is generally the preferred method for numerical accuracy. $\endgroup$
Rfunctions, but, definitely, eigen-decomposition can be used in n<p case, too. $\endgroup$
stop("'princomp' can only be used with more units than variables")$\endgroup$