# Multiway Covariates Regression 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. ```r 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) ``` ## Arguments - `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`? ## Details 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. ## Returns - **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"`. ## References Smilde, A. K., & Kiers, H. A. L. (1999). Multiway covariates regression models, **Journal of Chemometrics, 13**, 31-48. ## Author(s) Nathaniel E. Helwig ## Note 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. ## See Also 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. ## Examples ```r ########## 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) ```

mcr function