R/testGlasso.R

Defines functions testGlasso

Documented in testGlasso

#' @name testGlasso
#' @title Test for for the equality of connectivity based on the Graphical lasso estimation.
#' @description This function utilizes Dynamic Connectivity Regression (DCR) algorithm proposed by Cribben el al. (2012) to test the equality of connectivity in two fMRI signals.
#' @param subY1 a sample of size length*dim
#' @param subY2 a sample of size length*dim
#' @param p  Gep(p) distribution controls the size of stationary bootstrap. The mean block length is 1/p
#' @param lambda two selections possible for optimal parameter of lambda. "bic" finds lambda from bic criteria, or user can directly input the penalty value.
#' @param nboot the number of bootstrap sample for pvalue. Default is 100.
#' @param n.cl number of cores in parallel computing. The default is (machine cores - 1)
#' @param bound bound of bic search in "bic" rule. Default is (.001, 1)
#' @param gridTF Utilize a grid search to optimize hyperparameters
#' @return \strong{pval} The empirical p-value for testing the equality of connectivity structure
#' @return \strong{rho} The sequence of penalty paramter based on the combined sample, subY1 and subY2.
#' @return \strong{fit0} Output of glasso for combied sample
#' @return \strong{fit1} Output of glasso for subY1
#' @return \strong{fit2} Output of glasso for subY2
#' @examples test1= testGlasso(testsim$X, testsim$Y, n.cl=1)
#' @export

testGlasso <- function(subY1, subY2, p, lambda = "bic", nboot = 100, n.cl, bound = c(.001, 1), gridTF=FALSE) {
  if (missing(n.cl)) {
    n.cl <- parallel::detectCores(logical = TRUE) - 1
  }
  if (missing(p)) {
    p <- .2
  }

  #library(glasso)
  #library(doParallel)

  ## Data dimension is Tt*q
  ### Auxilary functions
  glasso.bic2 <- function(y, khat, lambda, bound, debiasTF=TRUE) {
    st <- khat[1]
    end <- khat[2]
    id <- seq(st, end, by = 1)
    tj <- length(id)
    cy <- stats::cov(y[id, ])
    n <- nrow(y);

    
    if (lambda == "bic") {
      bic.rho <- function(rho, cy, tj) {
        B <- 0
        out <- glasso::glasso(cy, rho, maxit=400)
        ell1 <- 1 * (out$wi != 0)
        diag(ell1) <- rep(0, nrow(out$wi))
        ## tj is replacec by tj-1 by comparing result in the below.
        dwi <- sum(log(svd(out$wi)$d));
        B <- sum(diag((tj - 1) * cy %*% out$wi)) - tj * dwi + sum(ell1) * (log(tj))/2;
        return(B)
      }

      if(gridTF){
        ## Grid search?? Most stable..
        sqrho = seq(bound[1],  bound[2], length=20);
        est = lapply(sqrho, bic.rho, cy = cy, tj=tj);
        id = which.min(unlist(est)); rho = sqrho[id];
      } else{
        est <- stats::optim(par = .2, fn = bic.rho, cy = cy, tj = tj, method = "L-BFGS-B", lower = bound[1], upper = bound[2], control = list(maxit = 400, factr = 1e2, lmm = 4))
        rho <- est$par
      }

    } else {
      rho <- lambda
    } ## Fixed constant

    out <- glasso::glasso(cy, rho, maxit=400)
    idz <- which(out$wi == 0)

    ## If close to singular, do not debias
    if(det(out$w) < 1e-6){ debiasTF=FALSE; }

    ## Re-estimate precision matrix by imposing zero constraints
    ## It only valid when there is zero estimates in out()
    if ((length(idz) > 0)*debiasTF) {
      q <- dim(out$wi)[1]
      I <- matrix(rep(1:q, each = q), byrow = TRUE, ncol = q)
      J <- matrix(rep(1:q, each = q), ncol = q)
      zeroid <- cbind(I[idz], J[idz])
      out.zero <- glasso::glasso(cy, 0, zero = zeroid, trace = FALSE, maxit=10)
      ell1 <- 1 * (out.zero$wi != 0)
      diag(ell1) <- rep(0, nrow(out.zero$wi))
      det.value = sum(log(svd(out.zero$wi)$d));
      B = sum(diag((tj-1)*cy%*%out.zero$wi)) - tj*det.value + sum(ell1)*(log(n))/2
      out.zero$BIC <- B
    } else {
      ell1 <- 1 * (out$wi != 0)
      diag(ell1) <- rep(0, nrow(out$wi))
      dwi <- sum(log(svd(out$wi)$d));
      out$BIC <- sum(diag((tj - 1) * cy %*% out$wi)) - tj * dwi + sum(ell1) * (log(n))/2
      out.zero <- out
    }

    out.zero$debiasTF= debiasTF;
    out.zero$rho <- rho
    out.zero$wiold <- out$wi
    options(warn = 0)
    return(out.zero)
  }


  ### Speed up bootstrapping by using the same rho inside bootstrap sample

  boot.Cribben2 <- function(zz, tcp, p, rhos, n.cl, bound, nboot) {
    T <- dim(zz)[1]
    q <- dim(zz)[2]
    B <- nboot

    # building a stationary bootstrap sample
    if (missing(p)) {
      p <- .2
    }

    # similar but for Cribben et al

    #library(doParallel)
    cl <- parallel::makeCluster(n.cl)
    doParallel::registerDoParallel(cl)
    Bstat <- foreach::foreach(b = 1:B, .combine = c, .export = c("glasso.bic2"), .packages = c("glasso")) %dopar% {
      # building a stationary bootstrap sample
      size <- 0
      z_b <- NULL
      while (size < T) {
        b <- 1 + stats::rgeom(1, p)
        i <- sample(1:T, 1)
        if (i + b - 1 <= T) {
          z_b <- rbind(z_b, zz[i:(i + b - 1), ])
        } else {
          z_b <- rbind(z_b, zz[i:T, ], zz[(1:min(i - 1, i + b - T)), ])
        }
        size <- size + b
      }
      z_b <- z_b[(1:T), ]

      out <- glasso.bic2(z_b, khat = c(1, T), lambda = rhos[1], bound = bound)
      out1 <- glasso.bic2(z_b, khat = c(1, tcp), lambda = rhos[2], bound = bound)
      out2 <- glasso.bic2(z_b, khat = c(tcp + 1, T), lambda = rhos[3], bound = bound)

      dif <- out$BIC - out1$BIC - out2$BIC
    }
    parallel::stopCluster(cl)
    return(Bstat)
  }

  #################################
  ### Main function starts here
  # Combine sample
  #################################
  y <- rbind(subY1, subY2)
  n <- dim(y)[1]
  khat <- nrow(subY1)

  fit0 <- glasso.bic2(y, khat = c(1, n), lambda = lambda, bound = bound)
  fit1 <- glasso.bic2(y, khat = c(1, khat), lambda = lambda, bound = bound)
  fit2 <- glasso.bic2(y, khat = c(khat + 1, n), lambda = lambda, bound = bound)
  diff <- fit0$BIC - (fit1$BIC + fit2$BIC)
  rhos <- c(fit0$rho, fit1$rho, fit2$rho)
  bsample <- boot.Cribben2(y, tcp = khat, p, rhos = rhos, n.cl = n.cl, bound = bound, nboot = nboot)
  pval <- sum(bsample >= diff) / length(bsample)

  out <- list()
  out$bsample = bsample;
  out$fit0 <- fit0
  out$fit1 <- fit1
  out$fit2 <- fit2
  out$diff <- diff
  out$pval <- pval
  out$rho <- rhos

  return(out)
}

Try the detectR package in your browser

Any scripts or data that you put into this service are public.

detectR documentation built on Feb. 8, 2021, 5:06 p.m.