R: Cluster based on similar linear relationships

Question

I'm looking for an unsupervised clustering technique available in R that will allow me to combine repeated measures I have taken at many independent sites, to form subgroups that have similar linear relationships (in terms of slope not intercept) with a particular covariate.

I've looked far and wide for something like this, but the closest approach I've seen is to generate correlation coefficients for each independent site, then use something like k-means clustering to group the sites based on the r value. However, this isn't exactly what I'm looking for because I'd like to derive linear regression estimates for each of the groups of sites and validate them using cross validation (or something similar).

The workaround that I came up with is to use either lmtree or lmertree from the glmertree package in R (which I was already using for another purpose) to group sites based on their individual correlation coefficients, and simultaneously estimate linear models for each group, followed by out of sample testing (not included in code below). However, I'm wondering if there's something that wouldn't necessitate partitioning based on correlation coefficients, which may or may not be accurate/stat. significant due to low sample sizes and outliers at individual sites. Examples of how to validate the clusters out-of-sample would also be appreciated.

Here's some fake data and the code that shows my current approach that I'd like to improve on:

structure(list(log_abund_centered = c(-0.48964375524682, -0.354876147884037, 
-0.596639594245377, 1.04568618503356, -0.288648437155621, 0.593478814373263, 
-0.205300584060553, 0.50382845201984, -0.969216835508706, 0.761331902674453, 
-1.03336934363746, -1.29573360810495, 0.259931157461819, 0.774531341769983, 
0.863269756411074, 0.491013319299562, 0.0561519733294005, 0.479218742806726, 
-0.813307458860655, 0.218294119524494, 3.59555207741902, 0.388981289464139, 
-0.146778038109325, -0.967067505148195, 0.537585431212312, -1.19095152299013, 
0.678585647280611, 0.386307491079737, -2.24541402704531, -1.03680084316287, 
-1.5308603150916, 0.553507344924291, -0.133385304782127, 0.50750694478563, 
-0.0774266511340822, 1.20878425142883, 0.0969267360106958, 0.787570786352522, 
-0.407856573299332, -1.00476721919482, -0.194723213221809, -0.891774107262228, 
1.75522765341157, 0.941325666806081, -0.429295460359642, -0.262241375696476, 
-0.951792861673726, 0.086423716704644), Site.Code = structure(c(5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 6L, 6L), .Label = c("Site3", 
"Site4", "Site6", "Site7", "Site8", "T228"), class = "factor"), 
    pdsi_prev_sum = c(-2.73126133282979, 1.75857214132945, 1.90214856465658, 
    5.77923981348673, 5.08813714981079, 4.28079509735108, -2.56452918052673, 
    1.85949591795603, 3.4981784025828, 3.69054238001506, -2.30718366305033, 
    1.96718434492747, 2.31010007858276, 5.60476223627726, 4.96186192830404, 
    4.18374554316203, -2.5044375260671, 2.09257872899373, 4.00702516237895, 
    2.92732580502828, -2.16150712966919, 1.75016431013743, 2.05120766162872, 
    5.25892893473307, 4.31651020050049, 3.60424534479777, -2.47522457440694, 
    2.04869890213013, 4.7067833741506, 3.36243359247843, -2.78738840421041, 
    1.57434491316478, 1.45780511697134, 5.51815017064412, 3.75366536776225, 
    3.56783596674601, -2.98427677154541, 2.7864560286204, 2.7763071457545, 
    3.15730079015096, 5.25488535563151, 4.53062645594279, -2.46269559860229, 
    2.37915563583374, 4.26401273409526, 3.55170520146688, 1.97975746790568, 
    5.87917391459147), prev_summer_cors = c(0.428816345938213, 
    0.428816345938213, 0.428816345938213, 0.428816345938213, 
    0.428816345938213, 0.428816345938213, 0.428816345938213, 
    0.428816345938213, 0.428816345938213, 0.428816345938213, 
    0.4862086001943, 0.4862086001943, 0.4862086001943, 0.4862086001943, 
    0.4862086001943, 0.4862086001943, 0.4862086001943, 0.4862086001943, 
    0.4862086001943, 0.4862086001943, -0.771777574166274, -0.771777574166274, 
    -0.771777574166274, -0.771777574166274, -0.771777574166274, 
    -0.771777574166274, -0.771777574166274, -0.771777574166274, 
    -0.771777574166274, -0.771777574166274, 0.454383299746965, 
    0.454383299746965, 0.454383299746965, 0.454383299746965, 
    0.454383299746965, 0.454383299746965, 0.454383299746965, 
    0.454383299746965, 0.454383299746965, 0.454383299746965, 
    -0.895821520966809, -0.895821520966809, -0.895821520966809, 
    -0.895821520966809, -0.895821520966809, -0.895821520966809, 
    -0.328111860245374, -0.328111860245374)), .Names = c("log_abund_centered", 
"Site.Code", "pdsi_prev_sum", "prev_summer_cors"), row.names = c(232L, 
233L, 234L, 235L, 236L, 237L, 238L, 239L, 240L, 241L, 385L, 386L, 
387L, 388L, 389L, 390L, 391L, 392L, 393L, 394L, 423L, 424L, 425L, 
426L, 427L, 428L, 429L, 430L, 431L, 432L, 468L, 469L, 470L, 471L, 
472L, 473L, 474L, 475L, 476L, 477L, 634L, 635L, 636L, 637L, 638L, 
639L, 640L, 641L), class = "data.frame")

Code:

library(glmertree)

lmer.tree<-lmertree(log_abund_centered ~ pdsi_prev_sum|(1|Site.Code)|prev_summer_cors, data = toy_data ,bonferroni=F)

Answer 1

The model-based recursive partitioning approach underlying lmtree, lmertree, and friends, can be considered a form of "supervised clustering" because you need partitioning variables to form the clusters through recursive splits.

The corresponding approach for "unsupervised clustering" is often called model-based clustering. Under this label this often pertains to multivariate models with certain types of correlation matrices in each cluster. However, it is also possible to use regression models within each cluster. In this case the models are often called cluster-wise regression. And certain latent-class analysis models also allow using regressors in each cluster/class.

An overview of various R package on CRAN is given in the "Cluster" task view https://CRAN.R-project.org/view=Cluster. There are various models/packages that allow the inclusion of both regressors and random effects (to capture correlations in repeated measurements). I'm not sure which of these fits your needs best but it should give a useful starting point.

Answer 2

After thinking on this a bit, I realized this is essentially a latent variable analysis problem, where the latent variables are the estimated slopes (essentially from a regression process) and the clustering process. I would recommend using a program such as Mplus that can model the regression parameters (either slope or intercept or other higher order powers) combined with the ability to cluster on those latent variables (along with other variables you may have measured). My understanding is that there are currently no R packages for latent class analysis (though I'd be very happy to learn that I am wrong).