cpfa {cpfa}  R Documentation 
Classification with Parallel Factor Analysis
Description
Fits Richard A. Harshman's Parallel Factor Analysis1 (Parafac) model or Parallel Factor Analysis2 (Parafac2) model to a threeway or fourway data array. Allows for different constraint options on multiple tensor modes. Uses Parafac component weights from a single mode of this model as predictors to tune parameters for one or more classification methods via a kfold crossvalidation procedure. Predicts class labels and calculates multiple performance measures for binary or multiclass classification over some number of replications with different traintest splits. Provides descriptive statistics to pool output across replications.
Usage
cpfa(x, y, model = c("parafac", "parafac2"), nfac = 1, nrep = 5, ratio = 0.8,
nfolds = 10, method = c("PLR", "SVM", "RF", "NN", "RDA", "GBM"),
family = c("binomial", "multinomial"), parameters = list(),
type.out = c("measures", "descriptives"), foldid = NULL,
prior = NULL, cmode = NULL, seeds = NULL, plot.out = FALSE,
plot.measures = NULL, parallel = FALSE, cl = NULL, verbose = TRUE, ...)
Arguments
x 
A threeway or fourway data array. For Parafac2, can be a list of length 
y 
A vector containing at least two unique class labels. Should be a factor that contains two or more levels . For binary case, ensure the order of factor levels (left to right) is such that negative class is first and positive class is second. 
model 
Character designating the Parafac model to use, either 
nfac 
Number of components for each Parafac or Parafac2 model to fit. Default is 
nrep 
Number of replications to repeat the procedure. Default is 
ratio 
Split ratio for dividing data into train and test sets. Default is 
nfolds 
Numeric setting number of folds for kfold crossvalidation. Must be 2 or greater. Default is 
method 
Character vector indicating classification methods to use. Possible methods include penalized logistic regression (PLR); support vector machine (SVM); random forest (RF); feedforward neural network (NN); regularized discriminant analysis (RDA); and gradient boosting machine (GBM). If none are selected, default is to use all methods with 
family 
Character value specifying binary classification ( 
parameters 
List containing arguments related to classification methods. When specified, must contain one or more of the following:

type.out 
Type of output desired: 
foldid 
Integer vector containing fold IDs for kfold crossvalidation. If not provided, fold IDs are generated randomly for number of folds 
prior 
Prior probabilities of class membership. If unspecified, the class proportions for input 
cmode 
Integer value of 1, 2, or 3 (or 4 if 
seeds 
Random seeds to be associated with each replication. Default is 
plot.out 
Logical indicating whether to output one or more box plots of classification performance measures that are plotted across classification methods and number of components. 
plot.measures 
Character vector containing values that specify for plotting one or more of 11 possible classification performance measures. Only relevant when 
parallel 
Logical indicating if parallel computing should be implemented. If TRUE, the package parallel is used for parallel computing. For all classification methods except penalized logistic regression, the doParallel package is used as a wrapper. Defaults to FALSE, which implements sequential computing. 
cl 
Cluster for parallel computing, which is used when 
verbose 
If TRUE, progress is printed. 
... 
Additional arguments to be passed to function 
Details
Data are split into a training set and a testing set. After fitting a Parafac or Parafac2 model with the training set using package multiway (see parafac
or parafac2
in multiway for details), the estimated classification mode weight matrix is passed to one or several of six classification methods. The methods include: penalized logistic regression (PLR); support vector machine (SVM); random forest (RF); feedforward neural network (NN); regularized discriminant analysis (RDA); and gradient boosting machine (GBM).
Package glmnet fits models for PLR. PLR tunes penalty parameter lambda while the elastic net parameter alpha is set by the user (see the help file for function cv.glmnet
in package glmnet). For SVM, package e1071 is used with a radial basis kernel. Penalty parameter cost and radial basis parameter gamma are used (see svm
in package e1071). For RF, package randomForest is used and implements Breiman's random forest algorithm. The number of predictors sampled at each node split is set at the default of sqrt(R), where R is the number of Parafac or Parafac2 components. Two tuning parameters allowed are ntree, the number of trees to be grown, and nodesize, the minimum size of terminal nodes (see randomForest
in package randomForest). For NN, package nnet fits a singlehiddenlayer, feedforward neural network model. Penalty parameters size (i.e., number of hidden layer units) and decay (i.e., weight decay) are used (see nnet). For RDA, package rda fits a shrunken centroids regularized discriminant analysis model. Tuning parameters include rda.alpha, the shrinkage penalty for the withinclass covariance matrix, and delta, the shrinkage penalty of class centroids towards the overall dataset centroid. For GBM, package xgboost fits a gradient boosting machine model. Four tuning parameters are allowed: (1) eta, the learning rate; (2) max.depth, the maximum tree depth; (3) subsample, the fraction of samples per tree; and (4) nrounds, the number of boosting trees to build.
For all six methods, kfold crossvalidation is implemented to tune classification parameters where the number of folds is set by argument nfolds
. Separately, the trained Parafac or Parafac2 model is used to predict the classification mode's component weights using the testing set data. The predicted component weights and the optimized classification method are then used to predict class labels. Finally, classification performance measures are calculated. The process is repeated over a number of replications with different random splits of the input array and of the class labels at each replication.
Value
Returns an object of class wrapcpfa
either with a threeway array with classification performance measures for each model and for each replication, or with a list containing matrices with descriptive statistics for performance measures calculated across all replications. Specify type.out = "measures"
to output the array of performance measures. Specify type.out = "descriptives"
to output descriptive statistics across replications. In addition, for both options, the following are also provided:
predweights 
List of predicted classification weights for each Parafac or Parafac2 model and for each replication. 
train.weights 
List of lists of training weights for each Parafac or Parafac2 model and for each replication. 
opt.tune 
List of optimal tuning parameters for classification methods for each Parafac or Parafac2 model and for each replication. 
mean.opt.tune 
Mean across all replications of optimal tuning parameters for classification methods for each Parafac or Parafac2 model. 
X 
Threeway or fourway data array or list used in argument 
nfac 
Number of components used to fit each Parafac or Parafac2 model. 
model 
Character designating the Parafac model that was used, either 
method 
Classification methods used. 
const 
Constraints used in fitting Parafac or Parafac2 models. 
Note
If argument cmode
is not null, input array x
is reshaped with function aperm
such that the cmode
dimension of x
is ordered last. Estimated mode A and B (and mode C for a fourway array) weights that are outputted as Aweights
and Bweights
(and Cweights
) reflect this permutation. For example, if x
is a fourway array and cmode = 2
, the original input modes 1, 2, 3, and 4 will correspond to output modes 1, 3, 4, 2. Here, output A = input 1; B = 3, and C = 4 (i.e., the second mode specified by cmode
has been moved to the D mode/last mode). For model = "parafac2"
, classification mode is assumed to be the last mode (i.e., mode C for threeway array and mode D for fourway array).
In addition, note that the following combination of arguments will give an error: nfac = 1, family = "multinomial", method = "PLR"
. The issue arises from providing glmnet::cv.glmnet
input x
a matrix with a single column. The issue is resolved for family = "binomial"
because a column of 0s is appended to the single column, but this solution does not appear to work for the multiclass case. As such, this combination of arguments is not currently allowed. This issue will be resolved in a future update.
Author(s)
Matthew A. Snodgress <snodg031@umn.edu>
References
Breiman, L. (2001). Random forests. Machine Learning, 45(1), 532.
Chen, T., He, T., Benesty, M., Khotilovich, V., Tang, Y., Cho, H., Chen, K., Mitchell, R., Cano, I., Zhou, T., Li, M., Xie, J., Lin, M., Geng, Y., Li, Y., Yuan, J. (2024). xgboost: Extreme gradient boosting. R Package Version 1.7.7.1.
Cortes, C. and Vapnik, V. (1995). Supportvector networks. Machine Learning, 20(3), 273297.
Friedman, J. H. (2001). Greedy function approximation: a gradient boosting machine. Annals of Statistics, 29(5), 11891232.
Friedman, J. H. (1989). Regularized discriminant analysis. Journal of the American Statistical Association, 84(405), 165175.
Friedman, J. Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software, 33(1), 122.
Guo, Y., Hastie, T., and Tibshirani, R. (2007). Regularized linear discriminant analysis and its application in microarrays. Biostatistics, 8(1), 86100.
Guo Y., Hastie T., and Tibshirani, R. (2023). rda: Shrunken centroids regularized discriminant analysis. R Package Version 1.21.
Harshman, R. (1970). Foundations of the PARAFAC procedure: Models and conditions for an "explanatory" multimodal factor analysis. UCLA Working Papers in Phonetics, 16, 184.
Harshman, R. (1972). PARAFAC2: Mathematical and technical notes. UCLA Working Papers in Phonetics, 22, 3044.
Harshman, R. and Lundy, M. (1994). PARAFAC: Parallel factor analysis. Computational Statistics and Data Analysis, 18, 3972.
Helwig, N. (2017). Estimating latent trends in multivariate longitudinal data via Parafac2 with functional and structural constraints. Biometrical Journal, 59(4), 783803.
Helwig, N. (2019). multiway: Component models for multiway data. R Package Version 1.06.
Liaw, A. and Wiener, M. (2002). Classification and regression by randomForest. R News 2(3), 18–22.
Meyer, D., Dimitriadou, E., Hornik, K., Weingessel, A., and Leisch, F. (2023). e1071: Misc functions of the Department of Statistics, Probability Theory Group (Formerly: E1071), TU Wien. R Package Version 1.713.
Ripley, B. (1994). Neural networks and related methods for classification. Journal of the Royal Statistical Society: Series B (Methodological), 56(3), 409437.
Venables, W. and Ripley, B. (2002). Modern applied statistics with S. Fourth Edition. Springer, New York. ISBN 0387954570.
Zou, H. and Hastie, T. (2005). Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(2), 301320.
Examples
########## Parafac2 example with 4way array and multiclass response ##########
# set seed and specify dimensions of a fourway tensor
set.seed(5)
mydim < c(10, 11, 12, 100)
nf < 3
# create correlation matrix between response and fourth mode's weights
rho.dd < .35
rho.dy < .75
cormat.values < c(1, rho.dd, rho.dd, rho.dy, rho.dd, 1, rho.dd, rho.dy,
rho.dd, rho.dd, 1, rho.dy, rho.dy, rho.dy, rho.dy, 1)
cormat < matrix(cormat.values, nrow = (nf + 1), ncol = (nf + 1))
# sample from a multivariate normal with specified correlation structure
ymean < Dmean < 2
mu < as.matrix(c(Dmean, Dmean, Dmean, ymean))
eidecomp < eigen(cormat, symmetric = TRUE)
L.sqrt < diag(eidecomp$values^0.5)
cormat.sqrt < eidecomp$vectors %*% L.sqrt %*% t(eidecomp$vectors)
Z < matrix(rnorm(mydim[4] * (nf + 1)), nrow = mydim[4], ncol = (nf + 1))
Xw < rep(1, mydim[4]) %*% t(mu) + Z %*% cormat.sqrt
Dmat < Xw[, 1:nf]
# create a random fourway data tensor with D weights related to a response
Bmat < matrix(runif(mydim[2] * nf), nrow = mydim[2], ncol = nf)
Cmat < matrix(runif(mydim[3] * nf), nrow = mydim[3], ncol = nf)
nDd < rep(c(10, 12, 14), length.out = mydim[4])
Gmat < matrix(rnorm(nf * nf), nrow = nf)
Amat < vector("list", mydim[4])
X < Xmat < Emat < Amat
for (Dd in 1:mydim[4]) {
Amat[[Dd]] < matrix(nf * rnorm(nDd[Dd]), nrow = nDd[Dd], ncol = nf)
Amat[[Dd]] < svd(Amat[[Dd]], nv = 0)$u %*% Gmat
leftMat < Amat[[Dd]] %*% diag(Dmat[Dd,])
Xmat[[Dd]] < array(tcrossprod(leftMat, krprod(Cmat, Bmat)),
dim = c(nDd[Dd], mydim[2], mydim[3]))
Emat[[Dd]] < array(rnorm(nDd[Dd] * mydim[2] * mydim[3]),
dim = c(nDd[Dd], mydim[2], mydim[3]))
X[[Dd]] < Xmat[[Dd]] + Emat[[Dd]]
}
# create a multiclass response
stor < matrix(rep(1, nrow(Xw)), nrow = nrow(Xw))
stor[which(Xw[, (nf + 1)] < (ymean  0.4 * sd(Xw[, (nf + 1)])))] < 2
stor[which(Xw[, (nf + 1)] > (ymean + 0.4 * sd(Xw[, (nf + 1)])))] < 0
y < factor(stor)
# initialize
alpha < seq(0, 1, length = 2)
gamma < c(0, 1)
cost < c(0.1, 5)
ntree < c(200, 300)
nodesize < c(1, 2)
size < c(1, 2)
decay < c(0, 1)
rda.alpha < seq(0.1, 0.9, length = 2)
delta < c(0.1, 2)
eta < c(0.3, 0.7)
max.depth < c(1, 2)
subsample < c(0.75)
nrounds < c(100)
method < c("PLR", "SVM", "RF", "NN", "RDA", "GBM")
family < "multinomial"
parameters < list(alpha = alpha, gamma = gamma, cost = cost, ntree = ntree,
nodesize = nodesize, size = size, decay = decay,
rda.alpha = rda.alpha, delta = delta, eta = eta,
max.depth = max.depth, subsample = subsample,
nrounds = nrounds)
model < "parafac2"
nfolds < 3
nstart < 3
# constrain first mode weights to be orthogonal, fourth mode to be nonnegative
const < c("orthog", "uncons", "uncons", "nonneg")
# fit Parafac2 model and use fourth mode weights to tune classification
# methods, to predict class labels, and to return classificaiton
# performance measures pooled across multiple traintest splits
output < cpfa(x = X, y = y, model = model, nfac = nf, nrep = 2, ratio = 0.8,
nfolds = nfolds, method = method, family = family,
parameters = parameters, type.out = "descriptives",
seeds = NULL, plot.out = TRUE, parallel = FALSE, const = const,
nstart = nstart)
# print performance measure means across traintest splits
output$descriptive$mean