R/bootstrap.R

In davidycliao/redguards: A replication package for the paper of "Factionalism and the Red Guards under Mao’s China: Ideal Point Estimation Using Text Data."

#' @export bootstrap_poisIRT
#' @rdname bootstrap_poisIRT
#' @param emIRT.out an emIRT() object, which is output from a call to poisIRT
#' @param .data the data used to produce the emIRT object.  
#' @param .starts a list containing several matrices of starting values for the parameters, which include alpha A (J x 1) matrix , psi A (K x 1) matrix, beta A (J x 1) matrix and x An (NI x 1) matrix.
#' @param .priors list, containing several matrices of starting values for the parameters. .priors can be generated by create_prior() by default via mu = 0 and variance = 100.  
#' @importFrom emIRT poisIRT
#' @importFrom utils flush.console
#' @importFrom utils tail
#' @importFrom stats sd
#' @title  Parametric Bootstrap of EM Standard Errors
#' @description bootstrap_poisIRT take an poisIRT() object (from emIRT package) and implements a parametric bootstrap of the standard errors for the ideal points. 
#' It assumes you have already run the Poisson IRT estimation via EM , and takes the output from that estimate, 
#' along with the original dataset, word frequency matrix. This function will conduct a bootstrap by running the estimates from sub-sampled observations based on the a poisson distribution.
bootstrap_poisIRT <- function(emIRT.out, .data, .starts, .priors,
                              .control = {list(threads = 1, verbose = TRUE, thresh = 1e-6, maxit = 5000)},
                              set.seed = 1234,
                              Ntrials = 5){
  if (class(.data)[1] == "wfm"){
    if (class(emIRT.out)[1] == "poisIRT") {
      wfm = .data
      # austin-wordfish: theta document positions
      # emIRT-poisIRT:x An (NI x 1) matrix of point estimates for the actor ideal points x  
      # lout$means$x
      idealpts = emIRT.out$means$x
      # austin-wordfish: psi word fixed effects
      # beta A (J x 1) matrix of point estimates for the word discrimination parameter β.
      # lout$means$beta
      beta.poisIRT = emIRT.out$means$beta    
      lambda = exp(outer(as.vector(emIRT.out$vars$beta), idealpts))    
      poisIRT.trials = matrix(NA, nrow = length(idealpts), ncol = Ntrials)
      for (trial in 1:Ntrials) {
        # task is to draw random integers from the Poisson distribution with given λ
        set.seed(set.seed)
        rpois_mat = matrix(rpois(rep(1, length(idealpts) * length(beta.poisIRT)), 
                                 lambda = lambda), nrow = length(beta.poisIRT))
        rownames(rpois_mat) = rownames(wfm)
        colnames(rpois_mat) = colnames(wfm) 
        sink("emIRTjunk.kjz")
        poisIRT.trials[, trial] = poisIRT(.rc =rpois_mat, 
                                          i = 0:(ncol(rpois_mat)-1),
                                           NI=ncol(rpois_mat),
                                           .starts = .starts,
                                           .priors = .priors, 
                                           .control = .control)$means$x      
      sink()
      unlink("emIRTjunk.kjz")
      if (isFALSE(trial%%Ntrials == 0)) { 
        cat("\n\t Bootstrap Iteration", trial, "complete...")
        utils::flush.console()}
      else if (isTRUE(trial%%Ntrials == 0)){
        cat("\n\t Last Bootstrap Iteration", trial, "complete...")
      }
    }
      emIRT.out$bse$xtrial = poisIRT.trials
      emIRT.out$bse$x = apply(poisIRT.trials, 1, stats::sd)
    }else
          {stop("This is not poisIRT() object, please check it")}
    }else {
      stop(".data is not word frequency matrix. Please check the data sctructure of .data")
      }
  return(emIRT.out)
  }