R/bandle-function.R

In ococrook/bandle: An R package for the Bayesian analysis of differential subcellular localisation experiments

##' These function implement the bandle model for dynamic mass spectrometry based
##' spatial proteomics datasets using MCMC for inference
##' 
##' The `bandle` function generate the sample from the posterior distributions
##' (object or class `bandleParams`) based on an annotated quantitative spatial
##' proteomics datasets (object of class [`MSnbase::MSnSet`]). Both are then
##' passed to the `bandlePredict` function to predict the sub-cellular localisation
##' and compute the differential localisation probability of proteins. See the
##' vignette for examples
##' 
##' @title Differential localisation experiments using the bandle method
##' @param objectCond1 A list of [`MSnbase::MSnSet`]s where each is an experimental
##' replicate for the first condition, usually a control
##' @param objectCond2 A list of [`MSnbase::MSnSet`]s where each is an experimental
##' replicate for the second condition, usually a treatment
##' @param fcol The feature meta-data containing marker definitions. Default is
##' `markers`
##' @param hyperLearn Algorithm to learn posterior hyperparameters of the Gaussian
##' processes. Default is `LBFGS` and `MH` for metropolis-hastings is also implemented. 
##' @param numIter The number of iterations of the MCMC
##'     algorithm. Default is 1000. Though usually much larger numbers are used
##' @param burnin The number of samples to be discarded from the
##'     beginning of the chain. Default is 100.
##' @param thin The thinning frequency to be applied to the MCMC
##'     chain.  Default is 5.
##' @param u The prior shape parameter for Beta(u, v). Default is 2
##' @param v The prior shape parameter for Beta(u, v). Default is 10.
##' @param lambda Controls the variance of the outlier component. Default is 1.
##' @param gpParams Parameters from prior fitting of GPs to each niche to accelerate 
##' inference. If used, must be of class `gpParams`. Default is NULL.
##' @param hyperIter The frequency of MCMC interation to update the hyper-parameters
##' default is 20
##' @param hyperMean The prior mean of the log normal prior of the GP parameters.
##' Default is 0 for each. Order is length-scale, amplitude and noise variance
##' @param hyperSd The prior standard deviation of the log normal prior fo the GP
##' parameters. Default is 1 for each. Order is length-scale, ampliture and noise
##'  variance.
##' @param seed The random number seed. This is set internally if not provided.
##' For those doing comparison experiments using simulations must set this
##' parameter themselves, else the results will not be representative of the 
##' performance of the package.
##' @param pg `logical` indicating whether to use polya-gamma prior. Default is
##'  `FALSE`.
##' @param pgPrior A matrix generated by pgPrior function. If param pg is TRUE
##' but pgPrior is NULL then a pgPrior is generated on the fly.
##' @param tau The tau parameter for the polya-Gamma prior (is used). Defaults
##' to 0.2
##' @param dirPrior A matrix generated by dirPrior function. Default is NULL
##'  and dirPrior is generated on the fly. This should have dimension
##'  equal to the number of subcellular niches. We suggest referring
##'  to the vignette on direction to set this value.
##' @param maternCov `logical` indicated whether to use a matern or Gaussian
##' covariance. Default is True and matern covariance is used
##' @param PC `logical` indicating whether to use a penalised complexity prior.
##' Default is TRUE.
##' @param pcPrior `matrix` of length 3 indicating the lambda parameters for the
##' penalised complexity prior. Default is `matrix(c(0.5, 3, 100), nrow = 1)` and the order is 
##' length-scale, amplitude and variance.
##' @param nu `integter` indicating the smoothness of the matern prior. Default
##' is 2.
##' @param propSd If MH is used to learn posterior hyperparameters then the proposal
##' standard deviations. A Gaussian random-walk proposal is used.
##' @param numChains `integer` indicating the number of parallel chains to run.
##' Defaults to 4.
##' @param BPPARAM BiocParallel parameter. Defaults to machine default backend
##' using bpparam()
##' @return  `bandle` returns an instance of class `bandleParams`
##' @md
##' 
##' @examples
##' library(pRolocdata)
##' data("tan2009r1")
##' set.seed(1)
##' tansim <- sim_dynamic(object = tan2009r1, 
##'                       numRep = 6L,
##'                       numDyn = 100L)
##' gpParams <- lapply(tansim$lopitrep, function(x) 
##' fitGPmaternPC(x, hyppar = matrix(c(0.5, 1, 100), nrow = 1)))
##' d1 <- tansim$lopitrep
##' control1 <- d1[1:3]
##' treatment1 <- d1[4:6]
##' mcmc1 <- bandle(objectCond1 = control1,
##'  objectCond2 = treatment1, gpParams = gpParams,
##'  fcol = "markers", numIter = 5L, burnin = 1L, thin = 2L,
##'  numChains = 1, BPPARAM = SerialParam(RNGseed = 1))
##' 
##' @rdname bandle
bandle <- function(objectCond1,
                   objectCond2,
                   fcol = "markers",
                   hyperLearn = "LBFGS",
                   numIter = 1000,
                   burnin = 100L,
                   thin = 5L,
                   u = 2,
                   v = 10,
                   lambda = 1,
                   gpParams = NULL,
                   hyperIter = 20,
                   hyperMean = c(0, 0, 0),
                   hyperSd = c(1, 1,1),
                   seed = NULL,
                   pg = FALSE,
                   pgPrior = NULL,
                   tau = 0.2,
                   dirPrior = NULL,
                   maternCov = TRUE,
                   PC = TRUE,
                   pcPrior = matrix(c(0.5,3,100), nrow = 1),
                   nu = 2,
                   propSd = c(0.3, 0.1, 0.05),
                   numChains = 4L,
                   BPPARAM = BiocParallel::bpparam()){
    
    
    # Checks
    stopifnot(exprs = {
                "ObjectCond1 must be a list of MSnSet"=is(objectCond1[[1]], "MSnSet")
                "ObjectCond2 must be a list of  MSnSet"=is(objectCond2[[1]], "MSnSet")
                "hyperLearn must be either MH or LBFGS"=hyperLearn %in% c("MH", "LBFGS")
                "numIter must be a numeric"=is(numIter, "numeric")
                "burnin must be an integer"=is(burnin, "integer")
                "thin must be an integer"=is(thin, "integer")
                "burnin must be less than numIter"= burnin < numIter
                "u must be numeric"=is(u, "numeric")
                "v must be numeric"=is(v, "numeric")
                "lambda must be numeric"=is(lambda, "numeric")
                "hyperIter must be numeric"=is(hyperIter, "numeric")
                "hyperMean must be numeric"=is(hyperMean, "numeric")
                "hyperSd must be numeric"=is(hyperSd, "numeric")
                "must provide 3 values for hyperMean"=length(hyperMean) == 3
                "must provide 3 values for hyperSd" = length(hyperSd) == 3
                "tau must be numeric" = is(tau, "numeric")
                "nu must be numeric" = is(nu, "numeric")
                "propSd must be numeric"=is(propSd, "numeric")
                "Must provide 3 values for propSd"=length(propSd) == 3})
    
    # if dirPrior is not NULL
    if(!is.null(dirPrior)){
        K <- length(getMarkerClasses(objectCond1[[1]], fcol = fcol))
        stopifnot("dirPrior must have dimensions equal to the number of
                  niches"=dim(dirPrior)==c(K, K))
    }
  
    # if gpParams if not NULL
    if(!is.null(gpParams)){
      nP <- length(gpParams)
      nD <- length(objectCond1) + length(objectCond2)
      stopifnot("gpParams object must be the same length as the total numer of
                  datasets"= nP == nD)
    }

    # valid experiment
    ## same number of rows
    validrow <- length(unique(vapply(c(objectCond1, objectCond2),
                                     function(x) nrow(x),
                                     FUN.VALUE = numeric(1)))) > 1
    if (isTRUE(validrow)){
        stop("Number of rows do not match, you may wish to subset your proteins")
    }
    ## same number of columns
    validcol <- length(unique(vapply(c(objectCond1, objectCond2),
                                     function(x) ncol(x),
                                     FUN.VALUE = numeric(1)))) > 1
    if (isTRUE(validcol)){
        stop("Number of columns do not match, this analysis is not currently
             implemented in bandle. Subset fractions so they match")
    }
    validrow2 <- var(apply(vapply(c(objectCond1, objectCond2),
                                  function(x) rownames(x),
                                  character(nrow(objectCond1[[1]]))), 1,
                           function(x) length(unique(x)))) > 0
    if (isTRUE(validrow2)){
        stop("rownames do not match across experiments. Analysis
             will lead to confusing results. Subset your proteins
             so they match")
    }
    ## valid fcol
    if (isFALSE(all(vapply(c(objectCond1, objectCond2), function(x)
        fcol %in% colnames(fData(x)), logical(1))))){
        stop("fcol is not in all the datasets. Check fcol is present")
    }
    
    ## chains run in parallel, repeating number of iterations
    .res <- BiocParallel::bplapply(rep(numIter, numChains),
                                   FUN = diffLoc,
                                   objectCond1 = objectCond1,
                                   objectCond2 = objectCond2,
                                   fcol = fcol,
                                   hyperLearn = hyperLearn,
                                   #numIter = numIter, #parrallelised 
                                   burnin = burnin,
                                   thin = thin,
                                   u = u,
                                   v = v,
                                   lambda = lambda,
                                   gpParams = gpParams,
                                   hyperIter = hyperIter,
                                   hyperMean = hyperMean,
                                   hyperSd = hyperSd,
                                   seed = seed,
                                   pg = pg,
                                   pgPrior = pgPrior,
                                   tau = tau,
                                   dirPrior = dirPrior,
                                   maternCov = maternCov,
                                   PC = PC,
                                   pcPrior = pcPrior,
                                   nu = nu,
                                   propSd = propSd,
                                   BPPARAM = BPPARAM)
    
    K <- length(getMarkerClasses(objectCond1[[1]], fcol = fcol))
    
    # construct empirical Bayes Polya-Gamma prior
    if (is.null(pgPrior)) {
        pgPrior <- pg_prior(objectCond1, objectCond2, K = K, pgPrior = NULL, fcol = fcol)
    }
    
    #or Dirichlet weights Prior, unless dirPrior is provided by default
    if(pg == FALSE){
        if (is.null(dirPrior)){
            dirPrior <- diag(rep(1, K)) + matrix(0.05, nrow = K, ncol = K)
            rownames(dirPrior) <- colnames(dirPrior) <- getMarkerClasses(objectCond1[[1]],
                                                                       fcol = fcol)
        }
    }
    
    # name objects
    colnames(pcPrior) <- c("length-scale", "amplitude", "sigma")
    
    ## Construct class bandleChains
    .ans <- .bandleChains(chains = .res)
    
    ## Construct class bandleParams
    .out <- .bandleParams(method = "bandle",
                          priors = list(dirPrior = dirPrior,
                                        pgPrior = pgPrior,
                                        pcPrior = pcPrior,
                                        hyperMean = hyperMean,
                                        hyperSd = hyperSd), 
                          seed = as.integer(seed), 
                          summary = .bandleSummaries(
                                        summaries = list(.bandleSummary(),
                                                         .bandleSummary())),
                          chains = .ans)
    
    return(.out)
}

##' @title Make predictions from a bandle analysis
##' @param objectCond1 A list of instances of class [`MSnbase::MSnSet`]s
##' where each is an experimental replicate for the first condition, usually a control
##' @param objectCond2 A list of instance of class [`MSnbase::MSnSet`]s 
##' where each is an experimental replicate for the second condition, usually a treatment
##' @param params An instance of class `bandleParams`, as generated by
##'     [`bandle()`].
##' @param fcol A feature column indicating the markers. Defaults to "markers"    
##' @return `bandlePredict` returns an instance of class
##'     [`MSnbase::MSnSet`] containing the localisation predictions as
##'     a new `bandle.allocation` feature variable. The allocation
##'     probability is encoded as `bandle.probability`
##'     (corresponding to the mean of the distribution
##'     probability). In addition the upper and lower quantiles of
##'     the allocation probability distribution are available as
##'     `bandle.probability.lowerquantile` and
##'     `bandle.probability.upperquantile` feature variables. The
##'     Shannon entropy is available in the `bandle.mean.shannon`
##'     feature variable, measuring the uncertainty in the allocations
##'     (a high value representing high uncertainty; the highest value
##'     is the natural logarithm of the number of classes). An additional variable
##'     indicating the differential localization probability is also added as
##'     `bandle.differential.localisation`
##'     
##' @examples 
##' library(pRolocdata)
##' data("tan2009r1")
##' set.seed(1)
##' tansim <- sim_dynamic(object = tan2009r1, 
##'                     numRep = 6L,
##'                    numDyn = 100L)
##' gpParams <- lapply(tansim$lopitrep, function(x) 
##' fitGPmaternPC(x, hyppar = matrix(c(0.5, 1, 100), nrow = 1)))
##' d1 <- tansim$lopitrep
##' control1 <- d1[1:3]
##' treatment1 <- d1[4:6]
##' mcmc1 <- bandle(objectCond1 = control1, objectCond2 = treatment1, gpParams = gpParams,
##'                                      fcol = "markers", numIter = 5L, burnin = 1L, thin = 2L,
##'                                      numChains = 1, BPPARAM = SerialParam(RNGseed = 1))
##' mcmc1 <- bandleProcess(mcmc1)
##' out <- bandlePredict(objectCond1 = control1, objectCond2 = treatment1, params = mcmc1)
##' @md
##' @rdname bandle-predict
bandlePredict <- function(objectCond1,
                          objectCond2,
                          params,
                          fcol = "markers"){

    stopifnot(inherits(params, "bandleParams"))
    ## Checks for object and bandleParams match
    stopifnot(featureNames(unknownMSnSet(objectCond1[[1]], fcol = fcol))
          == rownames(params@summary@summaries[[1]]@posteriorEstimates))
    stopifnot(unlist(lapply(objectCond1, function(zz) inherits(zz, "MSnSet"))))
    stopifnot(unlist(lapply(objectCond2, function(zz) inherits(zz, "MSnSet"))))

    for (j in c(1, 2)){
    
        object <- list(objectCond1, objectCond2)[[j]][[1]]
    
        ## Create marker set and size
        markerSet <- markerMSnSet(object, fcol = fcol)
        markers <- getMarkerClasses(object, fcol = fcol)
        M <- nrow(markerSet)
        K <- chains(params)[[1]]@K
    
    
        ## Get Summary object from bandleParams maybe better to check
        ## columns exist/pass which objects we need
        .bandle.allocation <- c(as.character(summaries(params)[[j]]@posteriorEstimates[,"bandle.allocation"]),
                              as.character(fData(markerSet)[, fcol]))
        .bandle.probability <- c(summaries(params)[[j]]@posteriorEstimates[,"bandle.probability"],
                               rep(1, M)) ## set all probabilities of markers to 1.
        .bandle.probability.lowerquantile <- c(summaries(params)[[j]]@posteriorEstimates[,"bandle.probability.lowerquantile"],
                                             rep(1, M)) ## set all probabilities of markers to 1.
        .bandle.probability.upperquantile <- c(summaries(params)[[j]]@posteriorEstimates[,"bandle.probability.upperquantile"],
                                             rep(1, M)) ## set all probabilities of markers to 1.
        .bandle.mean.shannon <- c(summaries(params)[[j]]@posteriorEstimates[,"bandle.mean.shannon"],
                                rep(0, M)) ## set all probabilities of markers to 0
        .bandle.differential.localisation <- c(summaries(params)[[j]]@posteriorEstimates[, "bandle.differential.localisation"],
                                               rep(0, M)) ## set differential localisation of markers to 0
        
        ## Create data frame to store new summaries
        .bandle.summary <- DataFrame(bandle.allocation = .bandle.allocation ,
                                    bandle.probability = .bandle.probability,
                                    bandle.probability.lowerquantile = .bandle.probability.lowerquantile,
                                    bandle.probability.upperquantile = .bandle.probability.upperquantile,
                                    bandle.mean.shannon = .bandle.mean.shannon,
                                    bandle.differential.localisation = .bandle.differential.localisation)
        ## add outlier information
        .bandle.summary$bandle.outlier <- c(summaries(params)[[j]]@posteriorEstimates[, "bandle.outlier"],
              rep(0, M)) ## set all probabilities of markers to 0
        
        
        ## Check number of rows match and add feature names
        stopifnot(nrow(.bandle.summary) == nrow(object))
        rownames(.bandle.summary) <- c(rownames(summaries(params)[[j]]@posteriorEstimates),
                                     rownames(markerSet))
        
        ## Append data to fData of MSnSet
        fData(object) <- cbind(fData(object), .bandle.summary[rownames(fData(object)),])
        
        ## create allocation matrix for markers
        .probmat <- matrix(0, nrow = nrow(markerSet), ncol = K)
        .class <- fData(markerSet)[, fcol]
        for (k in seq_len(nrow(markerSet))) {
        ## give markers prob 1
            .probmat[k, as.numeric(factor(.class), seq(1, length(unique(.class))))[k]] <- 1
        }
        colnames(.probmat) <- markers
        rownames(.probmat) <- rownames(markerSet)
        .joint <- rbind(summaries(params)[[j]]@bandle.joint, .probmat)
        fData(object)$bandle.joint <- .joint[rownames(fData(object)), ]
    
        if(j == 1){
            objectCond1[[1]] <- object
        } else {
            objectCond2[[1]] <- object
        }
    }
    
    .out <- list(objectCond1 = objectCond1,
                 objectCond2 = objectCond2)
    
    return(.out)
}
##' @title process bandle results
##' @param params An object of class `bandleParams`
##' @return `bandleProcess` returns an instance of class
##'     `bandleParams` with its summary slot populated.
##'     
##' @examples 
##' library(pRolocdata)
##' data("tan2009r1")
##' set.seed(1)
##' tansim <- sim_dynamic(object = tan2009r1, 
##'                     numRep = 6L,
##'                    numDyn = 100L)
##' gpParams <- lapply(tansim$lopitrep, function(x) 
##' fitGPmaternPC(x, hyppar = matrix(c(0.5, 1, 100), nrow = 1)))
##' d1 <- tansim$lopitrep
##' control1 <- d1[1:3]
##' treatment1 <- d1[4:6]
##' mcmc1 <- bandle(objectCond1 = control1, objectCond2 = treatment1, gpParams = gpParams,
##'                                      fcol = "markers", numIter = 5L, burnin = 1L, thin = 2L,
##'                                      numChains = 1, BPPARAM = SerialParam(RNGseed = 1))
##' mcmc1 <- bandleProcess(mcmc1)
##' 
##' @md
##' @rdname bandle-process
bandleProcess <- function(params) {
    
    stopifnot(is(params, "bandleParams"))
    
    ## get require slots
    myChain <- chains(params)[[1]]
    numChains <- length(chains(params))
    N <- myChain@N
    K <- myChain@K
    n <- myChain@n
    
    ## storage
    meanComponentProbCond1 = vector("list", length = numChains)
    meanComponentProbCond2 = vector("list", length = numChains)
    meanOutlierProb = vector("list", length = numChains)
    bandle.jointCond1 <- matrix(0, nrow = N, ncol = K)
    bandle.jointCond2 <- matrix(0, nrow = N, ncol = K)
    bandle.outlierCond1 <- matrix(0, nrow = N, ncol = 1)
    bandle.outlierCond2 <- matrix(0, nrow = N, ncol = 1)
    .bandle.quantilesCond1 <- array(0, c(N, K, 2))
    .bandle.quantilesCond2 <- array(0, c(N, K, 2))
    pooled.ComponentCond1 <- matrix(0, nrow = N, ncol = n * numChains)
    pooled.ComponentProbCond1 <- array(0, c(N, n * numChains, K ))
    pooled.OutlierCond1 <- matrix(0, nrow = N, ncol = n * numChains)    
    pooled.ComponentCond2 <- matrix(0, nrow = N, ncol = n * numChains)
    pooled.ComponentProbCond2 <- array(0, c(N, n * numChains, K ))
    pooled.OutlierCond2 <- matrix(0, nrow = N, ncol = n * numChains)
    pooled.differential.localisation <- matrix(0, nrow = N, ncol = numChains)
    
    # Calculate basic quantities
    for (j in seq_len(numChains)) {
        
        mc <- chains(params)[[j]]
        meanComponentProbCond1[[j]] <- t(apply(mc@nicheProb[[1]][, , ], 1, colMeans))
        meanComponentProbCond2[[j]] <- t(apply(mc@nicheProb[[2]][, , ], 1, colMeans))
        bandle.jointCond1 <- bandle.jointCond1 + meanComponentProbCond1[[j]]
        bandle.jointCond2 <- bandle.jointCond2 + meanComponentProbCond2[[j]]
        bandle.outlierCond1 <- bandle.outlierCond1 + 1 - rowMeans(mc@outlier[[1]])
        bandle.outlierCond2 <- bandle.outlierCond2 + 1 - rowMeans(mc@outlier[[2]])
        ## Pool chains
        pooled.ComponentCond1[, n * (j - 1) + seq.int(n)] <- mc@niche[[1]]
        pooled.ComponentProbCond1[, n * (j - 1) + seq.int(n), ] <- mc@nicheProb[[1]]
        pooled.OutlierCond1[, n * (j - 1)+ seq.int(n)] <- mc@outlier[[1]]
        pooled.ComponentCond2[, n * (j - 1) + seq.int(n)] <- mc@niche[[2]]
        pooled.ComponentProbCond2[, n * (j - 1) + seq.int(n), ] <- mc@nicheProb[[2]]
    }
    ## take means across chains
    bandle.jointCond1 <- bandle.jointCond1/numChains
    bandle.outlierCond1 <- bandle.outlierCond1/numChains
    bandle.probabilityCond1 <- apply(bandle.jointCond1, 1, max)
    bandle.allocationCond1 <- colnames(bandle.jointCond1)[apply(bandle.jointCond1, 1, which.max)]
    bandle.jointCond2 <- bandle.jointCond2/numChains
    bandle.outlierCond2 <- bandle.outlierCond2/numChains
    bandle.probabilityCond2 <- apply(bandle.jointCond2, 1, max)
    bandle.allocationCond2 <- colnames(bandle.jointCond2)[apply(bandle.jointCond2, 1, which.max)]
    bandle.differential.localisation <- diffLocalisationProb(params = params)
    
    ## Calculate quantiles
    for (i in seq_len(N)) {
        for (j in seq_len(K)) {
            .bandle.quantilesCond1[i, j, ] <- quantile(pooled.ComponentProbCond1[i, , j],
                                                    probs = c(0.025, 0.975))
            .bandle.quantilesCond2[i, j, ] <- quantile(pooled.ComponentProbCond2[i, , j],
                                                       probs = c(0.025, 0.975))
        }
    }
    
    ## Store quantiles
    bandle.probability.lowerquantile.cond1 <- .bandle.quantilesCond1[cbind(seq.int(N),
                                                                apply(bandle.jointCond1, 1, which.max), rep(1, N))]
    bandle.probability.upperquantile.cond1 <- .bandle.quantilesCond1[cbind(seq.int(N),
                                                                apply(bandle.jointCond1, 1, which.max), rep(2, N))]
    bandle.probability.lowerquantile.cond2 <- .bandle.quantilesCond2[cbind(seq.int(N),
                                                                apply(bandle.jointCond2, 1, which.max), rep(1, N))]
    bandle.probability.upperquantile.cond2 <- .bandle.quantilesCond2[cbind(seq.int(N),
                                                                apply(bandle.jointCond2, 1, which.max), rep(2, N))]
    
    ## Compute Shannon Entropy
    bandle.shannonCond1 <- -apply(pooled.ComponentProbCond1 * log(pooled.ComponentProbCond1), c(1,2), sum)
    bandle.shannonCond2 <- -apply(pooled.ComponentProbCond2 * log(pooled.ComponentProbCond2), c(1,2), sum)
    bandle.shannonCond1[is.na(bandle.shannonCond1)] <- 0
    bandle.mean.shannonCond1 <- rowMeans(bandle.shannonCond1)
    bandle.mean.shannonCond2 <- rowMeans(bandle.shannonCond2)
    
    ## Name entries
    names(bandle.mean.shannonCond1) <- names(bandle.probability.lowerquantile.cond1) <-
        names(bandle.probability.upperquantile.cond1) <- names(bandle.mean.shannonCond2) <-
        names(bandle.probability.lowerquantile.cond1) <-
        names(bandle.probability.upperquantile.cond2) <- rownames(myChain@niche[[1]])
    rownames(bandle.jointCond1) <- names(bandle.probabilityCond1) <-
        names(bandle.allocationCond1) <- rownames(bandle.jointCond2) <- names(bandle.probabilityCond2) <-
        names(bandle.allocationCond2) <- rownames(myChain@niche[[1]])
    colnames(bandle.jointCond1) <- colnames(bandle.jointCond2) <- dimnames(myChain@nicheProb[[1]])[[3]]
    
    ## Outlier colNames # not use

    ## Compute convergence diagnostics
    weights <- vector("list", length = numChains)
    r_hat <- vector(mode = "list", length = K)
    
    for(j in seq_len(numChains)) {
        mc <- chains(params)[[j]]
        weights[[j]] <- coda::mcmc(mc@weights[1,1,])
    }
    
    ## Summary of posterior estimates
    bandle.summary1 <- DataFrame(bandle.allocation = bandle.allocationCond1,
                                 bandle.probability = bandle.probabilityCond1,
                                 bandle.outlier = bandle.outlierCond1,
                                 bandle.probability.lowerquantile = bandle.probability.lowerquantile.cond1,
                                 bandle.probability.upperquantile = bandle.probability.upperquantile.cond1,
                                 bandle.mean.shannon = bandle.mean.shannonCond1,
                                 bandle.differential.localisation = bandle.differential.localisation)
    
    colnames(bandle.summary1) <- c("bandle.allocation", "bandle.probability",
                                   "bandle.outlier","bandle.probability.lowerquantile", "bandle.probability.upperquantile",
                                   "bandle.mean.shannon", "bandle.differential.localisation")
    
    
    
    bandle.summary2 <- DataFrame(bandle.allocation = bandle.allocationCond2,
                                 bandle.probability = bandle.probabilityCond2,
                                 bandle.outlier = bandle.outlierCond2,
                                 bandle.probability.lowerquantile = bandle.probability.lowerquantile.cond2,
                                 bandle.probability.upperquantile = bandle.probability.upperquantile.cond2,
                                 bandle.mean.shannon = bandle.mean.shannonCond2,
                                 bandle.differential.localisation = bandle.differential.localisation)
    
    colnames(bandle.summary2) <- c("bandle.allocation", "bandle.probability",
                                   "bandle.outlier","bandle.probability.lowerquantile", "bandle.probability.upperquantile",
                                   "bandle.mean.shannon", "bandle.differential.localisation")
    
    .diagnostics <- matrix(NA, 1, 1)
    ## Compute convergence diagnostics
    if(numChains > 1){
        weights <- coda::as.mcmc.list(weights)
        gd <- coda::gelman.diag(x = weights, autoburnin = FALSE)
        .diagnostics <- gd$psrf
        rownames(.diagnostics) <- c("weights")
        
    }

    ## Constructor for summary
    params@summary@summaries[[1]] <- .bandleSummary(posteriorEstimates = bandle.summary1,
                                                    diagnostics = .diagnostics,
                                                    bandle.joint = bandle.jointCond1)
    params@summary@summaries[[2]] <- .bandleSummary(posteriorEstimates = bandle.summary2,
                                                    diagnostics = .diagnostics,
                                                    bandle.joint = bandle.jointCond2)
    
    return(params)
}