### R code from vignette source 'pcaMethods.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: pcaMethods.Rnw:102-114 ################################################### library(pcaMethods) x <- c(-4,7); y <- c(-3,4) distX <- rnorm(100, sd=0.3)*3 distY <- rnorm(100, sd=0.3) + distX * 0.3 mat <- cbind(distX, distY) res <- pca(mat, nPcs=2, method="svd", center=F) loading <- loadings(res)[1,] grad <- loading[2] / loading[1] if (grad < 0) grad <- grad * -1 lx <- c(-4,7) ly <- c(grad * -4, grad * 7) ################################################### ### code chunk number 2: pcaMethods.Rnw:118-125 ################################################### par(mar=c(2, 3, 2, 2)) plot(x,y, type="n", xlab="", ylab="") abline(v=0, col="dark gray", lwd = 2); abline(h=0, col = "dark gray", lwd = 2) points(distX, distY, type = 'p', col = "blue") lines(lx,ly, lwd = 2) points(-1, -1 * grad + 0.5, pch = 19, col = "red", lwd=4) points(6, 6 * grad + 0.5, pch = 19, col = "red", lwd=4) ################################################### ### code chunk number 3: pcaMethods.Rnw:253-255 ################################################### library(lattice) library(pcaMethods) ################################################### ### code chunk number 4: pcaMethods.Rnw:258-261 ################################################### library(pcaMethods) data(metaboliteData) data(metaboliteDataComplete) ################################################### ### code chunk number 5: pcaMethods.Rnw:264-266 ################################################### md <- prep(metaboliteData, scale="none", center=TRUE) mdC <- prep(metaboliteDataComplete, scale="none", center=TRUE) ################################################### ### code chunk number 6: pcaMethods.Rnw:271-277 ################################################### resPCA <- pca(mdC, method="svd", center=FALSE, nPcs=5) resPPCA <- pca(md, method="ppca", center=FALSE, nPcs=5) resBPCA <- pca(md, method="bpca", center=FALSE, nPcs=5) resSVDI <- pca(md, method="svdImpute", center=FALSE, nPcs=5) resNipals <- pca(md, method="nipals", center=FALSE, nPcs=5) resNLPCA <- pca(md, method="nlpca", center=FALSE, nPcs=5, maxSteps=300) ################################################### ### code chunk number 7: pcaMethods.Rnw:293-296 ################################################### sDevs <- cbind(sDev(resPCA), sDev(resPPCA), sDev(resBPCA), sDev(resSVDI), sDev(resNipals), sDev(resNLPCA)) matplot(sDevs, type = 'l', xlab="Eigenvalues", ylab="Standard deviation of PC", lwd=3) legend(x="topright", legend=c("PCA", "PPCA", "BPCA", "SVDimpute","Nipals PCA","NLPCA"), lty=1:6, col=1:6, lwd=3) ################################################### ### code chunk number 8: pcaMethods.Rnw:308-311 ################################################### par(mfrow=c(1,2)) plot(loadings(resBPCA)[,1], loadings(resPCA)[,1], xlab="BPCA", ylab="classic PCA", main = "Loading 1") plot(loadings(resBPCA)[,2], loadings(resPCA)[,2], xlab="BPCA", ylab="classic PCA", main = "Loading 2") ################################################### ### code chunk number 9: pcaMethods.Rnw:335-337 ################################################### q2SVDI <- Q2(resSVDI, mdC, fold=10) q2PPCA <- Q2(resPPCA, mdC, fold=10) ################################################### ### code chunk number 10: pcaMethods.Rnw:339-349 ################################################### # PPCA does not converge / misestimate a value in very rare cases. # This is a workaround to avoid that such a case will break the # diagram displayed in the vignette. # From the 2.0 release of bioconductor on, the convergence threshold # for PPCA was lowert to 1e-5, this should make the method much more # stable. So this workaround might be obsolete now... # [nope it is not, ppca is unstable] while( sum((abs(q2PPCA)) > 1) >= 1 ) { q2PPCA <- Q2(resPPCA, mdC, fold=10) } ################################################### ### code chunk number 11: pcaMethods.Rnw:353-356 ################################################### q2 <- data.frame(Q2=c(drop(q2PPCA), drop(q2SVDI)), method=c("PPCA", "SVD-Impute")[gl(2, 5)], PC=rep(1:5, 2)) print(xyplot(Q2~PC|method, q2, ylab=expression(Q^2), type="h", lwd=4)) ################################################### ### code chunk number 12: pcaMethods.Rnw:389-390 ################################################### errEsti <- kEstimate(md, method = "ppca", evalPcs=1:5, nruncv=1, em="nrmsep") ################################################### ### code chunk number 13: pcaMethods.Rnw:396-397 ################################################### barplot(drop(errEsti$eError), xlab="Loadings", ylab="NRMSEP (Single iteration)") ################################################### ### code chunk number 14: pcaMethods.Rnw:420-421 ################################################### barplot(drop(errEsti$variableWiseError[, which(errEsti$evalPcs == errEsti$bestNPcs)]), xlab="Incomplete variable Index", ylab="NRMSEP") ################################################### ### code chunk number 15: pcaMethods.Rnw:445-446 ################################################### slplot(resPCA) ################################################### ### code chunk number 16: pcaMethods.Rnw:456-457 ################################################### plotPcs(resPPCA, pc=1:3, type="score") ################################################### ### code chunk number 17: pcaMethods.Rnw:469-475 ################################################### pc <- pca(iris) irdf <- merge(iris, scores(pc), by=0) library(ggplot2) ggplot(irdf, aes(PC1, PC2, colour=Species)) + geom_point() + stat_ellipse()