Likelihood-based clustering and aggregation of GLM response channels

Cluster columns of a scalar-family response matrix, select a tree cut by an information criterion, and aggregate the original, unbinned channels. Native OpenMP kernels calculate distances, Bregman merges and profile losses.

Usage

glm_cluster(
  Y,
  weights = NULL,
  family = poisson(identity),
  method = if (family$family == "poisson") "poissondist" else "deviance",
  size_fac = colSums(Y),
  rowbins = 1L,
  nthreads = NULL,
  linkage = NULL,
  criterion = NULL,
  k = NULL,
  k.max = 32L,
  dispersion = 1
)

glm_aggregate(Y, weights = NULL, family = poisson(identity), clusterids)

Arguments

Y: Finite numeric observations-by-channels matrix. Sparse matrices are accepted; response calculations use a dense matrix internally.
weights: Nonnegative case/precision weights: a scalar, a vector with one entry per observation, or an observations-by-channels matrix. NULL uses unit weights. Binomial observations are proportions and weights give the trial totals.
family: Scalar GLM family object, default poisson(identity). Response values must lie in the family's mean/deviance domain.
method: "deviance" (default) or "poissondist".
size_fac: Nonnegative channel size factors for "poissondist", default colSums(Y). The same factors apply to both Ward and Bregman clustering when method="poissondist"; method="deviance" uses raw channel means. Supplying rep(1,ncol(Y)) disables size adjustment. A zero factor is allowed for an uninformative zero-count channel. Counts remain on their original scale.
rowbins: Positive integer, default 1L. Values above one sum adjacent rows for clustering only; use this acceleration for ordered Poisson scans with unit weights. Returned aggregates retain every original row.
nthreads: Number of native worker threads. NULL uses the detected physical cores; 1L requests serial execution. Results have the same within-pair summation order at every thread count.
linkage: "ward" or "bregman". The size-adjusted Poisson distance defaults to the Ward hierarchy. Regular scalar GLM likelihoods default to Bregman merges; signal_families use Ward likelihood-ratio distances.
criterion: "gic" (default), "aic", "bic", "ebic", "mbic", or "dynamic" for dynamicTreeCut on a Ward hierarchy. Quasibinomial and custom unit-deviance families default to dynamic cutting. For quasipoisson, supply a common dispersion to obtain quasi-information criteria. Supplying k requests that fixed number of groups.
k: Optional fixed number of groups, between one and ncol(Y).
k.max: Largest number of groups considered by criterion selection, default 32L, capped at ncol(Y).
dispersion: Known positive scalar or per-channel dispersion for Gaussian, Gamma and inverse Gaussian clustering. Default 1 gives unit dispersion. It scales clustering precision and is held fixed across cuts. Poisson and binomial use dispersion one; NB uses its fixed family shape.
clusterids: One positive integer label per channel. Labels need not be consecutive; glm_aggregate() puts groups in increasing label order.

Value

glm_cluster() returns a list with D (packed distances for Ward), hc (hierarchy; NULL for one channel), clusterids, Y_agg, weights_agg, history (cut losses, degrees of freedom and criterion), centres on the clustering row grid, size_fac, linkage, criterion and rowbins. Incompatible count weights give Y_agg=NULL and an explanatory aggregation_note; clustering itself remains available. glm_aggregate() returns Y_agg and weights_agg.

Details

For a partition $\mathcal G$, the raw-mean clustering model assigns a free mean $\mu_{ig}$ at each row $i$ of group $g$. Its profile distortion is $$L(\mathcal G)=\sum_g\sum_{j\in G_g}\sum_i \frac{w_{ij}}{\phi_j}d(y_{ij},\widehat\mu_{ig},1),\qquad \widehat\mu_{ig}=\frac{\sum_{j\in G_g}w_{ij}y_{ij}/\phi_j} {\sum_{j\in G_g}w_{ij}/\phi_j}.$$ Bregman clustering greedily merges the pair minimising $L(A\cup B)-L(A)-L(B)$. The implementation retains sufficient statistics and a priority queue, recalculating the true likelihood-merge cost after each merge. This is agglomerative Bregman information loss (Banerjee et al., 2005; Telgarsky and Dasgupta, 2012).

For method="poissondist", clustering profiles global channel intensity using fixed size factors: $\mu_{ij}=s_j c_{ig}$ and $\widehat c_{ig}=\sum_j w_{ij}y_{ij}/\sum_j w_{ij}s_j$. The Ward route uses Witten's smoothed pairwise distance followed by fastcluster::hclust(method="ward.D"). The Bregman route uses the unsmoothed Poisson likelihood-merge cost for this size-factor model. Ward's recurrence is exact for squared-Euclidean sum-of-squares geometry; general deviance dissimilarities give a computational Ward hierarchy. Both routes score their candidate cuts under the specified profile model. Thus the same size factors adjust both linkages. The factors are nuisance exposures used to compare profile shapes; the regression stages continue to use raw counts, and the returned aggregates sum those raw counts.

For each candidate cut the criterion is $$L(\mathcal G)-2\sum_g m_g\log(m_g/m)+c(n,n)q_{\mathcal G}.$$ Here $m_g$ is the number of channels in group $g$, the second term is the fitted classification-mixture allocation loss, and $q_{\mathcal G}$ counts informative row/group means plus $K-1$ mixing proportions. $c$ is the package IC penalty and $n$ is the original number of rows (at least three for scoring). The saturated likelihood constant is common to all cuts and omitted. For GIC, $c(n,n)=\log(n)\log\log(n)$. This selects a cut within the fitted hierarchy; the regression likelihood and predictor count used by L0glm.fit() are separate. Dispersion and NB shape are fixed here rather than re-estimated at every cut. Gaussian/Gamma/inverse Gaussian dispersion should be supplied when its known scale differs from one.

Native scalar kernels support Gaussian, Poisson, binomial proportions, Gamma, inverse Gaussian and fixed-shape negative binomial families. Joint multinomial and Cox responses have different likelihood structures. Quasipoisson and quasibinomial have native working-deviance kernels and support fixed cuts or Ward/dynamic cutting. Custom scalar unit deviances remain available through glm_dist() and the Ward/dynamic route. NB1, common-scale Gamma and censored Gamma use the fixed-parameter exact likelihoods in signal_families. Their pooled means are optimized with a native scalar score solve, and the saturated reference uses separate optimized means. They support Ward clustering with IC, fixed or dynamic cuts.

Aggregation summarises the original columns. Poisson counts and their means add. Shared Poisson case weights stay unchanged; differing weights within a group cannot be combined into one exact weighted Poisson model. Gaussian means use precision-weighted averaging and summed precisions, including the supplied dispersion in those precisions. Binomial proportions use trial-weighted proportions. For Gamma, inverse Gaussian and NB, the returned summaries do not generally have the same fitted GLM family/link/dispersion. Clustering support therefore does not imply that automatic response aggregation is valid for every GLM.

Functions

glm_aggregate(): Summarise the original channels using supplied groups.

References

Banerjee, A. et al. (2005). Clustering with Bregman Divergences. JMLR 6, 1705–1749. https://jmlr.org/papers/v6/banerjee05b.html.

Telgarsky, M. and Dasgupta, S. (2012). Agglomerative Bregman Clustering. ICML. https://icml.cc/2012/papers/758.pdf.

Witten, D. M. (2011). Classification and Clustering of Sequencing Data using a Poisson Model. AOAS 5(4), 2493–2518. doi:10.1214/11-AOAS493 .

Examples

# Native likelihood clustering of response channels
set.seed(12)
Y <- matrix(rpois(120*12, rep(c(2, 12), each=120*6)), 120, 12)
cl <- glm_cluster(Y, k.max=6L, nthreads=2L)
cl$history
cl$clusterids
dim(cl$Y_agg)

# A fixed-variance Gaussian example uses exact Bregman likelihood merges.
Z <- cbind(matrix(rnorm(60*5, -2, .5),60),
           matrix(rnorm(60*5,  2, .5),60))
normal <- glm_cluster(Z, family=gaussian(), dispersion=.25, k.max=5L)
normal$clusterids

# \donttest{
data(GCMS_data2)
Y <- as.matrix(GCMS_data2)
# Ordered scans allow coarsening for the clustering stage.
ward <- glm_cluster(Y, method="poissondist", linkage="ward",
                    size_fac=colSums(Y), rowbins=5L, k.max=32L)
bregman <- glm_cluster(Y, method="poissondist", linkage="bregman",
                       size_fac=colSums(Y), rowbins=5L, k.max=32L)
matplot(ward$Y_agg^.2,type="l",lty=1,xlab="Scan",ylab="Summed intensity ^ 0.2")
# }