MixedTypeData"
In BCClong: Bayesian Consensus Clustering for Multiple Longitudinal Features

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

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:

ls("package:BCClong")

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(\bm{\y}, \bm{\Theta}) 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.

library(BCClong)
library(mixAK)
data(PBC910)

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 Ryzen$^{TM}$ 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.

# pre-compiled result
data(PBCseqfit)
fit.BCC2 <- PBCseqfit

Printing Summary Statistics for key model parameters

To print the BCC model

print(fit.BCC2)

To print the summary statistics for all parameters

summary(fit.BCC2)

To print the proportion \pi for each cluster (mean, sd, 2.5% and 97.5% percentile) geweke statistics (geweke.stat) between -2 and 2 suggests the parameters converge

fit.BCC2$summary.stat$PPI

The code below prints out all major parameters

summary(fit.BCC2)

Visualize Clusters

Generic plot can be used on BCC object, all relevant plots will be generate one by one using return key

plot(fit.BCC2)

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" )