R/random_biLCM_data.R

In insongkim/polnet: Mapping Political Communities: A Statistical Analysis of Lobbying Networks in Legislative Politics

#'@param m The number of group1
#'@param n The number of group2
#'@param k The number of link communities
#'@param kappa_weight The base weight. Expecting a k length vector. Assumed to be drawn from a gamma distribution with a shape of a and rate of b. An optional argument, but must specifiy this or a and b.
#'@param a The shape parameter to generate kappa. Expecting a positive real number. An optional argument, required if missing kappa_weight.
#'@param b The rate parameter to generate kappa. Expecting a positive real number. An optional argument, required if missing kappa_weight.
#'@param alpha_membership The desired latent mixed-membership of the group1. Expecting a matrix with m rows and k columns. Assumed to be drawn from a dirichlet distribution with concentration parameters alpha_c. An optional argument, but must specify this or alpha_c.
#'@param alpha_c The concentration vector to generate alpha_membership. Expecting a vector of length m. An optional argument, required if missing alpha_membership.
#'@param beta_membership The desired latent mixed-membership of the group2. Expecting a matrix with n rows and k columns. Assumed to be drawn from a dirichlet distribution with concentration parameters beta_c. An optional argument, but must specify this or beta_c.
#'@param beta_c The concentration vector to generate beta_membership. Expecting a vector of length n. An optional argument, required if missing beta_membership.
#'@param non_zero If non_zero = TRUE, repeat the data generating process untill each member has at least one interaction with others.
#'@return A list of kappa_weight, alpha_membership, beta_membership, and then the randomly generated poisson count matrix A.

#'@export random_biLCM_data

random_biLCM_data <- function(m, 
                              n, 
                              k,
                              kappa_weight = NULL,
                              a = NULL,
                              b = NULL,
                              alpha_membership = NULL,
                              alpha_c = NULL,
                              beta_membership = NULL,
                              beta_c = NULL,
                              non_zero = FALSE) {
  
  repeat {
    # Generating kappa
    if(is.null(kappa_weight)){
      if(is.null(a) | is.null(b)){
        stop("Invalid Kappa Parameters.\n")
      }else{
        kappa_weight <- rgamma(k, shape = a, rate = b)
      }
    }
    
    # Generating alpha
    if(is.null(alpha_membership)){
      if(is.null(alpha_c)){
        stop("Invalid Alpha Parameters.\n")
      }else{
        alpha_membership  <- t(rdirichlet(k, alpha_c))
      }
    }
    
    # Generating beta
    if(is.null(beta_membership)){
      if(is.null(beta_c)){
        stop("Invalid Alpha Parameters.\n")
      }else{
        beta_membership  <- t(rdirichlet(k, beta_c))
      }
    }
    
    mu <- sweep(alpha_membership, 2, kappa_weight, "*") %*% t(beta_membership)
    A_mat <- matrix(sapply(mu, function(x) rpois(1, x)), m, n)
    
    if (((!0 %in% colSums(A_mat)) & (!0 %in% rowSums(A_mat))) | !non_zero) break
  }
  
  
  
  return(list(kappa = kappa_weight, alpha = alpha_membership, beta = beta_membership, A = A_mat))
}


rdirichlet = function(n, alpha) {
  l <- length(alpha)
  sample <- matrix(rgamma(l*n, alpha), ncol = l, byrow = TRUE)
  sample <- sample/rowSums(sample)
  
  return(sample)
}

#'@param biLCM_Object A trained object of class biLCM
#'@param simulated_data A list of synthetic data generated by the function random_biLCM_data
#'@param group1 If group1 = TRUE, show the result of a member of group1 (i.e., alpha_iz's). Otherwise that of group2 (i.e., beta_jz's)
#'@param nth Show the result of 'n'th member's community distribution
#'@return Two pie plots: true and estimated community distribution (either alpha_iz's or beta_jz's)

#'@export plot.compare.biLCM

plot.compare.biLCM <- function(biLCM_Object,
                               simulated_data,
                               group1 = TRUE,
                               nth) {
  if(class(biLCM_Object)!="biLCM") stop("'biLCM_object' is not of class 'biLCM'.\n")
  
  par(mfrow=c(1,2))
  
  k1 <- length(simulated_data$kappa)
  cols1 <- gg_color_hue(k1, alpha = 0.7)
  
  k2 <- length(biLCM_Object$kappa)
  cols2 <- gg_color_hue(k2, alpha = 0.7)
  
  if (group1) {
    pie(sort(simulated_data$alpha[nth,]), order(simulated_data$alpha[nth,]), col = cols1, clockwise=TRUE,
        main = paste0("True\nCommunity Distribution (i=",nth,")"))
    pie(sort(biLCM_Object$alpha[nth,]), order(biLCM_Object$alpha[nth,]), col = cols2, clockwise=TRUE,
        main = paste0("Estimated\nCommunity Distribution (i=",nth,")"))
  } else {
    pie(sort(simulated_data$beta[nth,]), order(simulated_data$beta[nth,]), col = cols1, clockwise=TRUE,
        main = paste0("True\nCommunity Distribution (j=",nth,")"))
    pie(sort(biLCM_Object$beta[nth,]), order(biLCM_Object$beta[nth,]), col = cols2, clockwise=TRUE,
        main = paste0("Estimated\nCommunity Distribution (j=",nth,")"))
  }
  
  par(mfrow=c(1,1))
  
}