Clusters of effects curves and waveforms with clustEff

1. The dissimilarity

Sottile and Adelfio (2019) proposed a dissimilarity between curves $\beta_1(p), \dots, \beta_q(p)$ evaluated on a common grid that combines

their shape: $d_{shape}(p) = I\{sign(\beta_i''(p)) = sign(\beta_{i'}''(p))\}$, the concordance of convexity;
their closeness: $d_{distance}(p) = I\{|\beta_i(p) - \beta_{i'}(p)| \le f_\alpha(p)\}$, where $f_\alpha(p)$ is the $\alpha$-quantile of the distances between all pairs of curves at $p$;

as $d(i, i') = 1 - \int d_{shape}(p)\, d_{distance}(p)\,dp$. The dissimilarity is computed by distshape() and used by clustEff() for hierarchical clustering.

2. Waveforms

We simulate 16 noisy waveforms from three shapes, as in Section 4.2 of the paper, with a noise variance that changes over time.

set.seed(1234)
tt <- seq(0, 1, length.out = 200)
sdt <- 0.3 * sqrt(pmax(4 * pmax(tt - 0.2, 0) - 8 * pmax(tt - 0.5, 0) + 4 * pmax(tt - 0.8, 0), 0))
W <- cbind(sapply(1:6, function(i) sin(3 * pi * tt) + rnorm(200, 0, sdt)),
           sapply(1:6, function(i) cos(3 * pi * tt) + rnorm(200, 0, sdt)),
           sapply(1:4, function(i) sin(3 * pi * tt) * cos(pi * tt) + rnorm(200, 0, sdt)))

The shape component uses second derivatives, so noisy curves are smoothed first (smooth = TRUE, smoothing splines). With waveforms the default $\alpha$ is 0.5.

fit <- clustEff(W, p = tt, k = c(2, 6), smooth = TRUE)
fit
#> 
#> Call:
#> clustEff(curves = W, p = tt, k = c(2, 6), smooth = TRUE)
#> 
#> Clusters of curves: 3 clusters of sizes 4, 6, 6
#> Dissimilarity: step = "both", alpha = 0.5, derivative = 2, combine = "product"
#> Agglomeration method: ward.D2   criterion: height
summary(fit)
#> 
#> Call:
#> clustEff(curves = W, p = tt, k = c(2, 6), smooth = TRUE)
#> 
#> n. of curves: 16    n. of grid points: 200 
#> n. of clusters: 3 (selected by "height" in 2-6) 
#> 
#> Clusters:
#>   1: X13, X14, X15, X16
#>   2: X1, X2, X3, X4, X5, X6
#>   3: X7, X8, X9, X10, X11, X12
#> 
#> Within-cluster measures:
#>          size    area rho.dist silhouette dissimilarity
#> cluster1    4 0.02616 0.007818     0.6548       0.22948
#> cluster2    6 0.03019 0.004706     0.7633       0.15745
#> cluster3    6 0.01637 0.001423     0.9233       0.05729
#> Average    16 0.02400 0.004253     0.7962       0.13790
#> 
#> Criterion (height) by number of clusters:
#>       2       3       4       5       6 
#> 0.66260 1.05400 0.05500 0.02628 0.11090

oldpar <- par(mfrow = c(2, 2))
plot(fit, "dendrogram")
plot(fit, "clusters")
plot(fit, "silhouette")
plot(fit, "criterion")

par(oldpar)

3. Effects curves of quantile regression coefficients models

In quantile regression coefficients modeling (QRCM, Frumento and Bottai, 2016; package ‘qrcm’) each coefficient is a parametric function of the quantile order, $\beta_j(p \mid \theta) = \theta_j^\top b(p)$. extract_curves() extracts the effects curves, their pointwise confidence bands and their second derivatives (computed from the fitted coefficient functions) from one or more qrcm::iqr() fits.

Following the multivariate example of the paper, we fit one model for each of 9 responses and cluster the effects of the covariate x, which belong to three groups.

set.seed(1)
n <- 500
x <- runif(n, 0, 5)
slope <- list(function(u) 0.5 + 0.5 * u + u^2,
              function(u) -0.5 + 0.5 * u + u^2,
              function(u) 0.3 - 0.5 * u - u^2 + 2 * u^3)
fits <- lapply(1:9, function(r) {
  u <- runif(n)
  y <- 1 + qnorm(u) + slope[[(r - 1) %/% 3 + 1]](u) * x
  qrcm::iqr(y ~ x, formula.p = ~ I(p) + I(p^2) + I(p^3))
})
names(fits) <- paste0("y", 1:9)
ec <- extract_curves(fits, which = "x")
ec
#> 
#> Effects curves: 9 curves evaluated at 99 values of p in [ 0.01 , 0.99 ]
#> Pointwise bands at level 0.95 
#> 
#>  model coefficient
#>     y1           x
#>     y2           x
#>     y3           x
#>     y4           x
#>     y5           x
#>     y6           x
#>     y7           x
#>     y8           x
#>     y9           x
fit2 <- clustEff(ec, k = 3)
summary(fit2)
#> 
#> Call:
#> clustEff(curves = ec, k = 3)
#> 
#> n. of curves: 9    n. of grid points: 99 
#> n. of clusters: 3 
#> 
#> Clusters:
#>   1: y1:x, y2:x, y3:x
#>   2: y4:x, y5:x, y6:x
#>   3: y7:x, y8:x, y9:x
#> 
#> Within-cluster measures:
#>          size    area rho.dist silhouette dissimilarity
#> cluster1    3 0.10194 0.014290     0.3223        0.6531
#> cluster2    3 0.08089 0.004999     0.2991        0.5629
#> cluster3    3 0.02676 0.025465     0.7302        0.2143
#> Average     9 0.06986 0.014918     0.4505        0.4768
oldpar <- par(mfrow = c(1, 2))
plot(fit2, "dendrogram")
plot(fit2, "clusters")

par(oldpar)

With confidence bands the default $\alpha$ is 0.25: two curves are close at $p$ if their distance is among the 25% smallest distances between pairs of curves. With 9 curves in three groups of three, the pairs within groups are exactly 25% of all pairs, so $\alpha$ should not be smaller than this proportion. Here the number of clusters is fixed after inspecting the dendrogram; with only nine noisy estimated curves the automatic criteria are less reliable. fit2$signif indicates, for each cluster, where its mean band excludes zero.

4. Choosing the number of clusters

When k is a range, criterion selects the number of clusters:

"height" (default): the largest gap between consecutive heights of the dendrogram;
"pi.out": average proportion of curves outside the mean bands of their cluster (Sottile and Adelfio, 2019, Section 3.1), elbow rule;
"dist.rel": $l\,\max_j$ of the mean $L_1$ distance from the cluster mean (Section 3.1), first number of clusters at which it stops decreasing;
"silhouette": average silhouette width (default of ‘clustEff’ 0.3.x).

A small simulation with two to four groups of polynomial effects curves gives an idea of their behaviour (proportion of replicates in which the true number of groups is selected):

p <- seq(0.01, 0.99, 0.01)
th <- list(c(.5, .5, 1, 2), c(-3, .5, 1, .5), c(.3, -.5, -1, 2), c(2.5, .5, 1, 2))
crits <- c("height", "pi.out", "dist.rel", "silhouette")
set.seed(3)
res <- t(sapply(1:12, function(r) {
  G <- 2 + r %% 3
  B <- do.call(cbind, lapply(1:G, function(g) sapply(1:8, function(i) {
    t <- th[[g]] + rnorm(4, 0, 0.3); t[1] + t[2] * p + t[3] * p^2 + t[4] * p^3 })))
  sapply(crits, function(cr) clustEff(B, B - 0.6, B + 0.6, p = p, k = c(2, 7), criterion = cr)$k == G)
}))
colMeans(res)
#>     height     pi.out   dist.rel silhouette 
#> 0.91666667 0.33333333 0.83333333 0.08333333

5. FPCAC

The FPCAC algorithm of Adelfio et al. (2011) represents the curves on a B-spline basis, computes their functional principal component scores and applies a trimmed k-means to the scores; the trimmed curves (cluster 0) are potential outliers.