Type:	Package
Title:	Non-Negative Matrix Factorization with Kernel Covariates
Version:	0.6.7
Date:	2026-03-30
Maintainer:	Kenichi Satoh <kenichi-satoh@biwako.shiga-u.ac.jp>
URL:	https://github.com/ksatohds/nmfkc, https://ksatohds.github.io/nmfkc/
BugReports:	https://github.com/ksatohds/nmfkc/issues
Description:	Performs Non-negative Matrix Factorization (NMF) with Kernel Covariates. Given an observation matrix and kernel covariates, it optimizes both a basis matrix and a parameter matrix. Notably, if the kernel matrix is an identity matrix, the method simplifies to standard NMF. Also provides NMF with Random Effects (NMF-RE) via nmfre(), which estimates a mixed-effects model combining covariate-driven scores with unit-specific random effects together with wild bootstrap inference, and NMF-based Structural Equation Modeling (NMF-SEM) via nmf.sem(), which fits a two-block input-output model for blind source separation and path analysis. References: Satoh (2025) <doi:10.48550/arXiv.2403.05359>; Satoh (2025) <doi:10.48550/arXiv.2510.10375>; Satoh (2025) <doi:10.48550/arXiv.2512.18250>; Satoh (2026) <doi:10.48550/arXiv.2603.01468>; Satoh (2026) <doi:10.1007/s42081-025-00314-0>.
License:	MIT + file LICENSE
Imports:	stats, graphics, utils, grDevices
Encoding:	UTF-8
Language:	en-US
RoxygenNote:	7.2.3
ByteCompile:	true
VignetteBuilder:	knitr
Suggests:	knitr, rmarkdown, testthat (≥ 3.0.0), mclust, palmerpenguins, quanteda, vars, DiagrammeR, MASS, nlme, lavaan
Config/testthat/edition:	3
NeedsCompilation:	no
Packaged:	2026-04-08 21:27:06 UTC; ksato
Author:	Kenichi Satoh [aut, cre]
Repository:	CRAN
Date/Publication:	2026-04-15 13:00:14 UTC

Extract coefficients from NMF models

Description

Returns the coefficients data frame from a fitted NMF model that has been passed through an inference function (nmfkc.inference, nmfae.inference, nmfre.inference).

If inference has not been run, returns the parameter matrix C (\Theta) directly.

For nmf.sem objects, returns C_2 (exogenous block) as fallback.

Usage

## S3 method for class 'nmf'
coef(object, ...)

## S3 method for class 'nmf.sem'
coef(object, ...)

Arguments

Value

A data frame of coefficients (if inference was performed), or the parameter matrix C.

See Also

nmfkc.inference, nmfae.inference, nmfre.inference, nmf.sem.inference

Examples

Y <- matrix(cars$dist, nrow = 1)
A <- rbind(1, cars$speed)
result <- nmfkc(Y, A, rank = 1)
coef(result)  # returns C matrix

result2 <- nmfkc.inference(result, Y, A)
coef(result2)  # returns coefficients data frame

Extract fitted values from NMF models

Description

Returns the reconstructed matrix \hat{Y} = X B from a fitted NMF model.

For nmf.sem objects, returns the equilibrium prediction \hat{Y}_1 = M_{model} Y_2 if available. Supply Y1 and Y2 to get the direct reconstruction X (C_1 Y_1 + C_2 Y_2) instead.

Usage

## S3 method for class 'nmf'
fitted(object, ...)

## S3 method for class 'nmfae'
fitted(object, ...)

## S3 method for class 'nmf.sem'
fitted(object, ...)

Arguments

Value

The fitted matrix X B.

See Also

nmfkc, nmfae, nmfre, nmf.sem, residuals.nmf

Examples

result <- nmfkc(matrix(runif(50), 5, 10), rank = 2)
fitted(result)

NMF-SEM Main Estimation Algorithm

Description

Fits the NMF-SEM model

Y_1 \approx X \bigl( \Theta_1 Y_1 + \Theta_2 Y_2 \bigr)

under non-negativity constraints with orthogonality and sparsity regularization. The function returns the estimated latent factors, structural coefficient matrices, and the implied equilibrium (input–output) mapping.

At equilibrium, the model can be written as

Y_1 \approx (I - X \Theta_1)^{-1} X \Theta_2 Y_2 \equiv M_{\mathrm{model}} Y_2,

where M_{\mathrm{model}} = (I - X \Theta_1)^{-1} X \Theta_2 is a Leontief-type cumulative-effect operator in latent space.

Internally, the latent feedback and exogenous loading matrices are stored as C1 and C2, corresponding to \Theta_1 and \Theta_2, respectively.

Usage

nmf.sem(
  Y1,
  Y2,
  rank = NULL,
  X.init = NULL,
  X.L2.ortho = 100,
  C1.L1 = 1,
  C2.L1 = 0.1,
  epsilon = 1e-06,
  maxit = 20000,
  seed = 123,
  ...
)

Arguments

Value

A list with components:

References

Satoh, K. (2025). Applying non-negative matrix factorization with covariates to structural equation modeling for blind input-output analysis. arXiv:2512.18250. https://arxiv.org/abs/2512.18250

See Also

nmf.sem.inference, nmf.sem.cv, nmf.sem.split, nmf.sem.DOT, summary.nmf.sem

Examples

# Simple NMF-SEM with iris data (non-negative)
Y <- t(iris[, -5])
Y1 <- Y[1:2, ]  # Sepal
Y2 <- Y[3:4, ]  # Petal
result <- nmf.sem(Y1, Y2, rank = 2, maxit = 500)
result$MAE

Generate a Graphviz DOT Diagram for an NMF-SEM Model

Description

Creates a Graphviz DOT script that visualizes the structural network estimated by nmf.sem. The resulting diagram displays:

• endogenous observed variables (Y_1),

• exogenous observed variables (Y_2),

• latent factors (F_1, ..., F_Q),

together with the non-negative path coefficients whose magnitudes exceed a user-specified threshold.

Directed edges represent estimated relationships:

• Y_2 \rightarrow F_q: entries of C2 (exogenous loadings),

• F_q \rightarrow Y_1: rows of X (factor-to-endogenous mappings),

• Y_1 \rightarrow F_q: entries of C1 (feedback paths).

Edge widths are scaled by coefficient magnitude, and nodes are placed in optional visual clusters. Only variables participating in edges above the threshold are displayed, while latent factors are always shown.

Usage

nmf.sem.DOT(
  result,
  weight_scale = 5,
  weight_scale_y2f = weight_scale,
  weight_scale_fy1 = weight_scale,
  weight_scale_feedback = weight_scale,
  threshold = 0.01,
  rankdir = "LR",
  fill = TRUE,
  cluster.box = c("normal", "faint", "invisible", "none"),
  cluster.labels = NULL,
  hide.isolated = TRUE,
  sig.level = 0.1
)

Arguments

Value

A character string representing a valid Graphviz DOT script.

See Also

nmf.sem, nmf.sem.inference, plot.nmfkc.DOT

Examples

Y <- t(iris[, -5])
Y1 <- Y[1:2, ]
Y2 <- Y[3:4, ]
result <- nmf.sem(Y1, Y2, rank = 2, maxit = 500)
dot <- nmf.sem.DOT(result)
cat(dot)

Cross-Validation for NMF-SEM

Description

Performs K-fold cross-validation to evaluate the equilibrium mapping of the NMF-SEM model.

For each fold, nmf.sem is fitted on the training samples, yielding an equilibrium mapping \hat Y_1 = M_{\mathrm{model}} Y_2. The held-out endogenous variables Y_1 are then predicted from Y_2 using this mapping, and the mean absolute error (MAE) over all entries in the test block is computed. The returned value is the average MAE across folds.

This implements the hyperparameter selection strategy described in the paper: hyperparameters are chosen by predictive cross-validation rather than direct inspection of the internal structural matrices.

Usage

nmf.sem.cv(
  Y1,
  Y2,
  rank = NULL,
  X.init = NULL,
  X.L2.ortho = 100,
  C1.L1 = 1,
  C2.L1 = 0.1,
  epsilon = 1e-06,
  maxit = 20000,
  ...
)

Arguments

Value

A numeric scalar: mean MAE across CV folds.

See Also

nmf.sem

Examples

Y <- t(iris[, -5])
Y1 <- Y[1:2, ]
Y2 <- Y[3:4, ]
mae <- nmf.sem.cv(Y1, Y2, rank = 2, maxit = 500, nfolds = 3)
mae

Statistical inference for the exogenous parameter matrix C2

Description

nmf.sem.inference performs statistical inference on the exogenous parameter matrix C_2 from a fitted nmf.sem model, conditional on the estimated basis matrix \hat{X} and the endogenous parameter matrix \hat{C}_1.

Under the working model R = Y_1 - X C_1 Y_1 \approx X C_2 Y_2 + \varepsilon, inference on C_2 is conducted via sandwich covariance estimation and one-step wild bootstrap with non-negative projection.

Usage

nmf.sem.inference(object, Y1, Y2, wild.bootstrap = TRUE, ...)

Arguments

Value

The input object with additional inference components:

References

Satoh, K. (2025). Applying non-negative matrix factorization with covariates to structural equation modeling for blind input-output analysis. arXiv:2512.18250. https://arxiv.org/abs/2512.18250

See Also

nmf.sem, nmf.sem.DOT

Examples

Y <- t(iris[, -5])
Y1 <- Y[1:2, ]; Y2 <- Y[3:4, ]
res <- nmf.sem(Y1, Y2, rank = 2)
res2 <- nmf.sem.inference(res, Y1, Y2)
res2$coefficients

Heuristic Variable Splitting for NMF-SEM

Description

Infers a heuristic partition of observed variables into exogenous (Y_2) and endogenous (Y_1) blocks for use in NMF-SEM. The method is based on positive-SEM logic, causal ordering, and optional sign alignment using the first principal component (PC1).

The procedure:

• internally standardizes variables (mean 0, sd 1),

• optionally flips signs so that most variables align positively with PC1,

• infers a causal ordering by repeatedly regressing each variable on the remaining ones and selecting the variable with the largest minimum standardized coefficient,

• determines an exogenous block by scanning the ordering from upstream and stopping at the first variable whose strongest parent coefficient exceeds threshold.

If n.exogenous is supplied, it overrides the automatic threshold rule.

Usage

nmf.sem.split(
  x,
  n.exogenous = NULL,
  threshold = 0.1,
  auto.flipped = TRUE,
  verbose = FALSE
)

Arguments

Value

A list with:

See Also

nmf.sem

Examples

# Infer exogenous/endogenous split from iris
sp <- nmf.sem.split(iris[, -5], n.exogenous = 2)
sp$endogenous.variables
sp$exogenous.variables

Three-Layer Non-negative Matrix Factorization (NMF-AE)

Description

nmfae fits a three-layer nonnegative matrix factorization model Y_1 \approx X_1 \Theta X_2 Y_2, where X_1 is a decoder basis (column sum 1), \Theta is a bottleneck parameter matrix, X_2 is an encoder basis (row sum 1), and Y_2 is the input matrix.

When Y2 = Y1, the model acts as a non-negative autoencoder. When Y1 != Y2, it acts as a heteroencoder.

Initialization uses a three-step NMF procedure via nmfkc: (1) nmfkc(Y1, rank=Q) to obtain X_1, (2) nmfkc(Y1, A=Y2, rank=Q) with fixed X_1 to obtain C = \Theta X_2, (3) nmfkc(Y2, rank=R) to factor C into \Theta and X_2.

Usage

nmfae(
  Y1,
  Y2 = Y1,
  rank = 2,
  rank.encoder = rank,
  epsilon = 1e-04,
  maxit = 5000,
  verbose = FALSE,
  ...
)

Arguments

Value

An object of class "nmfae", a list with components:

Lifecycle

This function is experimental. The interface may change in future versions.

Source

Satoh, K. (2025). Applying Non-negative Matrix Factorization with Covariates to Multivariate Time Series. Japanese Journal of Statistics and Data Science.

References

Lee, D. D. and Seung, H. S. (2001). Algorithms for Non-negative Matrix Factorization. Advances in Neural Information Processing Systems, 13.

Saha, S. et al. (2022). Hierarchical Deep Learning Neural Network (HiDeNN): An Artificial Intelligence (AI) Framework for Computational Science and Engineering. Computer Methods in Applied Mechanics and Engineering, 399.

See Also

nmfae.inference, predict.nmfae, nmfae.ecv, nmfae.DOT, nmfkc

Examples

# Autoencoder example
Y <- matrix(c(1,0,1,0, 0,1,0,1, 1,1,0,0), nrow=3, byrow=TRUE)
res <- nmfae(Y, rank=2, rank.encoder=2)
res$r.squared

# Heteroencoder example
Y1 <- matrix(c(1,0,0,1), nrow=2)
Y2 <- matrix(runif(8), nrow=4)
res2 <- nmfae(Y1, Y2, rank=2, rank.encoder=2)

DOT graph visualization for nmfae objects

Description

nmfae.DOT generates a DOT language string for visualizing the structure of a three-layer NMF model. Two graph types are supported: "XCX" shows encoder factors, \Theta, and decoder factors; "YXCXY" shows the full structure from Y_2 through X_2, \Theta, X_1 to Y_1.

Edge widths are proportional to matrix element values, and edges below threshold are omitted for clarity.

Usage

nmfae.DOT(
  result,
  type = c("XCX", "YXCXY"),
  threshold = 0.01,
  sig.level = 0.1,
  rankdir = "LR",
  fill = TRUE,
  weight_scale = 5,
  weight_scale_x1 = weight_scale,
  weight_scale_theta = weight_scale,
  weight_scale_x2 = weight_scale,
  Y1.label = NULL,
  X1.label = NULL,
  X2.label = NULL,
  Y2.label = NULL,
  Y1.title = "Output (Y1)",
  X1.title = "Decoder (X1)",
  X2.title = "Encoder (X2)",
  Y2.title = "Input (Y2)",
  hide.isolated = TRUE
)

Arguments

Value

A character string containing the DOT graph specification.

See Also

nmfae

Examples

set.seed(1)
Y <- matrix(runif(20), nrow = 4)
res <- nmfae(Y, rank = 2)
dot <- nmfae.DOT(res)

Sample-wise k-fold Cross-Validation for nmfae

Description

nmfae.cv performs k-fold cross-validation by splitting columns (samples) of Y_1 and Y_2 into div folds. For each fold, the model Y_1 \approx X_1 \Theta X_2 Y_2 is fitted on the training samples and predictive performance is evaluated on the held-out samples.

When Y2 is a kernel matrix created by nmfkc.kernel (detected via attributes), the symmetric kernel splitting convention is used: Y2[train, train] for training and Y2[train, test] for prediction.

Usage

nmfae.cv(Y1, Y2 = Y1, rank = 2, rank.encoder = rank, ...)

Arguments

Value

A list with components:

See Also

nmfae, nmfae.ecv, nmfae.kernel.beta.cv, nmfkc.cv

Examples

Y <- t(iris[1:30, 1:4])
res <- nmfae.cv(Y, rank = 2, rank.encoder = 2, nfolds = 5, maxit = 500)
res$sigma

Element-wise Cross-Validation for nmfae (Wold's CV)

Description

nmfae.ecv performs k-fold element-wise cross-validation by randomly holding out individual elements of Y_1, assigning them a weight of 0 via Y1.weights, and evaluating the reconstruction error on those held-out elements.

This method (also known as Wold's CV) is suitable for determining the optimal rank pair (Q, R) in three-layer NMF. Both rank and rank.encoder accept vector inputs. When rank.encoder = NULL (default), rank.encoder is set equal to rank and pairs are evaluated element-wise (i.e., (Q_1, R_1), (Q_2, R_2), \dots). When rank.encoder is explicitly specified, all combinations of rank and rank.encoder are evaluated via expand.grid.

Usage

nmfae.ecv(Y1, Y2 = Y1, rank = 1:2, rank.encoder = NULL, ...)

Arguments

Value

A list with components:

See Also

nmfae, nmfkc.ecv

Examples

Y <- t(iris[1:30, 1:4])
# Default: rank.encoder=NULL -> paired rank=rank.encoder
res <- nmfae.ecv(Y, rank = 1:3, nfolds = 3, maxit = 500)
res$sigma
# Explicit rank.encoder: full grid
res2 <- nmfae.ecv(Y, rank = 1:3, rank.encoder = 1:3, nfolds = 3, maxit = 500)
res2$sigma

Heatmap visualization of nmfae factor matrices

Description

nmfae.heatmap displays the three factor matrices X_1, \Theta, and X_2 as side-by-side heatmaps. This provides an alternative to DOT graph visualization, especially when Y_2 has many variables (e.g., kernel matrix).

Usage

nmfae.heatmap(
  x,
  Y1.label = NULL,
  X1.label = NULL,
  X2.label = NULL,
  Y2.label = NULL,
  palette = NULL,
  ...
)

Arguments

Value

Invisible NULL. Called for its side effect (plot).

See Also

nmfae, plot.nmfae, nmfae.DOT

Examples

set.seed(1)
Y <- matrix(runif(20), nrow = 4)
res <- nmfae(Y, rank = 2)
nmfae.heatmap(res)

Statistical Inference for NMF-AE Parameter Matrix

Description

Performs post-estimation inference for \Theta in the three-layer NMF model Y_1 \approx X_1 \Theta X_2 Y_2, conditional on (\hat{X}_1, \hat{X}_2). Uses sandwich covariance estimation and one-step wild bootstrap with non-negative projection.

Usage

nmfae.inference(object, Y1, Y2 = Y1, wild.bootstrap = TRUE, ...)

Arguments

Value

An object of class c("nmfae.inference", "nmfae"), inheriting all components from the input object, with additional inference components:

See Also

nmfae, summary.nmfae.inference

Examples

Y <- matrix(c(1,0,1,0, 0,1,0,1, 1,1,0,0), nrow=3, byrow=TRUE)
res <- nmfae(Y, rank=2, rank.encoder=2)
res2 <- nmfae.inference(res, Y)
summary(res2)

Optimize kernel beta for nmfae by cross-validation

Description

nmfae.kernel.beta.cv selects the optimal beta parameter of the kernel function by evaluating nmfae.cv for each candidate value. The kernel matrix A = K(U, V; \beta) replaces Y_2 in the three-layer NMF model.

When beta = NULL, candidate values are automatically generated via nmfkc.kernel.beta.nearest.med.

Usage

nmfae.kernel.beta.cv(
  Y1,
  rank = 2,
  rank.encoder = rank,
  U,
  V = NULL,
  beta = NULL,
  plot = TRUE,
  ...
)

Arguments

Value

A list with components:

See Also

nmfae.cv, nmfkc.kernel, nmfkc.kernel.beta.cv

Examples

Y <- matrix(cars$dist, nrow = 1)
U <- matrix(cars$speed, nrow = 1)
res <- nmfae.kernel.beta.cv(Y, rank = 1, rank.encoder = 1, U = U,
                             beta = c(0.01, 0.02, 0.05), nfolds = 5)
res$beta

Rename decoder and encoder bases

Description

Assigns user-specified names to the decoder (X1 columns) and encoder (X2 rows) bases of an nmfae object. The names propagate to \Theta, the coefficients table, and all downstream displays such as summary, nmfae.DOT, and nmfae.heatmap.

Usage

nmfae.rename(x, X1.colnames = NULL, X2.rownames = NULL)

Arguments

Value

A modified copy of x with updated names.

See Also

nmfae

Examples

set.seed(1)
Y <- matrix(runif(15), nrow = 3)
res <- nmfae(Y, rank = 2, rank.encoder = 2)
res <- nmfae.rename(res,
  X1.colnames = c("Basis1", "Basis2"),
  X2.rownames = c("Enc1", "Enc2"))
summary(res)

Optimize NMF with kernel covariates (Full Support for Missing Values)

Description

nmfkc fits a nonnegative matrix factorization with kernel covariates under the tri-factorization model Y \approx X C A = X B.

This function supports two major input modes:

1. Matrix Mode (Existing): nmfkc(Y=matrix, A=matrix, ...)

2. Formula Mode (New): nmfkc(formula=Y_vars ~ A_vars, data=df, rank=Q, ...)

The rank of the basis matrix can be specified using either the rank argument (preferred for formula mode) or the hidden Q argument (for backward compatibility).

Usage

nmfkc(
  Y,
  A = NULL,
  rank = NULL,
  data,
  epsilon = 1e-04,
  maxit = 5000,
  verbose = TRUE,
  ...
)

Arguments

Value

A list with components:

References

Satoh, K. (2024). Applying Non-negative Matrix Factorization with Covariates to the Longitudinal Data as Growth Curve Model. arXiv:2403.05359. https://arxiv.org/abs/2403.05359

Satoh, K. (2025). Applying non-negative matrix factorization with covariates to multivariate time series data as a vector autoregression model. Japanese Journal of Statistics and Data Science. arXiv:2501.17446. doi:10.1007/s42081-025-00314-0

Satoh, K. (2025). Applying non-negative matrix factorization with covariates to label matrix for classification. arXiv:2510.10375. https://arxiv.org/abs/2510.10375

Ding, C., Li, T., Peng, W., & Park, H. (2006). Orthogonal Nonnegative Matrix Tri-Factorizations for Clustering. In Proc. 12th ACM SIGKDD (pp. 126–135). doi:10.1145/1150402.1150420

See Also

nmfkc.cv, nmfkc.rank, nmfkc.kernel, nmfkc.ar, predict.nmfkc

Examples

# Example 1. Matrix Mode (Existing)
X <- cbind(c(1,0,1),c(0,1,0))
B <- cbind(c(1,0),c(0,1),c(1,1))
Y <- X %*% B
rownames(Y) <- paste0("P",1:nrow(Y))
colnames(Y) <- paste0("N",1:ncol(Y))
print(X); print(B); print(Y)
res <- nmfkc(Y,rank=2,epsilon=1e-6)
res$X
res$B

# Example 2. Formula Mode
set.seed(1)
dummy_data <- data.frame(Y1=rpois(10,5), Y2=rpois(10,10),
                         A1=abs(rnorm(10,5)), A2=abs(rnorm(10,3)))
res_f <- nmfkc(Y1 + Y2 ~ A1 + A2, data=dummy_data, rank=2)

# Example 3. Symmetric NMF (bi: Y ~ X X^T)
S <- matrix(c(3,0,2, 0,3,1, 2,1,2), nrow=3)
res_bi <- nmfkc(S, rank=2, Y.symmetric="bi")
res_bi$X   # basis matrix (no column normalization)
res_bi$XB  # reconstruction X %*% t(X)

# Example 4. Symmetric NMF (tri: Y ~ X C X^T)
res_tri <- nmfkc(S, rank=2, Y.symmetric="tri")
res_tri$C  # Q x Q cluster interaction matrix
res_tri$XB # reconstruction X %*% C %*% t(X)

Generate Graphviz DOT Scripts for NMF or NMF-with-Covariates Models

Description

Produces a Graphviz DOT script visualizing the structure of an NMF model (Y \approx X C A) or its simplified forms.

Supported visualization types:

• "YX" — Standard NMF view: latent factors X map to observations Y.

• "YA" — Direct regression view: covariates A map directly to Y using the combined coefficient matrix X C.

• "YXA" — Full tri-factorization: A \rightarrow C \rightarrow X \rightarrow Y.

Edge widths are scaled by coefficient magnitude, and nodes with no edges above the threshold are omitted from the visualization.

Usage

nmfkc.DOT(
  result,
  type = c("YX", "YA", "YXA"),
  threshold = 0.01,
  sig.level = 0.1,
  rankdir = "LR",
  fill = TRUE,
  weight_scale = 5,
  weight_scale_ax = weight_scale,
  weight_scale_xy = weight_scale,
  weight_scale_ay = weight_scale,
  Y.label = NULL,
  X.label = NULL,
  A.label = NULL,
  Y.title = "Observation (Y)",
  X.title = "Basis (X)",
  A.title = "Covariates (A)",
  hide.isolated = TRUE
)

Arguments

Value

A character string representing a Graphviz DOT script.

See Also

nmfkc, nmfae.DOT, nmf.sem.DOT, nmfkc.ar.DOT, plot.nmfkc.DOT

Examples

Y <- matrix(cars$dist, nrow = 1)
A <- rbind(1, cars$speed)
result <- nmfkc(Y, A, rank = 1)
dot <- nmfkc.DOT(result)
cat(dot)

Construct observation and covariate matrices for a vector autoregressive model

Description

nmfkc.ar generates the observation matrix and covariate matrix corresponding to a specified autoregressive lag order.

If the input Y is a ts object, its time properties are preserved in the "tsp_info" attribute, adjusted for the lag. Additionally, the column names of Y and A are set to the corresponding time points.

Usage

nmfkc.ar(Y, degree = 1, intercept = TRUE)

Arguments

Value

A list containing:

References

Satoh, K. (2025). Applying non-negative matrix factorization with covariates to multivariate time series data as a vector autoregression model. Japanese Journal of Statistics and Data Science. arXiv:2501.17446. doi:10.1007/s42081-025-00314-0

See Also

nmfkc, nmfkc.ar.degree.cv, nmfkc.ar.stationarity, nmfkc.ar.DOT

Examples

# Example using AirPassengers (ts object)
d <- AirPassengers
ar_data <- nmfkc.ar(d, degree = 2)
dim(ar_data$Y)
dim(ar_data$A)

# Example using matrix input
Y <- matrix(1:20, nrow = 2)
ar_data <- nmfkc.ar(Y, degree = 1)
ar_data$degree.max

Generate a Graphviz DOT Diagram for NMF-AR / NMF-VAR Models

Description

Produces a Graphviz DOT script for visualizing autoregressive NMF-with-covariates models constructed via nmfkc.ar + nmfkc.

The diagram displays three types of directed relationships:

• Lagged predictors: T_{t-k} \rightarrow X,

• Current latent factors: X \rightarrow T_t,

• Optional intercept effects: Const -> X.

Importantly, no direct edges from lagged variables to current outputs (T_{t-k} \rightarrow T_t) are drawn, in accordance with the NMF-AR formulation.

Each block of lagged variables is displayed in its own DOT subgraph (e.g., “T-1”, “T-2”, ...), while latent factor nodes and current-time outputs are arranged in separate clusters.

Usage

nmfkc.ar.DOT(
  result,
  degree = 1,
  intercept = any(colnames(result$C) == "(Intercept)"),
  threshold = 0.1,
  rankdir = "RL",
  fill = TRUE,
  weight_scale_xy = 5,
  weight_scale_lag = 5,
  weight_scale_int = 3,
  hide.isolated = TRUE
)

Arguments

Value

A character string representing a Graphviz DOT file.

See Also

nmfkc.ar, nmfkc, plot.nmfkc.DOT

Examples

d <- AirPassengers
ar_data <- nmfkc.ar(d, degree = 2)
result <- nmfkc(ar_data$Y, ar_data$A, rank = 1)
dot <- nmfkc.ar.DOT(result, degree = 2)
cat(dot)

Optimize lag order for the autoregressive model

Description

nmfkc.ar.degree.cv selects the optimal lag order for an autoregressive model by applying cross-validation over candidate degrees.

This function accepts both standard matrices (Variables x Time) and ts objects (Time x Variables). ts objects are automatically transposed internally.

Usage

nmfkc.ar.degree.cv(
  Y,
  rank = 1,
  degree = 1:2,
  intercept = TRUE,
  plot = TRUE,
  ...
)

Arguments

Value

A list with components:

See Also

nmfkc.ar, nmfkc.cv

Examples

# Example using ts object directly
d <- AirPassengers

# Selection of degree (using ts object)
# Note: Y is automatically transposed if it is a ts object
nmfkc.ar.degree.cv(Y=d, rank=1, degree=11:14)

Forecast future values for NMF-VAR model

Description

nmfkc.ar.predict computes multi-step-ahead forecasts for a fitted NMF-VAR model using recursive forecasting.

If the fitted model contains time series property information (from nmfkc.ar), the forecasted values will have appropriate time-based column names.

Usage

nmfkc.ar.predict(x, Y, degree = NULL, n.ahead = 1)

Arguments

Value

A list with components:

See Also

nmfkc, nmfkc.ar

Examples

# Forecast AirPassengers
d <- AirPassengers
ar_data <- nmfkc.ar(d, degree = 2)
result <- nmfkc(ar_data$Y, ar_data$A, rank = 1)
pred <- nmfkc.ar.predict(result, Y = matrix(d, nrow = 1), degree = 2, n.ahead = 3)
pred$pred

Check stationarity of an NMF-VAR model

Description

nmfkc.ar.stationarity assesses the dynamic stability of a VAR model by computing the spectral radius of its companion matrix. It returns both the spectral radius and a logical indicator of stationarity.

Usage

nmfkc.ar.stationarity(x)

Arguments

Value

A list with components:

See Also

nmfkc, nmfkc.ar

Examples

# Check stationarity of fitted AR model
d <- AirPassengers
ar_data <- nmfkc.ar(d, degree = 2)
result <- nmfkc(ar_data$Y, ar_data$A, rank = 1)
nmfkc.ar.stationarity(result)

Create a class (one-hot) matrix from a categorical vector

Description

nmfkc.class converts a categorical or factor vector into a class matrix (one-hot encoded representation), where each row corresponds to a category and each column corresponds to an observation.

Usage

nmfkc.class(x)

Arguments

Value

A binary matrix with one row per unique category and one column per observation. Each column has exactly one entry equal to 1, indicating the category of the observation.

See Also

nmfkc

Examples

# Example.
Y <- nmfkc.class(iris$Species)
Y[,1:6]

Compute model selection criteria for a fitted nmfkc model

Description

nmfkc.criterion computes information criteria (ICp, AIC, BIC), clustering quality measures (silhouette, CPCC, dist.cor), and soft-clustering statistics (B.prob entropy, max, sd) from a fitted nmfkc model.

This function can be called on a model that was fitted with detail = "fast" or detail = "minimal" to compute the full set of criteria afterwards.

Usage

nmfkc.criterion(object, Y, detail = c("full", "fast", "minimal"), ...)

Arguments

Value

A list with components:

r.squared

R-squared between Y and XB.

sigma

Residual standard deviation.

mae

Mean absolute error.

B.prob

Column-normalized coefficient matrix (soft-clustering probabilities).

B.cluster

Hard clustering labels (argmax of B.prob per column).

X.prob

Row-normalized basis matrix.

X.cluster

Hard clustering labels per row of X.

criterion

Named list: ICp, ICp1, ICp2, ICp3, AIC, BIC, B.prob.sd.min, B.prob.max.mean, B.prob.entropy.mean, silhouette, CPCC, dist.cor.

See Also

nmfkc, nmfkc.rank

Examples

Y <- t(iris[, -5])
res <- nmfkc(Y, rank = 3, detail = "fast")
crit <- nmfkc.criterion(res, Y)
crit$criterion$silhouette

Perform k-fold cross-validation for NMF with kernel covariates

Description

nmfkc.cv performs k-fold cross-validation for the tri-factorization model Y \approx X C A = X B, where

• Y(P,N) is the observation matrix,

• A(R,N) is the covariate (or kernel) matrix,

• X(P,Q) is the basis matrix,

• C(Q,R) is the parameter matrix, and

• B(Q,N) is the coefficient matrix (B = C A).

Given Y (and optionally A), X and C are fitted on each training split and predictive performance is evaluated on the held-out split.

Usage

nmfkc.cv(Y, A = NULL, rank = 2, data, ...)

Arguments

Value

A list with components:

objfunc

Mean loss per valid entry over all folds (MSE for method="EU").

sigma

Residual standard error (RMSE). Available only if method="EU"; on the same scale as Y.

objfunc.block

Loss for each fold.

block

Vector of fold indices (1, …, div) assigned to each column of Y.

See Also

nmfkc, nmfkc.kernel.beta.cv, nmfkc.ar.degree.cv

Examples

# Example 1 (with explicit covariates):
Y <- matrix(cars$dist, nrow = 1)
A <- rbind(1, cars$speed)
res <- nmfkc.cv(Y, A, rank = 1)
res$objfunc

# Example 2 (kernel A and beta sweep):
Y <- matrix(cars$dist, nrow = 1)
U <- matrix(c(5, 10, 15, 20, 25), nrow = 1)
V <- matrix(cars$speed, nrow = 1)
betas <- 25:35/1000
obj <- numeric(length(betas))
for (i in seq_along(betas)) {
  A <- nmfkc.kernel(U, V, beta = betas[i])
  obj[i] <- nmfkc.cv(Y, A, rank = 1, nfolds = 10)$objfunc
}
betas[which.min(obj)]

Denormalize a matrix from [0,1] back to its original scale

Description

nmfkc.denormalize rescales a matrix with values in [0,1] back to its original scale using the column-wise minima and maxima of a reference matrix.

Usage

nmfkc.denormalize(x, ref = x)

Arguments

Value

A numeric matrix with values transformed back to the original scale.

See Also

nmfkc.normalize

Examples

x <- nmfkc.normalize(iris[, -5])
x_recovered <- nmfkc.denormalize(x, iris[, -5])
apply(x_recovered - iris[, -5], 2, max)

Perform Element-wise Cross-Validation (Wold's CV)

Description

nmfkc.ecv performs k-fold cross-validation by randomly holding out individual elements of the data matrix (element-wise), assigning them a weight of 0 via Y.weights, and evaluating the reconstruction error on those held-out elements.

This method (also known as Wold's CV) is theoretically robust for determining the optimal rank (Q) in NMF. This function supports vector input for Q, allowing simultaneous evaluation of multiple ranks on the same folds.

When Y.symmetric = "bi" or "tri" is passed via ..., fold creation uses only the upper triangle (including the diagonal) to prevent information leakage through the symmetric entries Y_{ij} = Y_{ji}.

Usage

nmfkc.ecv(Y, A = NULL, rank = 1:3, data, ...)

Arguments

Value

A list with components:

See Also

nmfkc, nmfkc.cv

Examples

# Element-wise CV to select rank
Y <- t(iris[1:30, 1:4])
res <- nmfkc.ecv(Y, rank = 1:2, nfolds = 3)
res$objfunc

Statistical inference for the parameter matrix C (Theta)

Description

nmfkc.inference performs statistical inference on the parameter matrix C (\Theta) from a fitted nmfkc model, conditional on the estimated basis matrix \hat{X}.

Under the working model Y = X C A + \varepsilon where \varepsilon_{pn} \stackrel{iid}{\sim} N(0, \sigma^2), inference is conducted via sandwich covariance estimation and one-step wild bootstrap with non-negative projection.

Usage

nmfkc.inference(object, Y, A = NULL, wild.bootstrap = TRUE, ...)

Arguments

Value

An object of class c("nmfkc.inference", "nmfkc"), inheriting all components from the input object, with additional inference components:

References

Satoh, K. (2026). Wild Bootstrap Inference for Non-Negative Matrix Factorization with Random Effects. arXiv:2603.01468. https://arxiv.org/abs/2603.01468

See Also

nmfkc, summary.nmfkc.inference

Examples

Y <- matrix(cars$dist, nrow = 1)
A <- rbind(intercept = 1, speed = cars$speed)
result <- nmfkc(Y, A, rank = 1)
result2 <- nmfkc.inference(result, Y, A)
summary(result2)

Create a kernel matrix from covariates

Description

nmfkc.kernel constructs a kernel matrix from covariate matrices. It supports Gaussian, Exponential, Periodic, Linear, Normalized Linear, and Polynomial kernels.

Usage

nmfkc.kernel(
  U,
  V = NULL,
  kernel = c("Gaussian", "Exponential", "Periodic", "Linear", "NormalizedLinear",
    "Polynomial"),
  ...
)

Arguments

Value

Kernel matrix A(N,M).

Source

Satoh, K. (2024). Applying Non-negative Matrix Factorization with Covariates to the Longitudinal Data as Growth Curve Model. arXiv preprint arXiv:2403.05359. https://arxiv.org/abs/2403.05359

See Also

nmfkc.kernel.gaussian, nmfkc.cv

Examples

# Example.
Y <- matrix(cars$dist,nrow=1)
U <- matrix(c(5,10,15,20,25),nrow=1)
V <- matrix(cars$speed,nrow=1)
A <- nmfkc.kernel(U,V,beta=28/1000)
dim(A)
result <- nmfkc(Y,A,rank=1)
plot(as.vector(V),as.vector(Y))
lines(as.vector(V),as.vector(result$XB),col=2,lwd=2)

Optimize beta of the Gaussian kernel function by cross-validation

Description

nmfkc.kernel.beta.cv selects the optimal beta parameter of the kernel function by applying cross-validation over a set of candidate values.

Usage

nmfkc.kernel.beta.cv(Y, rank = 2, U, V = NULL, beta = NULL, plot = TRUE, ...)

Arguments

Value

A list with components:

See Also

nmfkc.kernel.gaussian, nmfkc.kernel.beta.nearest.med, nmfkc.kernel

Examples

# Example.
Y <- matrix(cars$dist,nrow=1)
U <- matrix(c(5,10,15,20,25),nrow=1)
V <- matrix(cars$speed,nrow=1)
nmfkc.kernel.beta.cv(Y,rank=1,U,V,beta=25:30/1000)
A <- nmfkc.kernel(U,V,beta=28/1000)
result <- nmfkc(Y,A,rank=1)
plot(as.vector(V),as.vector(Y))
lines(as.vector(V),as.vector(result$XB),col=2,lwd=2)

Estimate Gaussian/RBF kernel parameter beta from covariates (supports landmarks)

Description

Computes a data-driven reference scale for the Gaussian/RBF kernel from covariates using a robust "median nearest-neighbor (or nearest-landmark) distance" heuristic, and returns the corresponding kernel parameter \beta.

The Gaussian/RBF kernel is assumed to be written in the form

k(u,v) = \exp\{-\beta \|u-v\|^2\} = \exp\{-\|u-v\|^2/(2\sigma^2)\},

hence \beta = 1/(2\sigma^2). This function first estimates a typical distance scale \sigma_0 by the median of distances, then sets \beta_0 = 1/(2\sigma_0^2).

If Uk is NULL, \sigma_0 is estimated as the median of nearest-neighbor distances within U (excluding self-distance). If Uk is provided, \sigma_0 is estimated as the median of nearest-landmark distances from each sample in U to its closest landmark in Uk.

To control memory usage for large N (and M), distances are computed in blocks. Optionally, columns of U can be randomly subsampled via sample.size to reduce cost.

Usage

nmfkc.kernel.beta.nearest.med(
  U,
  Uk = NULL,
  block.size = 1000,
  block.size.Uk = 2000,
  sample.size = NULL,
  ...
)

Arguments

Details

Candidate grid: Along with beta, the function returns beta_candidates, a small logarithmic grid suitable for cross-validation.

In the landmark case (Uk provided), the grid is designed to be symmetric on the bandwidth scale \sigma around \sigma_0 over one decade:

\sigma = \sigma_0 \times 10^{t}, \quad t \in \{-1,-2/3,-1/3,0,1/3,2/3,1\}.

Using \beta = 1/(2\sigma^2), this corresponds to

\beta = \beta_0 \times 10^{-2t}.

When Uk is NULL, a shorter coarse grid may be returned (see Value).

Notes:

• When Uk is identical to U, the function detects this case and excludes self-distances (distance 0) to avoid \sigma_0=0.

• sample.size performs random subsampling without setting a seed. For reproducible results, set set.seed() before calling this function.

Value

A list with elements:

• beta: Estimated kernel parameter \beta_0 = 1/(2\sigma_0^2).

• beta_candidates: Numeric vector of candidate \beta values (logarithmic grid) intended for cross-validation.

• dist_median: The estimated distance scale \sigma_0 (median of nearest-neighbor or nearest-landmark distances).

• block.size.used: The effective block size(s) used. Either a scalar (no Uk) or a named vector c(U=..., Uk=...) when Uk is provided.

• sample.size.used: The number of columns of U actually used (after subsampling).

• uk_is_u: Logical flag indicating whether Uk was detected as identical to U (only returned when Uk is provided).

See Also

nmfkc.kernel.gaussian, nmfkc.kernel.beta.cv

Examples

# Basic (nearest-neighbor within U)
U <- matrix(runif(20), nrow = 2)
beta_info <- nmfkc.kernel.beta.nearest.med(U)
beta0 <- beta_info$beta
betas <- beta_info$beta_candidates

# With landmarks (nearest-landmark distances)
Uk <- matrix(runif(10), nrow = 2)

beta_info2 <- nmfkc.kernel.beta.nearest.med(U, Uk)

Create a Gaussian kernel matrix from covariates

Description

nmfkc.kernel.gaussian constructs a Gaussian (RBF) kernel matrix from covariate matrices. The kernel is defined as K(u,v) = \exp(-\beta \|u - v\|^2). When V contains NA values, two methods are available via na.method:

"pds"

Partial Distance Strategy. Computes the kernel using only observed (non-NA) rows, with beta adjusted by \beta_{adj} = \beta \times K / K_{obs} where K is the total number of rows and K_{obs} is the number of observed rows.

"egk"

Expected Gaussian Kernel (Mesquita et al., 2019). Uses a Gaussian Mixture Model (GMM) to estimate the conditional distribution of missing values given observed values, then computes the expected kernel value via a Gamma approximation. Requires gmm.means, gmm.sigmas, and gmm.weights passed through ....

Usage

nmfkc.kernel.gaussian(
  U,
  V = NULL,
  beta = 0.5,
  na.method = c("pds", "egk"),
  ...
)

Arguments

Value

Kernel matrix A(N,M).

Source

Mesquita, D., Gomes, J. P., & Rodrigues, L. R. (2019). Gaussian kernels for incomplete data. Applied Soft Computing, 77, 356–365.

See Also

nmfkc.kernel, nmfkc.kernel.beta.cv, nmfkc.kernel.beta.nearest.med

Examples

U <- matrix(c(5,10,15,20,25),nrow=1)
V <- matrix(1:25,nrow=1)
A <- nmfkc.kernel.gaussian(U,V,beta=28/1000)
dim(A)

# PDS example: V with NA in first row
U2 <- matrix(rnorm(20), nrow=2)
V2 <- matrix(rnorm(10), nrow=2)
V2[1, c(2,4)] <- NA
A2 <- nmfkc.kernel.gaussian(U2, V2, beta=0.5, na.method="pds")

Normalize a matrix to the range [0,1]

Description

nmfkc.normalize rescales the values of a matrix to lie between 0 and 1 using the column-wise minimum and maximum values of a reference matrix.

Usage

nmfkc.normalize(x, ref = x)

Arguments

Value

A matrix of the same dimensions as x, with each column rescaled to the [0,1] range.

See Also

nmfkc.denormalize

Examples

# Example.
x <- nmfkc.normalize(iris[,-5])
apply(x,2,range)

Rank selection diagnostics with graphical output

Description

nmfkc.rank provides diagnostic criteria for selecting the rank (Q) in NMF with kernel covariates. Several model selection measures are computed (e.g., R-squared, silhouette, CPCC, ARI), and results can be visualized in a plot.

By default (save.time = FALSE), this function also computes the Element-wise Cross-Validation error (Wold's CV Sigma) using nmfkc.ecv.

The plot explicitly marks the "BEST" rank based on two criteria:

1. Elbow Method (Red): Based on the curvature of the R-squared values (always computed if Q > 2).

2. Min RMSE (Blue): Based on the minimum Element-wise CV Sigma (only if detail="full").

Usage

nmfkc.rank(Y, A = NULL, rank = 1:2, detail = "full", plot = TRUE, data, ...)

Arguments

Value

A list containing:

References

Brunet, J.P., Tamayo, P., Golub, T.R., Mesirov, J.P. (2004). Metagenes and molecular pattern discovery using matrix factorization. Proc. Natl. Acad. Sci. USA, 101, 4164–4169. doi:10.1073/pnas.0308531101 Punera, K., & Ghosh, J. (2008). Consensus-based ensembles of soft clusterings. Applied Artificial Intelligence, 22(7–8), 780–810. doi:10.1080/08839510802170546

See Also

nmfkc, nmfkc.ecv

Examples

# Example.
Y <- t(iris[,-5])
# Full run (default)
nmfkc.rank(Y, rank=1:4)
# Fast run (skip ECV)
nmfkc.rank(Y, rank=1:4, detail="fast")

Plot Diagnostics: Original, Fitted, and Residual Matrices as Heatmaps

Description

This function generates a side-by-side plot of three heatmaps: the original observation matrix Y, the fitted matrix XB (from NMF), and the residual matrix E (Y - XB). This visualization aids in diagnosing whether the chosen rank Q is adequate by assessing if the residual matrix E appears to be random noise.

The axis labels (X-axis: Samples, Y-axis: Features) are integrated into the main title of each plot to maximize the plot area, reflecting the compact layout settings.

Usage

nmfkc.residual.plot(
  Y,
  result,
  fitted.palette = (grDevices::colorRampPalette(c("white", "orange", "red")))(256),
  residual.palette = (grDevices::colorRampPalette(c("blue", "white", "red")))(256),
  ...
)

Arguments

Value

NULL. The function generates a plot.

See Also

nmfkc, residuals.nmf

Examples

Y <- t(iris[1:30, 1:4])
result <- nmfkc(Y, rank = 2)
nmfkc.residual.plot(Y, result)

Non-negative Matrix Factorization with Random Effects

Description

Estimates the NMF-RE model

Y = X(\Theta A + U) + \mathcal{E}

where Y (P \times N) is a non-negative observation matrix, X (P \times Q) is a non-negative basis matrix learned from the data, \Theta (Q \times K) is a non-negative coefficient matrix capturing systematic covariate effects on latent scores, A (K \times N) is a covariate matrix, and U (Q \times N) is a random effects matrix capturing unit-specific deviations in the latent score space.

NMF-RE can be viewed as a mixed-effects latent-variable model defined on a reconstruction (mean) structure. The non-negativity constraint on X induces sparse, parts-based loadings, achieving measurement-side variable selection without an explicit sparsity penalty. Inference on \Theta provides covariate-side variable selection by identifying which covariates significantly affect which components.

Estimation alternates ridge-type BLUP-like closed-form updates for U with multiplicative non-negative updates for X and \Theta. The effective degrees of freedom consumed by U are monitored and a df-based cap can be enforced to prevent near-saturated fits.

When wild.bootstrap = TRUE, inference on \Theta is performed conditional on (\hat{X}, \hat{U}) via asymptotic linearization, a one-step Newton update, and a multiplier (wild) bootstrap, yielding standard errors, z-values, p-values, and confidence intervals without repeated constrained re-optimization.

Usage

nmfre(
  Y,
  A = NULL,
  rank = 2,
  df.rate = NULL,
  wild.bootstrap = TRUE,
  epsilon = 1e-05,
  maxit = 50000,
  ...
)

Arguments

Value

A list of class "nmfre" with components. The model is Y = X(\Theta A + U) + \mathcal{E}.

Core matrices

X

Basis matrix X (P \times Q), columns normalized to sum to 1.

C

Coefficient matrix \Theta (Q \times K).

U

Random effects matrix U (Q \times N).

Variance components

sigma2

Residual variance \hat{\sigma}^2.

tau2

Random effect variance \hat{\tau}^2.

lambda

Ridge penalty \lambda = \sigma^2 / \tau^2.

Convergence diagnostics

converged

Logical. Whether the algorithm converged.

stop.reason

Character string describing why iteration stopped.

iter

Number of iterations performed.

maxit

Maximum iterations setting used.

epsilon

Convergence tolerance used.

objfunc

Final objective function value \|Y - X(\Theta A + U)\|^2 + \lambda \|U\|^2.

rel.change.final

Final relative change in objective.

objfunc.iter

Numeric vector of objective values per iteration.

rss.trace

Numeric vector of \|Y - X(\Theta A + U)\|^2 per iteration.

Effective degrees of freedom (dfU) diagnostics

dfU

Final effective degrees of freedom \mathrm{df}_U = N \sum_q d_q / (d_q + \lambda), where d_q are eigenvalues of X'X.

dfU.cap

Upper bound imposed on \mathrm{df}_U.

dfU.cap.rate

Rate used to compute the cap.

dfU.cap.scan

Result of nmfre.dfU.scan, or NULL.

lambda.enforced

Final \lambda enforced to satisfy the cap.

dfU.hit.cap

Logical. Whether the cap was binding.

dfU.hit.iter

Iteration at which the cap first bound.

dfU.frac

\mathrm{df}_U / (NQ), fraction of maximum df.

dfU.cap.frac

\mathrm{df}_U^{\mathrm{cap}} / (NQ).

Fitted matrices

B

Fixed-effect scores \Theta A (Q \times N).

B.prob

Column-normalized probabilities from \max(\Theta A, 0).

B.blup

BLUP scores \Theta A + U (Q \times N).

B.blup.pos

Non-negative BLUP scores \max(\Theta A + U, 0) (Q \times N).

B.blup.prob

Column-normalized probabilities from \max(\Theta A + U, 0).

XB

Fitted values from fixed effects X \Theta A (P \times N).

XB.blup

Fitted values including random effects X(\Theta A + U) (P \times N).

Fit statistics

r.squared

Coefficient of determination R^2 for Y vs X(\Theta A + U) (BLUP).

r.squared.fixed

Coefficient of determination R^2 for Y vs X \Theta A (fixed effects only).

ICC

Trace-based Intraclass Correlation Coefficient. In the NMF-RE model, the conditional covariance of the n-th observation column is \mathrm{Var}(Y_n) = \tau^2 X X^\top + \sigma^2 I_P, a P \times P matrix. Unlike a standard random intercept model where the design matrix Z is a simple indicator (so the ICC reduces to \tau^2 / (\sigma^2 + \tau^2)), the basis matrix X plays the role of Z in a random slopes model, making the variance contribution of U depend on X. To obtain a scalar summary, we take the trace of each component:

\mathrm{ICC} = \frac{\tau^2 \, \mathrm{tr}(X^\top X)} {\tau^2 \, \mathrm{tr}(X^\top X) + \sigma^2 P}.

This equals the average (over P dimensions) proportion of per-column variance attributable to the random effects.

Inference on \Theta (wild bootstrap)

sigma2.used

\hat{\sigma}^2 used for inference.

C.vec.cov

Variance-covariance matrix for \mathrm{vec}(\Theta) (QK \times QK).

C.se

Standard error matrix for \Theta (Q \times K).

C.se.hess

Sandwich (Hessian-based) SE matrix for \Theta.

C.se.boot

Bootstrap SE matrix for \Theta.

coefficients

Data frame with columns Estimate, Std. Error, z value, Pr(>|z|), and confidence interval bounds for each element of \Theta.

C.ci.lower

Lower confidence interval matrix for \Theta (Q \times K).

C.ci.upper

Upper confidence interval matrix for \Theta (Q \times K).

C.boot.sd

Bootstrap standard deviation matrix for \Theta (Q \times K).

C.p.side

P-value sidedness used: "one.sided" or "two.sided".

References

Satoh, K. (2026). Wild Bootstrap Inference for Non-Negative Matrix Factorization with Random Effects. arXiv:2603.01468. https://arxiv.org/abs/2603.01468

See Also

nmfre.inference, nmfre.dfU.scan, nmfkc.DOT, summary.nmfre

Examples

# Example 1. cars data
Y <- matrix(cars$dist, nrow = 1)
A <- rbind(intercept = 1, speed = cars$speed)
res <- nmfre(Y, A, rank = 1, maxit = 5000)
summary(res)


# Example 2. Orthodont data (nlme)
if (requireNamespace("nlme", quietly = TRUE)) {
  Y <- matrix(nlme::Orthodont$distance, 4, 27)
  male <- ifelse(nlme::Orthodont$Sex[seq(1, 108, 4)] == "Male", 1, 0)
  A <- rbind(intercept = 1, male = male)

  # Scan dfU cap rates to choose an appropriate value
  nmfre.dfU.scan(1:10/10, Y, A, rank = 1)

  # Fit with chosen rate
  res <- nmfre(Y, A, rank = 1, df.rate = 0.2)
  summary(res)
}

Scan dfU cap rates for NMF-RE

Description

Fits the NMF-RE model across a range of dfU.cap.rate values and returns a diagnostic table showing the resulting effective degrees of freedom, variance components, and convergence diagnostics for each rate.

The dfU cap limits the effective degrees of freedom consumed by the random effects U. The cap is computed as rate * N * Q, where N is the number of observations and Q is the rank. A suitable rate is one where the final \mathrm{df}_U is below the cap (safeguard = TRUE) and the model has converged (converged = TRUE).

When called automatically by nmfre (i.e., dfU.cap.rate = NULL), the minimum rate satisfying both safeguard = TRUE and converged = TRUE is selected.

Usage

nmfre.dfU.scan(
  rates = (1:10)/10,
  Y,
  A,
  rank = NULL,
  X.init = NULL,
  C.init = NULL,
  U.init = NULL,
  print.trace = FALSE,
  ...
)

Arguments

Value

An object of class "nmfre.dfU.scan" with two components:

table

A data frame with the following columns:

rate

Cap rate used. The dfU cap is rate * N * Q.

dfU.cap

The dfU cap value (upper bound on effective degrees of freedom).

dfU

Final effective degrees of freedom for U at convergence.

safeguard

Logical. TRUE if the dfU cap is functioning as a safeguard (dfU / dfU.cap < 0.99): the cap prevents random-effects saturation without over-constraining U. FALSE if dfU is at or near the cap, indicating the cap is binding and the rate may be too small.

hit

Logical. TRUE if the cap was reached at least once during iteration, even if dfU later decreased below the cap.

converged

Logical. TRUE if the algorithm converged within the maximum number of iterations.

tau2

Final random effect variance \hat{\tau}^2.

sigma2

Final residual variance \hat{\sigma}^2.

ICC

Trace-based Intraclass Correlation Coefficient \tau^2 \, \mathrm{tr}(X^\top X) / (\tau^2 \, \mathrm{tr}(X^\top X) + \sigma^2 P). See nmfre for details.

cap.rate

Optimal cap rate selected automatically. If rows with safeguard = TRUE and hit = TRUE exist, the maximum rate among them is chosen (safeguard activated but giving U the most freedom). Otherwise, the minimum rate with safeguard = TRUE and hit = FALSE is chosen. NA if no suitable rate is found.

When printed, only the table is displayed. Access cap.rate directly from the returned object.

See Also

nmfre

Examples

# Example 1. cars data (small maxit for speed)
Y <- matrix(cars$dist, nrow = 1)
A <- rbind(intercept = 1, speed = cars$speed)
tab <- nmfre.dfU.scan(rates = c(0.1, 0.2), Y = Y, A = A, rank = 1, maxit = 1000)
print(tab)


# Example 2. Orthodont data (nlme)
if (requireNamespace("nlme", quietly = TRUE)) {
  Y <- matrix(nlme::Orthodont$distance, 4, 27)
  male <- ifelse(nlme::Orthodont$Sex[seq(1, 108, 4)] == "Male", 1, 0)
  A <- rbind(intercept = 1, male = male)
  nmfre.dfU.scan(1:10/10, Y, A, rank = 1)
}

Statistical inference for the coefficient matrix C from NMF-RE

Description

nmfre.inference performs statistical inference on the coefficient matrix C (\Theta) from a fitted nmfre model, conditional on the estimated basis matrix \hat{X} and random effects \hat{U}.

Under the working model Y^* = Y - X\hat{U} \approx X C A + \varepsilon, inference is conducted via sandwich covariance estimation and one-step wild bootstrap with non-negative projection.

The result is compatible with nmfkc.DOT for visualization (pass the result directly as x with type = "YXA").

Usage

nmfre.inference(object, Y, A = NULL, wild.bootstrap = TRUE, ...)

Arguments

Value

The input object with additional inference components:

References

Satoh, K. (2026). Wild Bootstrap Inference for Non-Negative Matrix Factorization with Random Effects. arXiv:2603.01468. https://arxiv.org/abs/2603.01468

See Also

nmfre, nmfkc.DOT, summary.nmfre

Examples

Y <- matrix(cars$dist, nrow = 1)
A <- rbind(intercept = 1, speed = cars$speed)
res <- nmfre(Y, A, rank = 1, wild.bootstrap = FALSE)
res2 <- nmfre.inference(res, Y, A)
res2$coefficients

plot.nmfae displays the convergence trajectory of the objective function across iterations. The title shows the achieved R^2.

Description

plot.nmfae displays the convergence trajectory of the objective function across iterations. The title shows the achieved R^2.

Usage

## S3 method for class 'nmfae'
plot(x, ...)

Arguments

Value

Invisible NULL. Called for its side effect (plot).

See Also

nmfae, nmfae.heatmap

Examples

set.seed(1)
Y <- matrix(runif(20), nrow = 4)
res <- nmfae(Y, rank = 2)
plot(res)

Plot method for nmfae.cv objects

Description

Displays a bar chart of per-fold cross-validation errors from nmfae.cv. The overall RMSE (sigma) is shown in the title.

Usage

## S3 method for class 'nmfae.cv'
plot(x, ...)

Arguments

Value

Invisible NULL. Called for its side effect (plot).

See Also

nmfae.cv

Plot method for nmfae.ecv objects

Description

Visualizes element-wise cross-validation results. When rank.encoder was NULL (paired), a line plot of sigma vs rank is drawn. When rank.encoder was explicitly specified (grid), a heatmap of sigma over the (rank, rank.encoder) grid is drawn.

Usage

## S3 method for class 'nmfae.ecv'
plot(x, ...)

Arguments

Value

Invisible NULL. Called for its side effect of producing a plot.

See Also

nmfae.ecv

Plot method for nmfae.kernel.beta.cv objects

Description

Displays the cross-validation objective function across candidate beta values (log scale). The optimal beta is highlighted in red.

Usage

## S3 method for class 'nmfae.kernel.beta.cv'
plot(x, ...)

Arguments

Value

Invisible NULL. Called for its side effect (plot).

See Also

nmfae.kernel.beta.cv

Plot method for objects of class nmfkc

Description

plot.nmfkc produces a diagnostic plot for the return value of nmfkc, showing the objective function across iterations.

Usage

## S3 method for class 'nmfkc'
plot(x, ...)

Arguments

Value

Called for its side effect (a plot). Returns NULL invisibly.

See Also

nmfkc, summary.nmfkc

Examples

Y <- matrix(cars$dist, nrow = 1)
A <- rbind(1, cars$speed)
result <- nmfkc(Y, A, rank = 1)
plot(result)

Plot method for nmfkc.DOT objects

Description

Renders a DOT graph string using DiagrammeR::grViz. If the DiagrammeR package is not installed, prints the DOT source to the console instead.

This method handles all DOT objects produced by the nmfkc package: nmfkc.DOT, nmfae.DOT, nmf.sem.DOT, and nmfkc.ar.DOT.

Usage

## S3 method for class 'nmfkc.DOT'
plot(x, ...)

Arguments

Value

Called for its side effect (rendering). Returns x invisibly.

See Also

nmfkc.DOT, nmfae.DOT, nmf.sem.DOT, nmfkc.ar.DOT

Plot convergence diagnostics for NMF models

Description

Plots the objective function value over iterations for nmfre and nmf.sem objects. (For nmfkc and nmfae, plot methods are defined in their respective source files.)

Usage

## S3 method for class 'nmfre'
plot(x, ...)

## S3 method for class 'nmf.sem'
plot(x, ...)

Arguments

Value

Invisible NULL.

See Also

nmfre, nmf.sem

Examples

set.seed(1)
Y <- matrix(runif(20), nrow = 4)
A <- diag(5)
res <- nmfre(Y, A, rank = 2, wild.bootstrap = FALSE)
plot(res)

Plot method for predict.nmfae objects

Description

For type = "response": if actual values Y_1 were stored, displays an observed-vs-predicted scatter plot with R^2 in the title. Otherwise, displays the predicted matrix as a heatmap.

For type = "class": if actual classes were stored, displays a confusion matrix heatmap with accuracy (ACC) in the title.

Usage

## S3 method for class 'predict.nmfae'
plot(x, ...)

Arguments

Value

Invisible NULL. Called for its side effect (plot).

See Also

predict.nmfae

Examples

set.seed(1)
Y <- matrix(runif(20), nrow = 4)
res <- nmfae(Y, rank = 2)
pred <- predict(res)
plot(pred)

Predict method for nmfae objects

Description

predict.nmfae computes fitted or predicted values from a three-layer NMF model. Without newY2, returns the in-sample fitted values X_1 \Theta X_2 Y_2. With newY2, computes out-of-sample predictions X_1 \Theta X_2 \cdot \mathrm{newY2}.

When type = "class", each column is classified to the row with the maximum predicted value (useful when Y_1 is a one-hot class matrix from nmfkc.class).

If Y1 (actual values) is provided, it is stored as an attribute so that plot.predict.nmfae can produce an observed-vs-predicted scatter plot (for type = "response") or a confusion matrix heatmap (for type = "class").

Usage

## S3 method for class 'nmfae'
predict(object, newY2 = NULL, Y1 = NULL, type = c("response", "class"), ...)

Arguments

Value

For type = "response": a matrix of class "predict.nmfae". For type = "class": a factor of class "predict.nmfae" with predicted class labels. If Y1 was provided, actual classes are stored in attr(result, "actual").

See Also

nmfae, plot.predict.nmfae, nmfkc.class

Examples

set.seed(1)
Y <- matrix(runif(20), nrow = 4)
res <- nmfae(Y, rank = 2)
pred <- predict(res)

Prediction method for objects of class nmfkc

Description

predict.nmfkc generates predictions from an object of class nmfkc, either using the fitted covariates or a new covariate matrix.

When the model was fitted using a formula (Formula Mode), a newdata data frame can be supplied instead of newA; the covariate matrix is then constructed automatically from the stored formula metadata.

Usage

## S3 method for class 'nmfkc'
predict(object, newA = NULL, newdata = NULL, type = "response", ...)

Arguments

Value

Depending on type: a numeric matrix ("response" or "prob") or a character vector of class labels ("class").

See Also

nmfkc, nmfkc.cv

Examples

# Prediction with newA
Y <- matrix(cars$dist, nrow = 1)
A <- rbind(1, cars$speed)
result <- nmfkc(Y, A, rank = 1)
newA <- rbind(1, c(10, 20, 30))
predict(result, newA = newA)

Print method for NMF inference objects

Description

Prints a summary of any NMF inference result object ("nmfkc.inference" or "nmfae.inference").

Usage

## S3 method for class 'nmf.inference'
print(x, ...)

Arguments

Value

Called for its side effect (printing). Returns x invisibly.

See Also

nmfkc.inference, nmfae.inference

Print method for summary.nmfae objects

Description

Prints a formatted summary of an NMF-AE model fit.

Usage

## S3 method for class 'summary.nmfae'
print(x, digits = max(3L, getOption("digits") - 3L), max.coef = 20, ...)

Arguments

Value

Called for its side effect (printing). Returns x invisibly.

See Also

summary.nmfae

Print method for summary.nmfae.inference objects

Description

Prints a formatted summary including the coefficients table.

Usage

## S3 method for class 'summary.nmfae.inference'
print(x, digits = max(3L, getOption("digits") - 3L), max.coef = 20, ...)

Arguments

Value

Called for its side effect (printing). Returns x invisibly.

See Also

summary.nmfae.inference

Print method for summary.nmfkc objects

Description

Prints a formatted summary of an nmfkc model fit.

Usage

## S3 method for class 'summary.nmfkc'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

Value

Called for its side effect (printing). Returns x invisibly.

Examples

Y <- matrix(cars$dist, nrow = 1)
A <- rbind(1, cars$speed)
result <- nmfkc(Y, A, rank = 1)
print(summary(result))

Print method for summary.nmfkc.inference objects

Description

Prints a formatted summary including the coefficients table.

Usage

## S3 method for class 'summary.nmfkc.inference'
print(x, digits = max(3L, getOption("digits") - 3L), max.coef = 20, ...)

Arguments

Value

Called for its side effect (printing). Returns x invisibly.

See Also

summary.nmfkc.inference

Extract residuals from NMF models

Description

Returns the residual matrix Y - \hat{Y} from a fitted NMF model. Requires the original observation matrix Y to be supplied.

For nmfre objects, residuals are computed from the BLUP reconstruction (Y - X(B_{blup})) by default. Set type = "fixed" to use fixed-effects only.

Usage

## S3 method for class 'nmf'
residuals(object, Y, ...)

## S3 method for class 'nmfae'
residuals(object, Y, ...)

## S3 method for class 'nmfre'
residuals(object, Y, type = c("blup", "fixed"), ...)

## S3 method for class 'nmf.sem'
residuals(object, Y, ...)

Arguments

Value

The residual matrix.

See Also

nmfkc, nmfae, nmfre, nmf.sem, fitted.nmf

Examples

Y <- matrix(runif(50), 5, 10)
result <- nmfkc(Y, rank = 2)
residuals(result, Y)

Summary method for nmf.sem objects

Description

Produces a formatted summary of a fitted NMF-SEM model, including matrix dimensions, convergence, stability diagnostics, fit statistics, and inference results (if available).

Usage

## S3 method for class 'nmf.sem'
summary(object, ...)

Arguments

Value

Invisible object.

See Also

nmf.sem, nmf.sem.inference

Examples

Y <- t(iris[, -5])
Y1 <- Y[1:2, ]; Y2 <- Y[3:4, ]
result <- nmf.sem(Y1, Y2, rank = 2, maxit = 500)
summary(result)

Summary method for nmfae objects

Description

summary.nmfae produces a summary of a fitted NMF-AE model, including dimensions, convergence status, goodness-of-fit statistics, and structure diagnostics (sparsity of factor matrices).

Usage

## S3 method for class 'nmfae'
summary(object, ...)

Arguments

Value

An object of class "summary.nmfae", a list with components:

See Also

nmfae, print.summary.nmfae

Summary method for nmfae.inference objects

Description

Produces a summary of a fitted NMF-AE model with inference results, including the coefficients table for \Theta.

Usage

## S3 method for class 'nmfae.inference'
summary(object, ...)

Arguments

Value

An object of class "summary.nmfae.inference".

See Also

nmfae.inference, summary.nmfae

Summary method for objects of class nmfkc

Description

Produces a summary of an nmfkc object, including matrix dimensions, runtime, fit statistics, and diagnostics.

Usage

## S3 method for class 'nmfkc'
summary(object, ...)

Arguments

Value

An object of class summary.nmfkc, containing summary statistics.

See Also

nmfkc, nmfkc.inference, plot.nmfkc

Examples

Y <- matrix(cars$dist, nrow = 1)
A <- rbind(1, cars$speed)
result <- nmfkc(Y, A, rank = 1)
summary(result)

Summary method for nmfkc.inference objects

Description

Produces a summary of a fitted NMF model with inference results, including the coefficients table for C (\Theta).

Usage

## S3 method for class 'nmfkc.inference'
summary(object, ...)

Arguments

Value

An object of class "summary.nmfkc.inference".

See Also

nmfkc.inference, summary.nmfkc

Summary method for objects of class nmfre

Description

Displays a concise summary of an NMF-RE model fit, including dimensions, convergence, variance components, and a coefficient table following standard R regression output conventions.

Usage

## S3 method for class 'nmfre'
summary(object, show_ci = FALSE, ...)

Arguments

Value

The input object, invisibly.

See Also

nmfre, nmfre.inference

Examples

Y <- matrix(cars$dist, nrow = 1)
A <- rbind(intercept = 1, speed = cars$speed)
res <- nmfre(Y, A, rank = 1, maxit = 5000)
summary(res)