###################################################
### chunk number 1: 
###################################################
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 = res@loadings[1,]
grad = loading[2] / loading[1]
if (grad < 0)
   grad = grad * -1
lx = c(-4,7)
ly = c(grad * -4, grad * 7)


###################################################
### chunk number 2: 
###################################################
par(mar=c(2, 3, 2, 2))
plot(x,y, type="n", cex.axis=2, cex.lab=2, 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)


###################################################
### chunk number 3: 
###################################################
library(lattice)
library(pcaMethods)


###################################################
### chunk number 4: 
###################################################
library(pcaMethods)
data(metaboliteData)
data(metaboliteDataComplete)


###################################################
### chunk number 5: 
###################################################
md  <- prep(metaboliteData, scale="none", center=TRUE)
mdC  <- prep(metaboliteDataComplete, scale="none", center=TRUE)


###################################################
### chunk number 6: 
###################################################
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)


###################################################
### chunk number 7: 
###################################################
par(mar=c(5, 6, 4, 2))
sDevs <- cbind(resPCA@sDev, resPPCA@sDev, resBPCA@sDev, resSVDI@sDev, resNipals@sDev, resNLPCA@sDev)
matplot(sDevs, type = 'l', xlab="Eigenvalues", ylab="size", cex.lab=2.5, cex.main=1.5, lwd=3, cex.axis=2)
legend(x="topright", legend=c("PCA", "PPCA", "BPCA", "SVDimpute","Nipals PCA","NLPCA"), lty=1:6, col=1:6, cex=2, lwd=3)


###################################################
### chunk number 8: 
###################################################
par(mar=c(5, 6, 4, 2))
par(mfrow=c(1,2))
plot(resBPCA@loadings[,1], resPCA@loadings[,1], xlab="BPCA", ylab="classic PCA", main = "Loading 1", cex.lab = 2.5, cex.main=2.5, cex.axis = 2)
plot(resBPCA@loadings[,2], resPCA@loadings[,2], xlab="BPCA", ylab="classic PCA", main = "Loading 2", cex.lab = 2.5, cex.main=2.5, cex.axis = 2)


###################################################
### chunk number 9: 
###################################################
#q2BPCA <- Q2(resBPCA, mdC, nPcs=5, nruncv=1, fold=10)
#q2Nipals <- Q2(resNipals, mdC, nPcs=5, nruncv=1, fold=10)
q2SVDI <- Q2(resSVDI, mdC, nPcs=5, nruncv=1, fold=10)
q2PPCA <- Q2(resPPCA, mdC, nPcs=5, nruncv=1, fold=10)


###################################################
### chunk number 10: 
###################################################
# 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...
while( sum((abs(q2PPCA)) > 1) >= 1 ) {
    q2PPCA <- Q2(resPPCA, mdC, nPcs=5, nruncv=1, fold=10)
}


###################################################
### chunk number 11: 
###################################################
#q2 <- data.frame(Q2=c(drop(q2BPCA), drop(q2Nipals), drop(q2PPCA), drop(q2SVDI)), 
#                 method=c("BPCA", "Nipals", "PPCA", "SVD-Impute")[gl(4, 5)], PC=rep(1:5, 4))
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))


###################################################
### chunk number 12: 
###################################################
errEsti <- kEstimate(md, method = "ppca", evalPcs=1:5, nruncv=1, em="nrmsep")


###################################################
### chunk number 13: 
###################################################
par(mar=c(5, 6, 4, 2))
barplot(drop(errEsti$eError), xlab="Loadings", ylab="NRMSEP (Single iteration)", cex.lab=2, cex.axis=1.5)


###################################################
### chunk number 14: 
###################################################
par(mar=c(5, 6, 4, 2))
barplot(drop(errEsti$variableWiseError[, which(errEsti$evalPcs == errEsti$bestNPcs)]), 
xlab="Incomplete variable Index", ylab="NRMSEP", cex.lab=2, cex.axis=1.5)


###################################################
### chunk number 15: 
###################################################
slplot(resPCA)


###################################################
### chunk number 16: 
###################################################
plotPcs(resPPCA, pc=1:3, scoresLoadings=c(TRUE, FALSE))