R/gen_pcors.R

In donaldRwilliams/IRCcheck: Irrepresentable Condition Check

#' Generate True Partial Correlation Matrix
#'
#' @param p number of variables (nodes)
#'
#' @param edge_prob connectivity
#'
#' @param lb lower bound for the partial correlations
#'
#' @param ub upper bound for the partial correlations
#'
#' @note The function checks for a valid matrix (positive definite),
#' but sometimes this will
#' still fail. For example, for larger \code{p}, to have
#' large partial correlations this requires a sparse GGM
#' (accomplished by setting \code{edge_prob} to a small value).
#'
#' @return A list with the true structure, adjacency matrix, and correlation matrix.
#'
#' @export
#'
#' @importFrom stats runif
#'
#' @examples
#'
#' true_net <- gen_net(p = 10)
gen_net <- function(p = 20, edge_prob = 0.3,
                    lb = 0.05,
                    ub = 0.3) {
  d <- -1
  
  trys <- 0
  
  while (d < 0) {
    
    trys <- trys + 1
    
    effects <- p * (p - 1) * 0.5
    
    mat <- matrix(1, p, p)
    
    prob_zero <- 1 - edge_prob
    
    pool <- c(rep(0, effects * prob_zero),
              runif(effects * edge_prob, lb, ub))
    
    if (length(pool) != effects) {
      pool <- c(0, pool)
    }
    
    mat[upper.tri(mat)] <- sample(pool, size = effects)
    pcs <- symm_mat(mat)
    pcs <- -pcs
    diag(pcs) <- -diag(pcs)
    d <- det(pcs)
    
  }
  
  cors <- cov2cor(solve(pcs))
  
  inv <- solve(cors)
  
  pcors <- cov2cor(inv) * -1
  
  diag(pcors) <- 1
  
  adj <- ifelse(pcs == 0, 0, 1)
  
  returned_object <- list(
    pcors = pcors * adj,
    cors = cors,
    trys = trys,
    adj = adj
  )
  
  returned_object
  
}