MixedTypeData
Introduction
BCClong is an R package for performing Bayesian
Consensus Clustering (BCC) model for clustering continuous, discrete and
categorical longitudinal data, which are commonly seen in many clinical
studies. This document gives a tour of BCClong package.
see help(package = "BCClong") for more information and
references provided by citation("BCClong")
To download BCClong, use the following commands:
require("devtools")
devtools::install_github("ZhiwenT/BCClong", build_vignettes = TRUE)
library("BCClong")To list all functions available in this package:
Components
Currently, there are 5 function in this package which are BCC.multi, BayesT, model.selection.criteria, traceplot, trajplot.
BCC.multi function performs clustering on mixed-type (continuous, discrete and categorical) longitudinal markers using Bayesian consensus clustering method with MCMC sampling and provide a summary statistics for the computed model. This function will take in a data set and multiple parameters and output a BCC model with summary statistics.
BayesT function assess the model goodness of fit by calculate the discrepancy measure T(, ) with following steps (a) Generate T.obs based on the MCMC samples (b) Generate T.rep based on the posterior distribution of the parameters (c) Compare T.obs and T.rep, and calculate the P values.
model.selection.criteria function calculates DIC and WAIC for the fitted model traceplot function visualize the MCMC chain for model parameters trajplot function plot the longitudinal trajectory of features by local and global clustering
Pre-process (Setting up)
In this example, the PBCseq data in the
mixAK package was used as it is a public data set. The
variables used here include lbili, platelet, and spiders. Of these three
variables, lbili and platelet are continuous variables, while spiders
are categorical variables.
Fit BCC Model Using BCC.multi Function
Here, We used a binomial distribution for spiders marker, a gaussian distribution for the lbili marker and poisson distribution for platelet, respectively. The number of clusters was set to 2. All hyper parameters were set to default.
We ran the model with 12,000 iterations, discard the first 2,000 sample, and kept every 10th sample. This resulted in 1,000 samples for each model parameter. The MCMC sampling process took about 30 minutes on an AMD RyzenTM 5 5600X desktop computer.
Since this program takes a long time to run, here we will use the pre-compile result in this example.
set.seed(89)
fit.BCC2 <- BCC.multi(
mydat = list(PBC910$lbili,PBC910$platelet,PBC910$spiders),
dist = c("gaussian","poisson","binomial"),
id = list(PBC910$id),
time = list(PBC910$month),
formula =list(y ~ time + (1|id),y ~ time + (1|id), y ~ time + (1|id)),
num.cluster = 2,
burn.in = 100,
thin = 10,
per = 10,
max.iter = 200) To run the pre-compiled result, you can use
data(PBCseqfit) to attach the fitted model.
Printing Summary Statistics for key model parameters
To print the BCC model
To print the summary statistics for all parameters
To print the proportion for each cluster (mean, sd, 2.5% and 97.5% percentile) geweke statistics (geweke.stat) between -2 and 2 suggests the parameters converge
The code below prints out all major parameters
summary(fit.BCC2)
#> Total number of individual:
#> [1] 260
#>
#> Number of features:
#> [1] 3
#>
#> Cluster proportions statistics for global clusters:
#> V1 V2
#> Min. :0.8496 Min. :0.005276
#> 1st Qu.:0.8762 1st Qu.:0.043312
#> Median :0.9371 Median :0.062932
#> Mean :0.9216 Mean :0.078430
#> 3rd Qu.:0.9567 3rd Qu.:0.123799
#> Max. :0.9947 Max. :0.150384
#>
#> Globle clusters table:
#>
#> 1
#> 260
#>
#> Adherence parameters statistics by feature:
#> V1 V2 V3
#> Min. :0.5007 Min. :0.5001 Min. :0.7772
#> 1st Qu.:0.5021 1st Qu.:0.5089 1st Qu.:0.7902
#> Median :0.5030 Median :0.5115 Median :0.7959
#> Mean :0.5042 Mean :0.5175 Mean :0.8078
#> 3rd Qu.:0.5041 3rd Qu.:0.5237 3rd Qu.:0.8275
#> Max. :0.5118 Max. :0.5583 Max. :0.8576
#>
#> Local clusters statistics by feature:
#> Cluster statistics for feature 1 :
#> , , 1
#>
#> [,1] [,2]
#> mean 1.46526216 -0.14926529
#> sd 0.03777137 0.01364611
#> 2.5% 1.40888462 -0.16946273
#> 97.5% 1.51895826 -0.13096229
#> geweke.stat -2.10391018 -1.75995576
#>
#> , , 2
#>
#> [,1] [,2]
#> mean 0.014619346 0.0034315243
#> sd 0.002271399 0.0016688712
#> 2.5% 0.010935934 0.0008741408
#> 97.5% 0.017225695 0.0056475790
#> geweke.stat 0.822043653 -1.7521732959
#>
#> Cluster statistics for feature 2 :
#> , , 1
#>
#> [,1] [,2]
#> mean 5.77898717 5.264695558
#> sd 0.01067935 0.007645266
#> 2.5% 5.75976591 5.252509447
#> 97.5% 5.79323839 5.276789296
#> geweke.stat -8.63291491 -6.594456992
#>
#> , , 2
#>
#> [,1] [,2]
#> mean -0.0016958877 -0.0077753262
#> sd 0.0002249505 0.0001821189
#> 2.5% -0.0020175319 -0.0080937406
#> 97.5% -0.0013641913 -0.0075888663
#> geweke.stat 0.2202520953 0.9780291395
#>
#> Cluster statistics for feature 3 :
#> , , 1
#>
#> [,1] [,2]
#> mean -2.4027089 1.7287061
#> sd 0.2711898 0.5153692
#> 2.5% -2.9158800 1.2182696
#> 97.5% -2.0870933 2.6400038
#> geweke.stat -6.4887491 -2.3314718
#>
#> , , 2
#>
#> [,1] [,2]
#> mean 0.030690608 0.024356269
#> sd 0.017911688 0.018026932
#> 2.5% 0.006128959 -0.003560592
#> 97.5% 0.056890173 0.050779581
#> geweke.stat 1.272508080 -0.934336018
#>
#>
#> Variance-covariance matrix statistics for random effects by feature:
#> Variance-covariance matrix statistics for feature 1 :
#> , , 1
#>
#> [,1]
#> mean 3.825797e-05
#> sd 4.405348e-06
#> 2.5% 3.241441e-05
#> 97.5% 4.558949e-05
#> geweke.stat -1.102672e+00
#>
#> , , 2
#>
#> [,1]
#> mean 2.373843e-05
#> sd 3.102229e-06
#> 2.5% 1.956378e-05
#> 97.5% 2.803884e-05
#> geweke.stat -9.137537e-01
#>
#> Variance-covariance matrix statistics for feature 2 :
#> , , 1
#>
#> [,1]
#> mean 0.0016138846
#> sd 0.0005110166
#> 2.5% 0.0010801860
#> 97.5% 0.0024609516
#> geweke.stat -3.6092736173
#>
#> , , 2
#>
#> [,1]
#> mean 0.011132440
#> sd 0.002629614
#> 2.5% 0.007708163
#> 97.5% 0.014974559
#> geweke.stat -7.067780574
#>
#> Variance-covariance matrix statistics for feature 3 :
#> , , 1
#>
#> [,1]
#> mean 4.569681e-05
#> sd 6.677545e-06
#> 2.5% 3.448587e-05
#> 97.5% 5.383768e-05
#> geweke.stat 5.774842e-01
#>
#> , , 2
#>
#> [,1]
#> mean 7.637431e-05
#> sd 2.080340e-05
#> 2.5% 5.194423e-05
#> 97.5% 1.147128e-04
#> geweke.stat -1.853548e+00
#>
#>
#> Residual variance of continuous features statistics by feature:
#> Residual variance statistics for feature 1 :
#> [,1] [,2]
#> mean 0.344304859 0.344304859
#> sd 0.018691680 0.018691680
#> 2.5% 0.314145269 0.314145269
#> 97.5% 0.365692696 0.365692696
#> geweke.stat -0.003239043 -0.003239043
#> Residual variance statistics for feature 2 :
#> NULL
#> Residual variance statistics for feature 3 :
#> NULL
#>
#> Local clusters tables by feature:
#> Clusters table for feature 1 :
#>
#> 1 2
#> 86 174
#> Clusters table for feature 2 :
#>
#> 1 2
#> 126 134
#> Clusters table for feature 3 :
#>
#> 1 2
#> 198 62Visualize Clusters
Generic plot can be used on BCC object, all relevant plots will be generate one by one using return key
We can use the traceplot function to plot the MCMC process and the trajplot function to plot the trajectory for each feature.
gp1 <- trajplot(fit=fit.BCC2,feature.ind=1,which.cluster = "local.cluster",
title= bquote(paste("Local Clustering (",hat(alpha)[1] ==.(round(fit.BCC2$alpha[1],2)),")")),
xlab="months",ylab="lbili",color=c("#00BA38", "#619CFF"))
gp2 <- trajplot(fit=fit.BCC2,feature.ind=2,which.cluster = "local.cluster",
title= bquote(paste("Local Clustering (",hat(alpha)[2] ==.(round(fit.BCC2$alpha[2],2)),")")),
xlab="months",ylab="platelet",color=c("#00BA38", "#619CFF"))
gp3 <- trajplot(fit=fit.BCC2,feature.ind=3,which.cluster = "local.cluster",
title= bquote(paste("Local Clustering (",hat(alpha)[3] ==.(round(fit.BCC2$alpha[3],2)),")")),
xlab="months",ylab="spiders",color=c("#00BA38", "#619CFF"))
gp4 <- trajplot(fit=fit.BCC2,feature.ind=1,which.cluster = "global.cluster",
title="Global Clustering",
xlab="months",ylab="lbili",color=c("#00BA38", "#619CFF"))
gp5 <- trajplot(fit=fit.BCC2,feature.ind=2,which.cluster = "global.cluster",
title="Global Clustering",
xlab="months",ylab="platelet",color=c("#00BA38", "#619CFF"))
gp6 <- trajplot(fit=fit.BCC2,feature.ind=3,which.cluster = "global.cluster",
title="Global Clustering",
xlab="months",ylab="spiders",color=c("#00BA38", "#619CFF"))
library(cowplot)
#dev.new(width=180, height=120)
plot_grid(gp1, gp2,gp3,gp4,gp5,gp6,
labels=c("(A)", "(B)", "(C)", "(D)", "(E)", "(F)"), ncol = 3, align = "v" )
Package References
sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] mixAK_5.8 lme4_1.1-35.5 Matrix_1.7-1 colorspace_2.1-1
#> [5] cowplot_1.1.3 ggplot2_3.5.1 joineRML_0.4.6 survival_3.7-0
#> [9] nlme_3.1-166 BCClong_1.0.3 rmarkdown_2.29
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 xfun_0.49 bslib_0.8.0
#> [4] lattice_0.22-6 vctrs_0.6.5 tools_4.4.2
#> [7] generics_0.1.3 stats4_4.4.2 parallel_4.4.2
#> [10] tibble_3.2.1 fansi_1.0.6 cluster_2.1.6
#> [13] pkgconfig_2.0.3 lifecycle_1.0.4 farver_2.1.2
#> [16] compiler_4.4.2 MatrixModels_0.5-3 mcmc_0.9-8
#> [19] munsell_0.5.1 mnormt_2.1.1 combinat_0.0-8
#> [22] fastGHQuad_1.0.1 codetools_0.2-20 SparseM_1.84-2
#> [25] quantreg_5.99 htmltools_0.5.8.1 sys_3.4.3
#> [28] buildtools_1.0.0 sass_0.4.9 evd_2.3-7.1
#> [31] yaml_2.3.10 gmp_0.7-5 pillar_1.9.0
#> [34] nloptr_2.1.1 jquerylib_0.1.4 MASS_7.3-61
#> [37] randtoolbox_2.0.5 truncdist_1.0-2 cachem_1.1.0
#> [40] iterators_1.0.14 foreach_1.5.2 boot_1.3-31
#> [43] abind_1.4-8 mclust_6.1.1 tidyselect_1.2.1
#> [46] digest_0.6.37 mvtnorm_1.3-2 LaplacesDemon_16.1.6
#> [49] dplyr_1.1.4 labeling_0.4.3 maketools_1.3.1
#> [52] splines_4.4.2 fastmap_1.2.0 grid_4.4.2
#> [55] cli_3.6.3 magrittr_2.0.3 cobs_1.3-8
#> [58] label.switching_1.8 utf8_1.2.4 withr_3.0.2
#> [61] Rmpfr_1.0-0 scales_1.3.0 rngWELL_0.10-10
#> [64] nnet_7.3-19 gridExtra_2.3 coda_0.19-4.1
#> [67] evaluate_1.0.1 lpSolve_5.6.22 knitr_1.49
#> [70] doParallel_1.0.17 mgcv_1.9-1 rlang_1.1.4
#> [73] MCMCpack_1.7-1 Rcpp_1.0.13-1 glue_1.8.0
#> [76] minqa_1.2.8 jsonlite_1.8.9 R6_2.5.1