R/bootERGM.R

In schmid86/bootERGM: Exponential Random Graph Models with Parametric Bootstrapped Pseudolikelihood

####

# include: checkdegeneracy function, gof

# dependencies: ergm, doParallel, foreach
#library(ergm)
#library(doParallel)
#library(foreach)


bootERGM <- function(formula, R= 500, control.mple = control.ergm(drop=TRUE),
                     control.sim = control.simulate.ergm(MCMC.interval=1024),
                     control.boot=control.ergm(drop=TRUE), nowonit=100, cpu.cores=1, constraints =~.){

  #control.boot <- control.mple
   # fit model using MPLE
   mple <- ergm(formula, estimate="MPLE" , control=control.mple, constraints = constraints)

   observed.statistics <- summary(formula)

   # get number of variables in model
   num.variables <- length(coef(mple))

   #simulate network from mple
   cat("Simulating ", R , " networks\n")
   sim.mple <- simulate(mple, nsim=R, control=control.sim ) # simulate network

   # create empty matrix to store mple of bootstrap samples
   boot.mple.mat <- matrix(0, R, num.variables)
   colnames(boot.mple.mat) <- names(coef(mple))

   # create empty matrix to store network statistics of bootstrap samples
   boot.stat.mat <- matrix(0, R, num.variables)
   colnames(boot.stat.mat) <- names(coef(mple))

   # replace response network in formula
   tt <- terms(formula)
   new.formula <- reformulate(attr(tt, "term.labels"), "sim.mple[[i]]")

   cat("Estimating MPLE of simulated networks\n")

   # if only one CPU core is used

   if(cpu.cores==1){
   for(i in 1:R){

     if(i %% nowonit == 0)
     {
       cat("Now on iteration ", i, "\n")
     } # end if

     boot.mple<- suppressMessages(ergm(new.formula, estimate="MPLE" , control=control.boot, constraints=constraints))

     boot.mple.mat[i,] <- coef(boot.mple)

     boot.stat.mat[i,] <- summary(new.formula)
   } # end for


   } else{ # parallel

    registerDoParallel(cpu.cores=cpu.cores)

    par.results <-foreach(i=1:R, .export='ergm', .combine=rbind)%dopar%{



    tt <- terms(formula)
    sim.mple.net <- sim.mple[[i]]
    par.formula <- reformulate(attr(tt, "term.labels"), "sim.mple.net")

    boot.mple<- suppressMessages(ergm(par.formula, estimate="MPLE" , control=control.boot))

    g1 <- coef(boot.mple)

    g2 <- summary(par.formula)

    g<- c(g1, g2)

    g

    } # end foreach

    # split par.results into boot.mple.mat and booot.stat.mat

    boot.mple.mat <- par.results[,1:num.variables ]
    boot.stat.mat <- par.results[,(num.variables+1):(2*num.variables) ]

    } # end else # parallel

    # create empty matrix
    Result.matrix <- matrix(0, num.variables, 3)
    rownames(Result.matrix) <- names(coef(mple))
    colnames(Result.matrix) <- c("MPLE", "2.5%", "97.5%")

    Result.matrix[,1] <- coef(mple)

    for(i in 1:num.variables){
      Result.matrix[i,2] <- quantile(boot.mple.mat[,i], probs = 0.025)
      Result.matrix[i,3] <- quantile(boot.mple.mat[,i], probs = 0.975)

    }

   return.list <- list(boot.interval = Result.matrix , boot.mple.matrix=boot.mple.mat , MPLE= mple, formula=formula,
                       observed.statistics=observed.statistics, bootstrap.statistics=boot.stat.mat)
   class(return.list) <- "bootERGM"

   return(return.list)
}