Two-stage (adjust-then-cluster) baseline for NMF-GMM

This function is experimental and still under development. The interface may change in future versions: argument names, defaults and the contents of the returned object are not yet stable. Code written against it today may need adjusting after an update. The rest of the package does not carry this caveat.

nmf.gmm.twostage runs the two-stage route that nmf.gmm is designed to improve on, as a matched baseline: (1) estimate least-squares scores on the initial basis, (2) regress the scores on the covariates blind to the class and keep the residuals, (3) reconstitute the residuals in observation space, shift them to non-negativity, and (4) refit an intercept-only nmf.gmm from the same basis initialization. Only the order of adjustment and clustering differs from the joint fit, so the pair isolates the displacement of the class means that two-stage adjustment incurs when the covariate is associated with the class (Satoh 2026, Proposition 4); when the covariate is (near-)mean-independent of the class the two routes agree.

Usage

nmf.gmm.twostage(Y, A = NULL, rank, K = 1, ...)

Arguments

Y

Data matrix \(Y\) (P x N).

A

Covariate matrix \(A\) (R x N) including an intercept row, or a one-sided formula evaluated in data (as in nmf.gmm). An intercept-only A is an error: there is nothing to adjust for.

rank

Integer rank \(Q\) of the basis.

K

Integer number of mixture components.

...

Additional arguments as in nmf.gmm (cov, X.init, nstart, maxit, seed, data, standardize, ...); they are applied to both stages.

Value

An object of class c("nmf.gmm.twostage", "nmf.gmm"): the stage-2 fit (all nmf.gmm fields and S3 methods apply), plus a twostage list with the non-negativity shift, the covariate matrix A that was removed, and the shared basis initialization X0.

See also

nmf.gmm (the joint route this baselines).

Examples

# \donttest{
set.seed(1)
Y <- matrix(abs(rnorm(12 * 40)) + 1, 12, 40)
A <- rbind(1, rnorm(40))
ts <- nmf.gmm.twostage(Y, A, rank = 2, K = 2, nstart = 2, maxit = 100)
table(ts$cluster)
#> 
#>  1  2 
#>  6 34 
# }