Merge branch 'develop'

DESCRIPTION

@@ -1,12 +1,12 @@
 Package: mdatools
 Title: Multivariate Data Analysis for Chemometrics
-Version: 0.7.0
-Date: 2015-11-30
+Version: 0.7.1
+Date: 2016-05-02
 Author: Sergey Kucheryavskiy
 Maintainer: Sergey Kucheryavskiy <[email protected]>
 Description: Package implements projection based methods for preprocessing,
     exploring and analysis of multivariate data used in chemometrics.
 License: GPL-3
-Imports: methods
-RoxygenNote: 5.0.0
-
+Imports:
+    methods
+RoxygenNote: 5.0.1

NAMESPACE

@@ -1,4 +1,4 @@
-# Generated by roxygen2 (4.1.1): do not edit by hand
+# Generated by roxygen2: do not edit by hand
 
 S3method(as.matrix,classres)
 S3method(as.matrix,ldecomp)

NEWS

@@ -1,3 +1,11 @@
+v.0.7.1
+=======
+* fixed an issue lead to plot.new() error in some cases
+* documentation was regenerated with new version of Roxygen
+* file People.RData was renamed to people.RData
+* NIPALS method for PCA has been added
+* code optimization to speed calculations up
+
 v.0.7.0
 =======
 * interval PLS variable selection (iPLS) is implemented

R/ldecomp.R

@@ -244,8 +244,7 @@ ldecomp.getResLimits = function(eigenvals, nobj, ncomp, alpha = 0.05)
    {   
       if (i < nvar)
       {         
-         evals = eigenvals[(i + 1):nvar]         
-
+         evals = eigenvals[(i + 1):nvar]   
          cl = 2 * conflim - 100
          t1 = sum(evals)
          t2 = sum(evals^2)

R/mdaplots.R

@@ -1151,6 +1151,13 @@ mdaplotg = function(data, type = 'p', pch = 16,  lty = 1, lwd = 1, bwd = 0.8,
    if (is.null(yticks) && !is.null(yticklabels))
       yticks = 1:length(yticklabels)
 
+   # check if plot.new() should be called first
+   tryCatch(
+      {par(new = TRUE)},
+      warning = function(w){plot.new()},
+      finally = {par(new = FALSE)}
+   )
+
    # calculate limits and get colors
    lim = mdaplot.getAxesLim(data, show.lines = show.lines, single.x = single.x, 
                             show.legend = show.legend, legend = legend, show.labels = show.labels,

R/pca.R

@@ -18,13 +18,11 @@
 #' @param alpha
 #' significance level for calculating limit for Q residuals.
 #' @param method
-#' method to compute principal components.
+#' method to compute principal components ('svd', 'nipals').
 #' @param info
 #' a short text line with model description.
 #'
 #' @details 
-#' So far only SVD (Singular Value Decompisition) method is available, more coming soon. 
-#'   
 #' By default \code{pca} uses number of components (\code{ncomp}) as a minimum of number of 
 #' objects - 1, number of variables and default or provided value. Besides that, there is also 
 #' a parameter for selecting an optimal number of components (\code{ncomp.selected}). The optimal 
@@ -171,7 +169,7 @@ pca = function(x, ncomp = 15, center = T, scale = F, cv = NULL, x.test = NULL,
 #' 
 getCalibrationData.pca = function(obj, ...)
 {
-   x = obj$calres$scores %*% t(obj$loadings) + obj$calres$residuals
+   x = tcrossprod(obj$calres$scores, obj$loadings) + obj$calres$residuals
 
    if (is.numeric(attr(x, 'prep:scale')))
       x = sweep(x, 2L, attr(x, 'prep:scale'), '*', check.margin = F)
@@ -326,7 +324,7 @@ pca.mvreplace = function(x, center = T, scale = F, maxncomp = 7,
       scoresp = scores
       loadings = res$loadings[, 1:ncomp]      
       scores = x.rep %*% loadings
-      x.new = scores %*% t(loadings)   
+      x.new = tcrossprod(scores, loadings)   
 
       # remove centering
       x.new = sweep(x.new, 2L, lmean, '+', check.margin = F)
@@ -366,15 +364,21 @@ pca.mvreplace = function(x, center = T, scale = F, maxncomp = 7,
 #' @param scale
 #' logical, do standardization or not
 #' @param method
-#' algorithm for compiting PC space (only 'svd' is supported so far)
+#' algorithm for compiting PC space (only 'svd' and 'nipals' are supported so far)
 #' 
 #' @return
 #' an object with calibrated PCA model
 #' 
 pca.cal = function(x, ncomp, center = T, scale = F, method = 'svd')
 {
    x = prep.autoscale(x, center = center, scale = scale)
-   model = pca.svd(x, ncomp)
+
+   if (method == 'svd')
+      model = pca.svd(x, ncomp)
+   else if(method == 'nipals')
+      model = pca.nipals(x, ncomp)
+   else
+      stop('Wrong value for PCA method!')
 
    model$tnorm = sqrt(colSums(model$scores ^ 2)/(nrow(model$scores) - 1));   
 
@@ -456,23 +460,25 @@ pca.nipals = function(x, ncomp)
 
       while (th > 0.000001)
       {      
-         p = (t(E) %*% t) / as.vector((t(t) %*% t))
-         p = p / as.vector(t(p) %*% p) ^ 0.5
-         t = (E %*% p)/as.vector(t(p) %*% p)
-         th = abs(tau - as.vector(t(t) %*% t))
-         tau = as.vector(t(t) %*% t)
+         p = crossprod(E, t) / as.vector(crossprod(t))
+         p = p / as.vector(crossprod(p)) ^ 0.5
+         t = (E %*% p)/as.vector(crossprod(p))
+         th = abs(tau - as.vector(crossprod(t)))
+         tau = as.vector(crossprod(t))
       }
-            
-      E = E - t %*% t(p)
+
+      E = E - tcrossprod(t, p)
       scores[, i] = t
       loadings[, i] = p
       eigenvals[i] = tau / (nobj - 1)
    }
-
+
+   s = svd(E)
+
    res = list(
       loadings = loadings,
       scores = scores,
-      eigenvals = eigenvals
+      eigenvals = c(eigenvals, (s$d[1:(nvar - ncomp + 1)]^2)/(nrow(x) - 1))
    )   
 }
 
@@ -564,7 +570,7 @@ predict.pca = function(object, x, cv = F, ...)
 {
    x = prep.autoscale(x, object$center, object$scale)
    scores = x %*% object$loadings
-   residuals = x - scores %*% t(object$loadings)
+   residuals = x - tcrossprod(scores, object$loadings)
 
    if (cv == F)
    {   

R/pls.R

@@ -448,13 +448,12 @@ pls.simpls = function(x, y, ncomp, cv = FALSE)
    for (n in 1:ncomp)
    {
       # get the dominate eigenvector of A'A
-      e = eigen(t(A) %*% A)
+      e = eigen(crossprod(A))
       q = e$vectors[1:nresp]
 
-
       # calculate and store weights
       w = A %*% q
-      c = t(w) %*% M %*% w
+      c = crossprod(w, (M %*% w))
       w = w/sqrt(as.numeric(c))
       W[, n] = w
 
@@ -463,18 +462,18 @@ pls.simpls = function(x, y, ncomp, cv = FALSE)
       P[, n] = p
 
       # calculate and store y loadings
-      q = t(A) %*% w
+      q = crossprod(A, w)
       Q[, n] = q
 
       v = C %*% p
-      v = v/sqrt(as.numeric(t(v) %*% v))
+      v = v/sqrt(as.numeric(crossprod(v)))
 
       # calculate and store regression coefficients
-      B[, n, ] = W[, 1:n, drop = FALSE] %*% t(Q[, 1:n, drop = FALSE])
+      B[, n, ] = tcrossprod(W[, 1:n, drop = FALSE], Q[, 1:n, drop = FALSE])
 
       # recalculate matrices for the next compnonent
-      C = C - v %*% t(v)
-      M = M - p %*% t(p)
+      C = C - tcrossprod(v)
+      M = M - tcrossprod(p)
       A = C %*% A      
 
       if (cv == F && e$value < 10^-12) {

README.md

@@ -10,10 +10,10 @@ For more details and examples read a [GitBook tutorial](http://svkucheryavski.gi
 How to install
 --------------
 
-The package now is available from CRAN by usual installing procedure.  However due to restrictions in CRAN politics regarding number of submissions (one in 3-4 month) only major releases will be published there. To get the latest release plase use GitHub sources. You can [download](https://github.com/svkucheryavski/mdatools/releases) a zip-file with source package and install it using the `install.packages` command, e.g. if the downloaded file is `mdatools_0.7.0.tar.gz` and it is located in a current working directory, just run the following:
+The package now is available from CRAN by usual installing procedure.  However due to restrictions in CRAN politics regarding number of submissions (one in 3-4 month) only major releases will be published there. To get the latest release plase use GitHub sources. You can [download](https://github.com/svkucheryavski/mdatools/releases) a zip-file with source package and install it using the `install.packages` command, e.g. if the downloaded file is `mdatools_0.7.1.tar.gz` and it is located in a current working directory, just run the following:
 
 ```
-install.packages('mdatools_0.7.0.tar.gz')
+install.packages('mdatools_0.7.1.tar.gz')
 ```
 
 If you have `devtools` package installed, the following command will install the latest release from the master branch of GitHub repository (do not forget to load the `devtools` package first):