Co-clustering two variable blocks with NMF-RRR (nmf.rrr)

Introduction

Many studies measure two blocks of variables on the same individuals – a block of covariates (inputs) and a block of responses (outputs) – and ask how groups of covariates relate to groups of responses. nmf.rrr() answers this by starting from the multivariate linear regression of the responses on the covariates and giving its non-negative regression coefficient matrix a tri-factorization

where (each column summing to one) softly clusters the response variables, (each row summing to one) softly clusters the covariate variables, and is a tested matrix of block correspondences. Because has rank , this is the non-negative, parts-based member of the reduced-rank regression (RRR) family – hence NMF-RRR – related to RRR as NMF is to PCA.

nmf.rrr() and its nmf.rrr.* helpers are the canonical interface; the names emphasise the reduced-rank-regression reading. The former nmfae() family names are retained as deprecated aliases for backward compatibility.

This vignette reproduces the Doubs example (community ecology), where a single dominant upstream–downstream gradient aligns both blocks, so the correspondence is a clean permutation.

library(nmfkc)

The Doubs data

The Doubs data (Verneaux, 1973) record 27 fish species and 11 environmental variables at 30 sites along a French river – a classic illustration of canonical (correspondence) analysis. We take the fish abundances as the response block and the environmental variables as the covariate block .

Each variable is mapped to by a per-variable min–max transform (nmfkc.normalize()), which makes the sign-free environmental variables non-negative, and both blocks are laid out as variables sites ().

data(doubs, package = "ade4")

# per-variable min-max to [0,1], then transpose to (variables x sites)
nz <- function(M) t(nmfkc.normalize(as.matrix(M)))
Y1 <- nz(doubs$fish)   # responses: 27 fish species x 30 sites
Y2 <- nz(doubs$env)    # covariates: 11 environment x 30 sites
dim(Y1)
#> [1] 27 30
dim(Y2)
#> [1] 11 30

Fitting NMF-RRR

Element-wise cross-validation (below) selects . Signed models and inference benefit from several k-means restarts, so we set nstart = 20 and a tight tolerance.

fit <- nmf.rrr(Y1, Y2, rank1 = 2, rank2 = 2,
               epsilon = 1e-8, nstart = 20, seed = 1)

# in-sample, column-centered R^2
Y1hat <- fit$X1 %*% fit$C %*% fit$X2 %*% Y2
R2 <- 1 - sum((Y1 - Y1hat)^2) / sum((Y1 - rowMeans(Y1))^2)
round(R2, 3)
#> [1] 0.435

The fit is R2 = 0.435, reproducing the classical longitudinal zonation of the river.

for (q in 1:ncol(fit$X1))
  cat(sprintf("Resp%d: %s\n", q,
      paste(rownames(Y1)[order(-fit$X1[, q])[1:6]], collapse = ", ")))
#> Resp1: Neba, Phph, Satr, Cogo, Thth, Teso
#> Resp2: Alal, Ruru, Gogo, Acce, Baba, Titi

for (r in 1:nrow(fit$X2))
  cat(sprintf("Cov%d: %s\n", r,
      paste(rownames(Y2)[order(-fit$X2[r, ])[1:5]], collapse = ", ")))
#> Cov1: dfs, flo, har, nit, pH
#> Cov2: oxy, alt, pH, slo, har

Choosing the two ranks

Because the attainable fit is bounded by , the in-sample fit cannot choose the ranks; we use element-wise cross-validation (nmf.rrr.ecv()), which holds out entries of and predicts them.

ecv <- nmf.rrr.ecv(Y1, Y2, rank1 = 1:2, rank2 = 1:2,
                   nfolds = 5, seed = 123)
#> Element-wise CV: 4 (Q,R) pairs, 5-fold, 20 tasks...
#>   Q=1, R=1: MSE=0.090948, sigma=0.3016
#>   Q=2, R=1: MSE=0.090948, sigma=0.3016
#>   Q=1, R=2: MSE=0.090954, sigma=0.3016
#>   Q=2, R=2: MSE=0.065768, sigma=0.2565
round(ecv$sigma, 4)
#> Q=1,R=1 Q=2,R=1 Q=1,R=2 Q=2,R=2 
#>  0.3016  0.3016  0.3016  0.2565

The smallest hold-out error is at .

Inference for the correspondence matrix

The entries of say how strongly each covariate group drives each response group. nmf.rrr.inference() attaches standard errors (Fisher + wild bootstrap) and a one-sided boundary test (each ).

inf <- nmf.rrr.inference(fit, Y1, Y2)
co  <- inf$coefficients
print(format(co[order(co$p_value), c("Basis","Covariate","Estimate","SE","z_value","p_value")],
             digits = 3))
#>   Basis Covariate Estimate    SE  z_value  p_value
#> 3 Resp1      Cov2 3.97e+00 0.505 7.86e+00 1.93e-15
#> 2 Resp2      Cov1 1.40e+01 1.837 7.65e+00 1.03e-14
#> 1 Resp1      Cov1 2.06e-48 0.603 3.42e-48 5.00e-01
#> 4 Resp2      Cov2 9.54e-27 1.004 9.51e-27 5.00e-01

is a near-permutation: the upstream guild is driven by the oxic gradient and the downstream guild by the nutrient gradient (both ), while the two off-diagonal paths are essentially zero ().

round(fit$C, 3)
#>         Cov1  Cov2
#> Resp1  0.000 3.971
#> Resp2 14.049 0.000

Visualising and the two co-clusterings

nmf.rrr.heatmap() shows the response basis , the correspondence , and the covariate basis together.

nmf.rrr.heatmap(fit)

Relation to other methods

Dropping non-negativity, at rank is ordinary reduced-rank regression (RRR): it attains a higher in-sample fit ( on Doubs) but returns signed loadings and no clusters. On these data the two share the dominant fitted direction (leading principal-angle cosine ); they differ in the basis of that subspace – non-negative parts versus signed singular directions – exactly as NMF relates to PCA. An unsupervised tri-NMF of the association recovers the same guilds and gradients here (the gradient is so dominant that supervised and unsupervised co-clusterings coincide), but, unlike NMF-RRR, cannot predict the community at a new site through .

When within-block and cross-block structure disagree – e.g. under – the non-negative, normalized parameterization of NMF-RRR stays well-behaved and exposes cross-structure (one response group driven by several covariate groups) that these baselines miss; see the paper for the nutrimouse and microbiome–metabolome examples.

References

• Satoh, K. & Tokuda, Y. Co-clustering of Response and Covariate Variables by Tri-Factorizing Their Non-negative Regression Coefficient Matrix (manuscript).
• Ding, C., Li, T., Peng, W. & Park, H. (2006). Orthogonal nonnegative matrix tri-factorizations for clustering. KDD.
• Verneaux, J. (1973). Cours d’eau de Franche-Comté. PhD thesis. (Doubs data, R package ade4.) ```