mcr: Multiway Covariates Regression
In multiway: Component Models for Multi-Way Data

Description Usage Arguments Details Value Note Author(s) References See Also Examples

Fits Smilde and Kiers's Multiway Covariates Regression (MCR) model to connect a 3-way predictor array and a 2-way response array that share a common mode. Parameters are estimated via alternating least squares with optional constraints.

mcr(X, Y, nfac = 1, alpha = 0.5, nstart = 10, 
    model = c("parafac", "parafac2", "tucker"),
    const = NULL, control = NULL, weights = NULL,
    Afixed = NULL, Bfixed = NULL, Cfixed = NULL, Dfixed = NULL,
    Astart = NULL, Bstart = NULL, Cstart = NULL, Dstart = NULL,
    Astruc = NULL, Bstruc = NULL, Cstruc = NULL, Dstruc = NULL,
    Amodes = NULL, Bmodes = NULL, Cmodes = NULL, Dmodes = NULL,
    maxit = 500, ctol = 1e-4, parallel = FALSE, cl = NULL, 
    output = c("best", "all"), verbose = TRUE, backfit = FALSE)

`X`	Three-way predictor array with `dim = c(I,J,K)`.
`Y`	Two-way response array with `dim = c(K,L)`.
`nfac`	Number of factors.
`alpha`	Tuning parameter between 0 and 1.
`nstart`	Number of random starts.
`model`	Model for `X`. Defaults to "parafac".
`const`	Character vector of length 4 giving the constraints for `A`, `B`, `C`, and `D` (defaults to unconstrained). See `const` for the 24 available options. Ignored if `model = "tucker"`.
`control`	List of parameters controlling options for smoothness constraints. This is passed to `const.control`, which describes the available options.
`weights`	Vector of length `K` giving non-negative weights for fitting via weighted least squares. Defaults to vector of ones.
`Afixed`	Used to fit model with fixed Mode A weights.
`Bfixed`	Used to fit model with fixed Mode B weights.
`Cfixed`	Used to fit model with fixed Mode C weights.
`Dfixed`	Used to fit model with fixed Mode D weights.
`Astart`	Starting Mode A weights. Default uses random weights.
`Bstart`	Starting Mode B weights. Default uses random weights.
`Cstart`	Starting Mode C weights. Default uses random weights.
`Dstart`	Starting Mode D weights. Default uses random weights.
`Astruc`	Structure constraints for Mode A weights. See Note.
`Bstruc`	Structure constraints for Mode B weights. See Note.
`Cstruc`	Structure constraints for Mode C weights. Ignored.
`Dstruc`	Structure constraints for Mode D weights. See Note.
`Amodes`	Mode ranges for Mode A weights (for unimodality constraints). See Note.
`Bmodes`	Mode ranges for Mode B weights (for unimodality constraints). See Note.
`Cmodes`	Mode ranges for Mode C weights (for unimodality constraints). Ignored.
`Dmodes`	Mode ranges for Mode D weights (for unimodality constraints). See Note.
`maxit`	Maximum number of iterations.
`ctol`	Convergence tolerance (R^2 change).
`parallel`	Logical indicating if `parLapply` should be used. See Examples.
`cl`	Cluster created by `makeCluster`. Only used when `parallel=TRUE`.
`output`	Output the best solution (default) or output all `nstart` solutions.
`verbose`	If `TRUE`, fitting progress is printed via `txtProgressBar`. Ignored if `parallel=TRUE`.
`backfit`	Should backfitting algorithm be used for `cmls`?

Given a predictor array X = array(x, dim=c(I,J,K)) and a response matrix Y = matrix(y, nrow=K, ncol=L), the multiway covariates regression (MCR) model assumes a tensor model for X and a bilinear model for Y, which are linked through a common C weight matrix. For example, using the Parafac model for X, the MCR model has the form

X[i,j,k] = sum A[i,r]*B[j,r]*C[k,r] + Ex[i,j,k]

and

Y[k,l] = sum C[k,r]*D[l,r] + Ey[k,l]

Parameter matrices are estimated by minimizing the loss function

LOSS = alpha * (SSE.X / SSX) + (1 - alpha) * (SSE.Y / SSY)

where

SSE.X = sum((X - Xhat)^2)

SSE.Y = sum((Y - Yhat)^2)

SSX = sum(X^2)

SSY = sum(Y^2)

When weights are input, SSE.X, SSE.Y, SSX, and SSY are replaced by the corresponding weighted versions.

`A`	Predictor A weight matrix.
`B`	Predictor B weight matrix.
`C`	Common C weight matrix.
`D`	Response D weight matrix.
`W`	Coefficients. See Note.
`LOSS`	Value of `LOSS` function.
`SSE`	Sum of Squared Errors for `X` and `Y`.
`Rsq`	R-squared value for `X` and `Y`.
`iter`	Number of iterations.
`cflag`	Convergence flag. See Note.
`model`	See argument `model`.
`const`	See argument `const`.
`control`	See argument `control`.
`weights`	See argument `weights`.
`alpha`	See argument `alpha`.
`fixed`	Logical vector indicating whether 'fixed' weights were used for each matrix.
`struc`	Logical vector indicating whether 'struc' constraints were used for each matrix.
`Phi`	Mode A crossproduct matrix. Only if `model = "parafac2"`.
`G`	Core array. Only if `model = "tucker"`.

When model = "parafac2", the arguments Afixed, Astart, and Astruc are treated as the arguments Gfixed, Gstart, and Gstruc from the parafac2 function.

Structure constraints should be specified with a matrix of logicals (TRUE/FALSE), such that FALSE elements indicate a weight should be constrained to be zero. Default uses unstructured weights, i.e., a matrix of all TRUE values. Structure constraints are ignored if model = "tucker".

When using unimodal constraints, the *modes inputs can be used to specify the mode search range for each factor. These inputs should be matrices with dimension c(2,nfac) where the first row gives the minimum mode value and the second row gives the maximum mode value (with respect to the indicies of the corresponding weight matrix).

C = Xc %*% W where Xc = matrix(aperm(X,c(3,1,2)),K)

Output cflag gives convergence information: cflag = 0 if algorithm converged normally, cflag = 1 if maximum iteration limit was reached before convergence, and cflag = 2 if algorithm terminated abnormally due to a problem with the constraints.

Nathaniel E. Helwig <helwig@umn.edu>

Smilde, A. K., & Kiers, H. A. L. (1999). Multiway covariates regression models, Journal of Chemometrics, 13, 31-48.

The fitted.mcr function creates the model-implied fitted values from a fit "mcr" object.

The resign.mcr function can be used to resign factors from a fit "mcr" object.

The rescale.mcr function can be used to rescale factors from a fit "mcr" object.

The reorder.mcr function can be used to reorder factors from a fit "mcr" object.

The cmls function (from CMLS package) is called as a part of the alternating least squares algorithm.

See parafac, parafac2, and tucker for more information about the Parafac, Parafac2, and Tucker models.

##########   multiway covariates regression   ##########

# create random data array with Parafac structure
set.seed(3)
mydim <- c(10, 20, 100)
nf <- 2
Amat <- matrix(rnorm(mydim[1]*nf), mydim[1], nf)
Bmat <- matrix(rnorm(mydim[2]*nf), mydim[2], nf)
Cmat <- matrix(rnorm(mydim[3]*nf), mydim[3], nf)
Xmat <- tcrossprod(Amat, krprod(Cmat, Bmat))
Xmat <- array(Xmat, dim = mydim)
EX <- array(rnorm(prod(mydim)), dim = mydim)
EX <- nscale(EX, 0, ssnew = sumsq(Xmat))   # SNR = 1
X <- Xmat + EX

# create response array
ydim <- c(mydim[3], 4)
Dmat <- matrix(rnorm(ydim[2]*nf), ydim[2], nf)
Ymat <- tcrossprod(Cmat, Dmat)
EY <- array(rnorm(prod(ydim)), dim = ydim)
EY <- nscale(EY, 0, ssnew = sumsq(Ymat))   # SNR = 1
Y <- Ymat + EY

# fit MCR model
mcr(X, Y, nfac = nf, nstart = 1)
mcr(X, Y, nfac = nf, nstart = 1, model = "parafac2")
mcr(X, Y, nfac = nf, nstart = 1, model = "tucker")



## Not run: 

##########   parallel computation   ##########

# create random data array with Parafac structure
set.seed(3)
mydim <- c(10, 20, 100)
nf <- 2
Amat <- matrix(rnorm(mydim[1]*nf), mydim[1], nf)
Bmat <- matrix(rnorm(mydim[2]*nf), mydim[2], nf)
Cmat <- matrix(rnorm(mydim[3]*nf), mydim[3], nf)
Xmat <- tcrossprod(Amat, krprod(Cmat, Bmat))
Xmat <- array(Xmat, dim = mydim)
EX <- array(rnorm(prod(mydim)), dim = mydim)
EX <- nscale(EX, 0, ssnew = sumsq(Xmat))   # SNR = 1
X <- Xmat + EX

# create response array
ydim <- c(mydim[3], 4)
Dmat <- matrix(rnorm(ydim[2]*nf), ydim[2], nf)
Ymat <- tcrossprod(Cmat, Dmat)
EY <- array(rnorm(prod(ydim)), dim = ydim)
EY <- nscale(EY, 0, ssnew = sumsq(Ymat))   # SNR = 1
Y <- Ymat + EY

# fit MCR-Parafac model (10 random starts -- sequential computation)
set.seed(1)
system.time({mod <- mcr(X, Y, nfac = nf)})
mod

# fit MCR-Parafac model (10 random starts -- parallel computation)
cl <- makeCluster(detectCores())
ce <- clusterEvalQ(cl, library(multiway))
clusterSetRNGStream(cl, 1)
system.time({mod <- mcr(X, Y, nfac = nf, parallel = TRUE, cl = cl)})
mod
stopCluster(cl)


## End(Not run)