Skip to contents

Fit NMF-GMM: a Gaussian-mixture latent-class extension of NMF with covariates

Source: R/nmf.gmm.R
nmf.gmm.Rd

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 fits the model $$\bm b_n\mid(z_n=k)\sim N_Q(C\bm a_n+\bm\mu_k,\Sigma_k),\qquad \bm y_n = X\bm b_n+\bm\varepsilon_n,\quad X\ge 0,$$ by a generalized EM algorithm: a \(K\)-component Gaussian mixture is placed on the latent NMF scores \(\bm b_n\), whose class-conditional mean combines a covariate regression \(C\bm a_n\) (\(C=\Theta\) in the paper, the \(Q\times R\) coefficient matrix) with a class shift \(\bm\mu_k\). The shared basis \(X\) stays non-negative and column-normalized. Clustering is model based, through the posterior responsibilities. \(K=1\) (tied covariance) reduces exactly to NMF with random effects (nmfre).

Like the rest of the package, nmf.gmm performs optimization only; use nmf.gmm.inference for standard errors and tests of \(C\), and nmf.gmm.select to choose \(K\).

Usage

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

Arguments

Y

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

A

Covariate matrix \(A\) (R x N). Default NULL uses an intercept-only \(1\times N\) matrix (a plain Gaussian mixture on the scores, with no covariate adjustment). A row of ones should be present for the identifiability constraint \(\sum_k \xi_k\bm\mu_k = 0\) (see intercept).

rank

Integer rank \(Q\) of the basis.

K

Integer number of mixture components. Default 1 (= NMF-RE).

...

Additional arguments:

  • cov: score-covariance structure. "tied" (default, shared diagonal \(\Sigma_k=\mathrm{diag}(\tau^2_1,\dots,\tau^2_Q)\), \(Q\) variances); "free" (per-class diagonal, \(Q\times K\)); or "scalar" (isotropic \(\Sigma_k=\tau^2 I_Q\), a single variance — the most parsimonious variant; at \(K=1\) it is the nmfre model, and tau2 is returned as one number).

  • X.init: initial basis. A \(P\times Q\) non-negative matrix, or an initialization method name passed to the shared initializer ("nndsvd" default, "kmeans++", "kmeans", ...).

  • intercept: index of the intercept row of A (default 1); the class-mean average is absorbed into that column of \(C\).

  • nstart: EM restarts (default 1 for K=1, else 8).

  • cores: run the nstart restarts in parallel (default getOption("mc.cores", 1L); PSOCK cluster on Windows, forking elsewhere). Each restart is an independent self-seeded EM and results are combined in order, so the returned fit is identical for any cores. (nmf.gmm.select parallelizes over K instead and runs each inner fit with cores = 1.)

  • maxit, tol: outer EM cap / tolerance (2e5, 1e-7). tol is the relative change of the marginal log-likelihood, the usual EM rule.

  • ptol: a second stopping rule, on the relative change of the reported parameters \(X\), \(C\) and \(\bm\mu\) (default 1e-9; NULL disables it). The log-likelihood rule alone can fail to fire when a nuisance variance drifts to its boundary – \(\tau^2\to 0\) on a sample that carries no random effect – so the objective keeps creeping up after the partition, the basis and the covariate effect have settled. The returned stop_by says which rule fired ("loglik", "parameters" or "maxit").

  • vfloor: lower bound on the variance parameters. The default (NULL) makes it scale-relative, \(10^{-12}\,\mathrm{mean}(Y^2)\), so the fit is equivariant under \(Y\to cY\); a fixed bound would bind on small-scale data.

  • fixX: hold the basis at its initial value instead of updating it (default FALSE). Used to compare the joint fit with nmf.gmm.twostage on the same \(X\), so that only the order of adjustment and clustering differs.

  • seed: RNG seed (default 1). prefix: basis-name prefix (default "Basis").

  • data: a data frame with one row per column of Y, required when A is a formula (see below).

  • standardize: when A is a formula, center and scale the constructed covariate columns (default TRUE). Factors are expanded to treatment indicators before standardization, so a \(5\)-level factor becomes four centered, scaled columns.

A may also be a one-sided formula (e.g. ~ size or ~ diet), evaluated in data: the design matrix is built with model.matrix, its intercept column is replaced by the package's own intercept row, and the remaining columns are standardized by default. The constructed numeric \(A\) is returned in the fit (fields A, A.formula, A.center, A.scale), so nmf.gmm.inference works unchanged.

Value

An object of class "nmf.gmm": a list with X (basis), C (\(=\Theta\), Q x R), mu (class means, Q x K), tau2, sigma2, xi (mixing proportions), gamma (responsibilities, N x K), cluster (hard labels), loglik, BIC, ICL, n.params, entropy, Yhat, objfunc(.iter), iter, converged, stop_by, K, rank, cov, dims and runtime.

See also

nmf.gmm.inference, nmf.gmm.select, nmfre (the \(K=1\) special case).

Examples

# \donttest{
set.seed(1)
Y <- matrix(abs(rnorm(20 * 60)) + 1, 20, 60)
A <- rbind(1, rnorm(60))                 # intercept + one covariate
fit <- nmf.gmm(Y, A, rank = 2, K = 2)
fit$BIC
#> [1] 2398.661
table(fit$cluster)
#> 
#>  1  2 
#> 14 46 
# }