old/example_virus.R

In kwstat/nipals: Principal Components Analysis using NIPALS or Weighted EMPCA, with Gram-Schmidt Orthogonalization

#

# tutorial
# https://nbviewer.jupyter.org/github/bertdv/AIP-5SSB0/blob/master/lessons/notebooks/11_Continuous-Latent-Variable-Models-PCA-and-FA.ipynb
# data source
# https://github.com/bertdv/AIP-5SSB0/blob/master/lessons/notebooks/datasets/virus3.dat
vir <- matrix(c(
  17.0,13.0,14.0,16.0,4.0,9.0,NA,1.0,13.0,0.0,11.0,13.0,NA,7.0,1.0,4.0,11.0,5.0,
  12.0,11.0,9.0,12.0,NA,5.0,12.0,1.0,9.0,1.0,7.0,12.0,5.0,6.0,0.0,NA,8.0,2.0,
  18.0,16.0,16.0,NA,8.0,6.0,14.0,1.0,NA,0.0,9.0,12.0,NA,8.0,0.0,2.0,11.0,NA,
  18.0,16.0,15.0,19.0,8.0,6.0,11.0,1.0,15.0,NA,7.0,13.0,5.0,8.0,0.0,2.0,9.0,3.0,
  17.0,13.0,13.0,NA,8.0,4.0,18.0,1.0,10.0,NA,8.0,11.0,7.0,6.0,1.0,NA,10.0,2.0,
  16.0,13.0,16.0,21.0,9.0,3.0,17.0,NA,10.0,NA,7.0,12.0,7.0,5.0,1.0,2.0,11.0,NA,
  22.0,NA,10.0,16.0,10.0,NA,18.0,NA,NA,2.0,8.0,11.0,6.0,8.0,0.0,1.0,8.0,NA,
  20.0,10.0,24.0,10.0,6.0,9.0,21.0,0.0,7.0,0.0,7.0,NA,4.0,9.0,1.0,NA,NA,2.0,
  20.0,21.0,12.0,15.0,NA,7.0,11.0,1.0,NA,3.0,8.0,14.0,6.0,7.0,NA,1.0,10.0,3.0,
  NA,21.0,NA,15.0,9.0,7.0,11.0,1.0,9.0,3.0,9.0,14.0,5.0,7.0,NA,1.0,10.0,3.0,
  NA,11.0,NA,10.0,9.0,6.0,19.0,0.0,12.0,NA,7.0,14.0,NA,NA,NA,4.0,9.0,1.0,
  20.0,12.0,23.0,10.0,8.0,5.0,20.0,0.0,13.0,0.0,NA,13.0,4.0,NA,NA,4.0,10.0,1.0,
  NA,NA,18.0,16.0,NA,4.0,12.0,0.0,12.0,0.0,10.0,15.0,8.0,6.0,1.0,1.0,12.0,1.0,
  17.0,NA,NA,15.0,8.0,6.0,14.0,1.0,14.0,0.0,9.0,12.0,NA,NA,0.0,3.0,NA,3.0,
  19.0,17.0,NA,NA,8.0,NA,14.0,1.0,14.0,0.0,8.0,12.0,NA,NA,0.0,2.0,12.0,NA,
  NA,NA,15.0,NA,8.0,5.0,14.0,1.0,14.0,NA,8.0,12.0,4.0,8.0,0.0,2.0,12.0,3.0,
  NA,15.0,16.0,16.0,8.0,6.0,14.0,1.0,15.0,NA,8.0,12.0,NA,8.0,NA,2.0,NA,3.0,
  17.0,17.0,16.0,19.0,NA,6.0,NA,1.0,NA,1.0,7.0,13.0,5.0,8.0,0.0,2.0,9.0,NA,
  18.0,17.0,NA,NA,8.0,6.0,11.0,1.0,NA,1.0,7.0,13.0,5.0,8.0,0.0,NA,9.0,3.0,
  22.0,NA,10.0,16.0,10.0,5.0,17.0,1.0,12.0,2.0,8.0,11.0,6.0,8.0,0.0,1.0,8.0,NA,
  18.0,NA,NA,18.0,NA,8.0,17.0,1.0,14.0,1.0,5.0,16.0,4.0,NA,0.0,NA,NA,2.0,
  NA,16.0,16.0,15.0,NA,6.0,13.0,1.0,14.0,1.0,8.0,12.0,4.0,8.0,1.0,2.0,12.0,3.0,
  20.0,21.0,12.0,NA,9.0,NA,11.0,1.0,10.0,3.0,NA,14.0,7.0,7.0,0.0,1.0,9.0,3.0,
  20.0,21.0,12.0,15.0,9.0,7.0,11.0,NA,10.0,3.0,9.0,14.0,5.0,7.0,0.0,1.0,10.0,3.0,
  18.0,12.0,NA,10.0,9.0,NA,NA,0.0,14.0,0.0,7.0,12.0,4.0,NA,0.0,NA,10.0,1.0,
  18.0,12.0,21.0,10.0,10.0,NA,18.0,0.0,13.0,0.0,8.0,12.0,4.0,NA,0.0,4.0,10.0,NA,
  17.0,12.0,22.0,10.0,8.0,5.0,18.0,0.0,NA,NA,5.0,13.0,4.0,10.0,0.0,3.0,NA,1.0,
  NA,16.0,16.0,NA,8.0,6.0,15.0,1.0,NA,NA,NA,12.0,4.0,8.0,0.0,2.0,11.0,NA,
  19.0,NA,15.0,17.0,7.0,NA,15.0,1.0,14.0,0.0,8.0,12.0,4.0,8.0,0.0,2.0,NA,3.0,
  NA,16.0,16.0,19.0,NA,6.0,11.0,1.0,15.0,1.0,NA,13.0,5.0,8.0,0.0,2.0,9.0,3.0,
  18.0,17.0,15.0,17.0,8.0,6.0,15.0,1.0,NA,0.0,8.0,12.0,4.0,8.0,0.0,3.0,9.0,3.0,
  15.0,12.0,14.0,23.0,8.0,3.0,17.0,1.0,9.0,NA,7.0,15.0,6.0,6.0,NA,2.0,11.0,2.0,
  13.0,11.0,14.0,22.0,NA,3.0,NA,1.0,10.0,4.0,8.0,13.0,6.0,6.0,1.0,3.0,11.0,2.0,
  16.0,11.0,15.0,23.0,10.0,4.0,18.0,1.0,10.0,3.0,7.0,12.0,6.0,NA,1.0,NA,9.0,3.0,
  14.0,NA,14.0,NA,11.0,3.0,19.0,2.0,10.0,2.0,7.0,12.0,6.0,5.0,NA,2.0,9.0,3.0,
  NA,NA,15.0,24.0,10.0,NA,18.0,1.0,11.0,1.0,7.0,NA,5.0,NA,2.0,3.0,11.0,2.0,
  15.0,NA,12.0,21.0,8.0,4.0,NA,1.0,10.0,3.0,7.0,15.0,7.0,6.0,1.0,3.0,10.0,3.0,
  15.0,11.0,NA,22.0,NA,3.0,19.0,1.0,8.0,3.0,4.0,14.0,6.0,5.0,1.0,2.0,10.0,NA),
  byrow=TRUE, ncol=18)
rownames(vir) <- 1:38

libs(nipals, pcaMethods)

# reasonable
m1 <- nipals(vir)
plot(m1$scores[,1], m1$scores[,2], type="n")
text(m1$scores[,1], m1$scores[,2], rownames(vir))

# lousy result.  Why?
m2 <- empca(vir, tol=1e-7)
m2$eig
plot(m2$scores[,1], m2$scores[,2], type="n")
text(m2$scores[,1], m2$scores[,2], rownames(vir))
# fewer PCs
m2 <- empca(vir, ncomp=10)
m2$eig
plot(m2$scores[,1], m2$scores[,2], type="n")
text(m2$scores[,1], m2$scores[,2], rownames(vir))
# need to compare to python results
m2 <- empca(vir, ncomp=10, maxiter=200, verbose=TRUE)
m2 <- empca(vir, ncomp=10, maxiter=200, verbose=TRUE)
plot(m2$scores[,1], m2$scores[,2], type="n")
text(m2$scores[,1], m2$scores[,2], rownames(vir))

# reasonable
m3 <- pca(vir, method="ppca")
biplot(m3)