R/Thin_Statistic_Samples.R

In GERGM: Estimation and Fit Diagnostics for Generalized Exponential Random Graph Models

#' @title A function to thin out simulated GOF statistics so that they have an
#' MCMC autocorrelation of less tahn 0.01.
#' @description Takes as input the `@simulated_statistics_for_GOF` field of the
#' GERGM object (a data.frame) and returns the same data frame but now thinned
#' to reduce autocorrelation in the samples. Useful for exactly replicating the
#' statistics used in the GOF() function.
#'
#' @param statistics A data.frame stored in the
#' `@simulated_statistics_for_GOF` field of the GERGM object.
#' @return A data.frame that has been thinned to reduce autocorrelation.
#' @export
Thin_Statistic_Samples <- function(statistics){

  if(nrow(statistics) > 1 & length(unique(statistics$ttriads)) > 1){
    # we are going to test using ttriads statistics
    ttriads <- statistics$ttriads
    ar1 <- stats::cor(ttriads[2:length(ttriads)],ttriads[1:(length(ttriads)-1)])
    if (is.na(ar1) | is.nan(ar1) | is.null(ar1)) {
      cat("There was perfect autocorrelation, unable to thin statistics...\n")
    } else {
      #print(ar1)
      thin <- 1
      while (ar1 > .01) {
        thin = thin + 1
        thinSeq <- round(seq(1,length(ttriads),by=thin))
        thinDens <- ttriads[thinSeq]
        ar1 <- stats::cor(thinDens[2:length(thinDens)],thinDens[1:(length(thinDens)-1)])
        # break if thinning produces perfect autocorrelation
        if (is.na(ar1) | is.nan(ar1) | is.null(ar1)) {
          thin <- 1
          break
        }
      }

      thinSeq <- round(seq(1,nrow(statistics),by=thin))

      if (length(thinSeq) > 99) {
        cat("Thinning statistics to correct for autocorrelation in calculating fit diagnostics...\n")
        cat("Statistics were thinned by a factor of ",thin,", resulting in ",length(thinSeq)," samples.\n", sep = "")
        statistics <- statistics[thinSeq,]
      }else{
        cat("Thinning statistics to correct for autocorrelation in calculating fit diagnostics resulted in two few samples (less than 100), consider increasing the number of iterations of MCMC sampling...\n")
        cat("Statistics were thinned by a factor of ",thin,", resulting in ",length(thinSeq)," samples. Sample size was manually increased to 100.\n", sep = "")
        thinSeq <- round(seq(1,nrow(statistics),length.out = 100))
        statistics <- statistics[thinSeq,]
      }
    }
  }else{
    cat("Could not thin statistics, no variation...\n")
  }

  return(statistics)
}