This document walks through example simulation code used to produce the results obtained in ‘Penalized model-based clustering of fMRI’ submitted to the journal Biostatistics. First, we install and load the version of the R package used for the submitted manuscript:

Simulating Data

We simulate an example data set for \(K=100\) total subjects belong to one of \(G=2\) equally sized clusters which share roughly \(20 \%\) of network connections. For each subject, we generate \(n=177\) observations for \(p=10\) variables with approximately \(\rho = 10\%\) of network connections differing from their true cluster-level network.

Tuning Parameter Selection

We now select the optimal sets of tuning parameters for 3 different analysis approaches: our random covariance clustering model (RCCM), a Ward clustering & group graphical lasso (GGL) approach, and a graphical lasso (GLasso) & K-means clustering approach. The scale of which tuning parameters work best for each respective method vary, although in practice we’ve found that the optimal lasso penalty parameter for RCCM, \(\lambda_1\), tends to be \(\approx 100\) times the magnitude of what works well for GGL and GLasso. Similarly, we found that the optimal tuning parameter for RCCM that encourages similarity within clusters, denoted by \(\lambda_2\), tends to be \(\approx 5000\) times the magnitude of what works well for GGL when \(p=10\), but this did vary by the number of variables. In other words, if \(\lambda_1 = 50\) works well for RCCM for a given data set, then \(\lambda_1 = 0.50\) tends to work well for GLasso and GGL, and if \(\lambda_2 = 5000\) works well for RCCM, then \(\lambda_2 = 1\) works fairly well for GGL.

Analysis

We then analyze our simulated data set using the optimally selected tuning parameters for each method.

Ward & GGL

# Using Ward & GGL approach

# Function for calculating pair-wise Frob norm differences and then clustering
MatClust <- function(mats, G) {
  K <- dim(mats)[3]
  combos <- expand.grid(s1 = 1:K, s2 = 1:K)
  distMat <- matrix(NA, nrow = K, ncol = K)
  for(r in 1:K) {
    for(s in 1:K) {
      distMat[r, s] <- norm(mats[, , r] - mats[, , s], type = 'F')
    }
  }

  cl0 <- cutree(hclust(d = as.dist(distMat), method = "ward.D"), k = G)

  wgk <- matrix(NA, nrow = G, ncol = K)

  for(i in 1:G) {
    for(j in 1:K) {
      wgk[i, j] <- ifelse(cl0[j] == i, 1, 0)
    }
  }
  return(wgk)
}

K <- length(myData$simDat)
Sl <- sapply(myData$simDat, cov, simplify = "array")
gl <- sapply(1:K, simplify = "array", FUN = function(x){glasso::glasso(Sl[, , x], rho = 0.001,
                                                                       penalize.diagonal = FALSE)$wi})

GGLResults <- lapply(X = c("stARS", "CV"), FUN = function(tune) {
  
lambda1 <- ifelse(tune == "stARS", optWardggl$lambda1[1], 
                    optTuneCV[which(optTuneCV$method == "GGL"), ]$lambda1)

lambda2 <- ifelse(tune == "stARS", optWardggl$lambda2[1], 
                  optTuneCV[which(optTuneCV$method == "GGL"), ]$lambda2)
  
# Estimating cluster memberships using Ward clustering based on dissimilarity matrix of Frob norm differences
GGLres <- list()
GGLres$weights <- MatClust(gl, G = G)

# Analyzing using GGL within each estimated cluster
GGLres$res <- unlist(lapply(FUN = function(g) {
  prec <- JGL::JGL(Y = myData$simDat[which(as.logical(GGLres$weights[g, ]))], penalty = "group",
                   penalize.diagonal = FALSE, lambda1 = lambda1 / 100,
                   lambda2 = lambda2 / 50000, return.whole.theta = TRUE)$theta
  return(setNames(prec, c(which(as.logical(GGLres$weights[g, ])))))}, X = 1:G), recursive = F)

GGLzHat <- apply(GGLres$weights, MARGIN = 2, FUN = which.max)

# Estimated group-level network for GGL. Edge present if >= 50% of subjects have it in group
  GGLres$Omegags <- array(NA, dim = c(p, p, G))
  for(g in 1:G) {
    GGLres$Omegags[, , g] <- round(apply(simplify2array(lapply(GGLres$res[which(GGLzHat == g)], FUN = adj)),
                                  MARGIN = 1:2, FUN = mean))
  }
  
  return(GGLres)})

Results

We summarize the performance of each method using both tuning parameter selection approaches. For evaluating clustering performance, we calculate the the rand index (RI) and adjusted rand index (RIadj). For evaluating performance in terms of edge detection we calculate the true positive rate (TPR), false positive rate (FPR), precision, and F1 scores for the subject-level networks, subscripted with \(k\), and cluster-level networks, subscripted with \(g\).

# Function for summarizing edge detection performances
performanceSummary <- function(Omegaks, Omegags, zs, Omega0ks = myData$Omegaks, Omega0gs = myData$Omega0s, z0s = myData$zgks) {

# Calculating Rand indexes and adjusted rand indexes
RI <- randCalc(zs, z0s)
RIadj <- mclust::adjustedRandIndex(zs, z0s)

# Calculating Precision Matrix Error, True Positive Rate, and False Positive Rates
subjSum <- sum(sapply(1:K, FUN = function(k){sum((adj(Omegaks[, , k]) +
                                                    adj(Omega0ks[, , k]))[lower.tri(Omega0ks[, , k], diag = FALSE)] == 2)}))
posskEdges <- sum(sapply(1:K, FUN = function(k){sum(adj(Omega0ks[, , k])[lower.tri(Omega0ks[, , k], diag = FALSE)] == 1)}))

TPRk <- subjSum / posskEdges
subjSum0 <- sum(sapply(1:K, FUN = function(k){sum((-1*adj(Omegaks[, , k]) +
                                                     adj(Omega0ks[, , k]))[lower.tri(Omega0ks[, , k], diag = FALSE)] == -1)}))
possk0s <- sum(sapply(1:K, FUN = function(k){sum(adj(Omega0ks[, , k])[lower.tri(Omega0ks[, , k], diag = FALSE)] == 0)}))

FPRk <- subjSum0 / possk0s

PrecisionK <- subjSum / (subjSum + subjSum0)

F1k <- 2*(PrecisionK*TPRk) / (PrecisionK + TPRk)

if(is.null(Omegags) == FALSE) {
  grpSum <- sum(sapply(1:G, FUN = function(g){sum((adj(Omegags[, , g]) +
                                                     adj(Omega0gs[, , g]))[lower.tri(Omega0gs[, , g], diag = FALSE)] == 2)}))
  possgEdges <- sum(sapply(1:G, FUN = function(g){sum(adj(Omega0gs[, , g])[lower.tri(Omega0gs[, , g], diag = FALSE)] == 1)}))
  TPRg <- grpSum / possgEdges
  grpSum0 <- sum(sapply(1:G, FUN = function(g){sum((-1*adj(Omegags[, , g]) +
                                                      adj(Omega0gs[, , g]))[lower.tri(Omega0gs[, , g], diag = FALSE)] == -1)}))
  possg0s <- sum(sapply(1:G, FUN = function(g){sum(adj(Omega0gs[, , g])[lower.tri(Omega0gs[, , g], diag = FALSE)] == 0)}))

  PrecisionG <- grpSum / (grpSum + grpSum0)

  F1g <- 2*(PrecisionG*TPRg) / (PrecisionG + TPRg)

  FPRg <- grpSum0 / possg0s
} else {
  TPRg <- NA
  FPRg <- NA
  PrecisionG <- NA
  F1g <- NA
}

return(data.frame(RI = RI, RIadj = RIadj, TPRk = TPRk, FPRk = FPRk,
                  PrecisionK = PrecisionK, F1k = F1k,
                  TPRg = TPRg, FPRg = FPRg, PrecisionG = PrecisionG, F1g = F1g))
}