R/hergm.postprocess.R

In hergm: Hierarchical Exponential-Family Random Graph Models

###########################################################################
# Copyright 2009 Michael Schweinberger                                    #
#                                                                         #
# This file is part of hergm.                                             #
#                                                                         # 
#    hergm is free software: you can redistribute it and/or modify        #
#    it under the terms of the GNU General Public License as published by #
#    the Free Software Foundation, either version 3 of the License, or    #
#    (at your option) any later version.                                  #
#                                                                         # 
#    hergm is distributed in the hope that it will be useful,             #
#    but WITHOUT ANY WARRANTY; without even the implied warranty of       #
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        #
#    GNU General Public License for more details.                         #
#                                                                         #
#    You should have received a copy of the GNU General Public License    #
#    along with hergm.  If not, see <http://www.gnu.org/licenses/>.       #
#                                                                         # 
###########################################################################

hergm.postprocess <- function(object,
                              burnin = 2000, 
                              thinning = 1,
                              relabel = 1,
                              number_runs = 1,
                              ...)
# input: MCMC sample generated by function hergm
# output: postprocessed MCMC sample
{
  object.hergm <- object

  # Extract 
  method <- object$method
  n <- object$n
  max_number <- object$max_number
  number_fixed <- object$number_fixed
  d1 <- object$d1
  d2 <- object$d2
  parallel <- object$parallel
  simulate <- object$simulate
  sample_size <- object$sample_size
  mcmc <- object$sample 
  predictions <- object$predictions
  verbose <- object$verbose

  # Check
  if ((d2 == 0) || (max_number == 1) || (number_fixed == n)) relabel <- 0 
  if (((sample_size - burnin) / thinning) < 1000)
    {
    burnin <- 0
    thinning <- 1
    }
 
  if (method == "bayes")
  {
  # Preprocess MCMC sample: delete burn-in iterations and transform vector into matrix, where rows correspond to MCMC draws
  d <- d1 + d2
  terms <- length_mcmc(d1, d2, max_number, n, predictions)  
  mcmc_sample <- NULL
  count <- 0
  for (i in 1:parallel) 
    {
    first <- count + (burnin * terms) + 1
    last <- count + (sample_size * terms)
    mcmc_sample <- append(mcmc_sample, mcmc[first:last])
    count <- count + (sample_size * terms)
    }
  mcmc_sample_size <- (parallel * (sample_size - burnin)) 

  mcmc_sample <- matrix(mcmc_sample, nrow = mcmc_sample_size, ncol = terms, byrow = TRUE)
  if (thinning > 1)
    {
    for (i in mcmc_sample_size:1) 
      {
      if (trunc(i / thinning) != (i / thinning)) mcmc_sample <- mcmc_sample[-i,] 
      }
    mcmc_sample_size <- nrow(mcmc_sample)
    }

  if (is.null(object$extract)) extract <- TRUE
  else extract <- object$extract

  if (extract == TRUE) # Extract MCMC sample from argument sample
  { # Begin: extract

  # Initialize arrays
  if (d1 > 0) 
    {
    object.hergm$ergm_theta <- matrix(data = 0, nrow = mcmc_sample_size, ncol = d1) 
    }
  if (d2 > 0)
    {
    object.hergm$eta_mean <- matrix(data = 0, nrow = mcmc_sample_size, ncol = d2)
    object.hergm$eta_precision <- matrix(data = 0, nrow = mcmc_sample_size, ncol = d2)
    object.hergm$hergm_theta <- matrix(data = 0, nrow = mcmc_sample_size, ncol = d2 * (max_number + 1))
    object.hergm$indicator <- matrix(data = 0, nrow = mcmc_sample_size, ncol = n)
    object.hergm$size <- matrix(data = 0, nrow = mcmc_sample_size, ncol = max_number)
    object.hergm$p_k <- matrix(data = 0, nrow = mcmc_sample_size, ncol = max_number)
    object.hergm$alpha <- matrix(data = 0, nrow = mcmc_sample_size, ncol = 1)
    }
  else 
    {
    object.hergm$indicator <- NULL
    object.hergm$hergm_theta <- NULL
    }
  if (predictions == TRUE) object.hergm$prediction <- matrix(data = 0, nrow = mcmc_sample_size, ncol = d)

  # Process MCMC sample
  for (row in 1:mcmc_sample_size)
    {
    column <- 0
    if (d1 > 0)
      {
      for (i in 1:d1) 
        {
        column <- column + 1
        object.hergm$ergm_theta[row,i] <- mcmc_sample[row,column]
        }
      }
    if (d2 > 0)
      {
      for (i in 1:d2) 
        {
        column <- column + 1
        object.hergm$eta_mean[row,i] <- mcmc_sample[row,column]
        }
      for (i in 1:d2) 
        {
        column <- column + 1
        object.hergm$eta_precision[row,i] <- mcmc_sample[row,column]
        }
      for (i in 1:(d2 * (max_number + 1))) 
        {
        column <- column + 1
        object.hergm$hergm_theta[row,i] <- mcmc_sample[row,column]
        }
      for (i in 1:n) 
        {
        column <- column + 1
        object.hergm$indicator[row,i] <- mcmc_sample[row,column]
        }
      for (i in 1:max_number) 
        {
        column <- column + 1
        object.hergm$size[row,i] <- mcmc_sample[row,column]
        }
      for (i in 1:max_number) 
        {
        column <- column + 1
        object.hergm$p_k[row,i] <- mcmc_sample[row,column]
        }
      column <- column + 1
      object.hergm$alpha[row,1] <- mcmc_sample[row,column]
      } 
    if (predictions == TRUE)
      {
      for (i in 1:d) 
        {
        column <- column + 1 
        object.hergm$prediction[row,i] <- mcmc_sample[row,column]
        }
      }
    }

  } # End: extract

  # Relabel sample
  if (relabel %in% c(1, 2)) # Minimizing posterior expected loss
    {
    if (relabel == 1) minimizer <- hergm.relabel_1(max_number=max_number, indicator=object.hergm$indicator, number_runs=number_runs, verbose=verbose) # Minimize posterior expected loss of Schweinberger and Handcock (2015)
    else minimizer <- hergm.relabel_2(max_number=max_number, indicator=object.hergm$indicator, verbose=verbose) # Minimize posterior expected loss of Peng and Carvalho (2015) 
    object.hergm$relabeled.indicator <- minimizer$indicator
    object.hergm$p_i_k <- minimizer$p
    object.hergm$relabeled.hergm_theta <- matrix(0, nrow=nrow(object.hergm$hergm_theta), ncol=ncol(object.hergm$hergm_theta))
    if (simulate == FALSE)
      {
      object.hergm$relabeled.hergm_theta <- matrix(0, nrow=nrow(object.hergm$hergm_theta), ncol=ncol(object.hergm$hergm_theta))
      index1 <- 1
      index2 <- max_number
      theta <- object.hergm$hergm_theta[,index1:index2]
      object.hergm$relabeled.hergm_theta[,index1:index2] <- hergm.permute_mcmc(theta, max_number, minimizer$min_permutations) # Within-block parameters 
      index2 <- index2 + 1 # Between-block parameters
      object.hergm$relabeled.hergm_theta[,index2] <- object.hergm$hergm_theta[,index2] # Between-block parameters
      if (d2 > 1)
        {
        for (h_term in 2:d2) # Relabel block-dependent parameters of block-dependent model terms one by one
          {
          index1 <- index2 + 1 # Increment starting index
          index2 <- index2 + max_number # Increment stopping index
          theta <- object.hergm$hergm_theta[,index1:index2] # Within-block parameters
          object.hergm$relabeled.hergm_theta[,index1:index2] <- hergm.permute_mcmc(theta, max_number, minimizer$min_permutations) # Within-block parameters 
          index2 <- index2 + 1 # Between-block parameters
          object.hergm$relabeled.hergm_theta[,index2] <- object.hergm$hergm_theta[,index2] # Copy between-block parameters
          }
        }
      if (verbose >= 0) cat("\n")
      }
    }
 
  # Store
  object.hergm$extract <- FALSE # If hergm.postprocess() is called hereafter, we will not extract the MCMC sample again, because otherwise we will get runtime errors
  object.hergm$n <- object$n
  object.hergm$network <- object$network
  object.hergm$model <- object$model 
  object.hergm$max_number <- object$max_number
  object.hergm$number_fixed <- object$number_fixed  
  object.hergm$d1 <- object$d1
  object.hergm$d2 <- object$d2
  object.hergm$hyper_prior <- object$hyper_prior
  object.hergm$parallel <- object$parallel
  object.hergm$simulate <- object$simulate
  object.hergm$sample_size <- mcmc_sample_size
  if (simulate == TRUE)
    {
    object.hergm$sample <- object$sample
    object.hergm$heads <- object$heads
    object.hergm$tails <- object$tails
    }
  object.hergm$relabel <- relabel
  object.hergm$predictions <- object$predictions
  object.hergm$verbose <- object$verbose
  if ((relabel == 0) && (!is.null(object$relabeled.hergm_theta))) 
    {
    object.hergm$relabel <- object$relabel # Restore original relabel choice
    object.hergm$relabeled.hergm_theta <- object$relabeled.hergm_theta # In other words, if we have been here before and have relabeled the MCMC sample, then store sample
    }
  }
  else
  {
  # Store
  object.hergm$extract <- FALSE # If hergm.postprocess() is called hereafter, we will not extract the MCMC sample again, because otherwise we will get runtime errors
  object.hergm$n <- object$n
  object.hergm$network <- object$network
  object.hergm$model <- object$model
  object.hergm$max_number <- object$max_number
  object.hergm$number_fixed <- object$number_fixed
  object.hergm$d1 <- object$d1
  object.hergm$d2 <- object$d2
  object.hergm$hyper_prior <- object$hyper_prior
  object.hergm$parallel <- object$parallel
  object.hergm$simulate <- object$simulate
  object.hergm$relabel <- relabel
  object.hergm$verbose <- object$verbose
  }

  object.hergm
}