Dissimilarities between scalar GLM response channels

Compute pairwise dissimilarities between columns of a response matrix. method="poissondist" uses Witten's size-adjusted, smoothed Poisson likelihood gain with a native OpenMP kernel. method="deviance" pools each pair at its precision-weighted mean and sums their unit deviances.

Usage

glm_dist(
  Y,
  weights = NULL,
  family = poisson(identity),
  method = c("deviance", "poissondist"),
  size_fac = colSums(Y),
  nthreads = NULL
)

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.

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.

Value

A labelled dist object in packed lower-triangular storage.

Details

For channels \(a,b\), the pooled-deviance dissimilarity is $$D_{ab}=\sum_i\{d(y_{ia},\bar\mu_i,w_{ia})+ d(y_{ib},\bar\mu_i,w_{ib})\},\qquad \bar\mu_i=\frac{w_{ia}y_{ia}+w_{ib}y_{ib}}{w_{ia}+w_{ib}}.$$ Native pooled-deviance kernels cover Gaussian, Poisson, binomial, Gamma, inverse Gaussian and fixed-shape NB families; custom scalar families use their R dev.resids function. NB shape is taken from the family object and held fixed. Rows with zero total weight contribute zero. For a regular scalar exponential family with convex mean potential \(F\), \(d(y,\mu,w)=2w B_F(y,\mu)\). The pooled loss is twice a weighted Jensen gap. This connects the dissimilarity to exponential-family likelihood and Bregman clustering (Banerjee et al., 2005). For Gaussian observations with unit weights it is \(\tfrac12\sum_i(y_{ia}-y_{ib})^2\). A generic pooled unit deviance is a dissimilarity; distributional Jensen–Shannon terminology and a universal metric claim do not apply. nb1(), Gamma_scale() and Gamma_censored() use the exact likelihood ratio between separate fitted means and a shared fitted mean. Their native scalar optimizer accounts for censoring and the specified fixed noise parameter. The fitted means need not equal the raw observations.

For the Poisson method let \(s_a,s_b\) be the size factors and define $$c_i=\frac{w_{ia}y_{ia}+w_{ib}y_{ib}} {w_{ia}s_a+w_{ib}s_b},\quad m_{ia}=s_a c_i,\quad r_{ia}=\frac{w_{ia}y_{ia}+1}{w_{ia}m_{ia}+1}.$$ The contribution for channel \(a\) is \(w_{ia}\{m_{ia}-m_{ia}r_{ia}+y_{ia}\log r_{ia}\}\), with the analogous contribution for \(b\). Unit weights recover Witten's Gamma(1,1)-smoothed likelihood-gain formula. The weighted version applies the same smoothing to the weighted count and expected count. Zero-weight and zero-exposure terms contribute zero. Common rescaling of all size factors leaves the model unchanged.

References

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

Banerjee, A., Merugu, S., Dhillon, I. S. and Ghosh, J. (2005). Clustering with Bregman Divergences. JMLR 6, 1705–1749. https://jmlr.org/papers/v6/banerjee05b.html.

See also

glm_cluster()

Examples

## example: calculating distance among different m/z mass traces in GC/MS datasets
data(GCMS_data2) # section of a GC/MS chromatogram of beer volatiles
Y <- GCMS_data2
dim(Y) # 1764 timepoints x 518 m/z mass traces
# we aggregate together every bin of 10 subsequent timepoints for extra performance:
Y2 <- t(glm_aggregate(t(Y), 
                    family = poisson(identity), 
                    clusterids = rep(1L:ceiling(nrow(Y)/10), each=10)[1:nrow(Y)])$Y_agg)
dim(Y2) # 177 time bins x 518 m/z mass channels
D1 <- glm_dist(Y = Y2, family  = poisson(identity), method="poissondist")
# identical to D1b <- PoiClaClu::PoissonDistance(t(Y2), type = "mle", transform = FALSE)$dd
# but faster