NSUM main analysis

Set up

Libraries

# get packages needed
fpackage.check <- function(packages) { # (c) Jochem Tolsma
  lapply(packages, FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  })
}
packages = c("haven", "NSUM","matrixStats", "parallel", 
             "MASS", "doParallel", "tictoc", "beepr", "dplyr", "purrr",
             "readr")
fpackage.check(packages)

Import the NSUM prep files.

#import nsum input
load(file = "data_analysis/data/data_processed/nsum_input/analysis_files_selection.rds")

#extract mat and known vector
mat <- analysis_files_selection[[1]]
known_vector <- analysis_files_selection[[2]]

NSUM analysis

We will initialize the NSUM estimation. First we will set the seed to make it repricable. Second, we set the number of iterations on 40000, the burin at 1000, the number of retains results at 4000, the population size, and the number of models we need to run (holdouts)

#set seed
set.seed(220914)

#------------------------------------------------------------------------------------------------
# some info to work with
iters <- 40000# number of iterations (40k)
burns <- 1000 # burnin size  (1k)
retain <- 4000 # how many chains do we want to retain? (4k)
popsize <- 17407585 # population size
holdouts <- 1:length(known_vector)  

#------------------------------------------------------------------------------------------------

We will use parallel computing to speed up the NSUM calculations. My work lapop has 8 cores so we can use 7 cores for overnight computation. There are 16 models to estimate so we need 3 runs.

# paralellize the estimation
numCores <- detectCores()
if (numCores > length(holdouts-1)) {
  registerDoParallel(cores = length(holdouts))
} else {    
  print("Too few cores to do in one run!") 
} 

[1] “Too few cores to do in one run!”

# so we need to do in seperate runs
registerDoParallel(core=numCores-1)

Main analysis loop. We use the killworth start function to create starting values. We then use these in the analysis loop. We calculate for each respondent the average and sd over the 4000 retained chains. These are stored in separate datafiles.

#------------------------------------------------------------------------------------------------
# fill up empty lists
kds <- list()
kdssd <- list()
data_list <- list()

tic()
foreach(i = holdouts[1:16]) %dopar% { # change the number here on own max core you want to use #i = 1
  # 16 models to estimate. I use 7 cores, each core solves one per night.
  # so 16/7 = 2.3 = 3 nights.
  
  cd <- mat[, -c(i)] # take out pop of interest and make a new mat
  
  # calculate starting values
  z <- NSUM::killworth.start(cbind(cd, mat[, c(i)]), # paste the "takenout" at the END of matrix
                             known_vector[-c(i)], # this is the known pop, WITHOUT unknown pop
                             popsize) # population size
  
    file_name_data <- paste0("data_analysis/results/nsum_output/main/data/estimates_holdout", i, ".txt") #create file.name for data
    file_name_model <- paste0("data_analysis/results/nsum_output/main/model/estimates_holdout", i, ".rds") #create file.name for model
    if(!file.exists(file_name_model)){#check for file, if it does not exist run the NSUM model
      # run the Bayesian estimation
      degree <- NSUM::nsum.mcmc(cbind(cd, mat[, c(i)]), # gets pasted at the last column
                                known_vector[-c(i)], # here we again take out the "holdout", or artificial "unknown" pop
                                popsize, model = "combined", # combined control for transmission and recall errors
                                indices.k = length(holdouts), # notice that "holdout" gets pasted as last column
                                iterations = iters, burnin = burns, size = retain, # 40k iterations, retain 4k chains
                                d.start = z$d.start, mu.start = z$mu.start, # starting values from simple estimator
                                sigma.start = z$sigma.start, NK.start = z$NK.start)
      
      #store model_results
      save(degree, file = file_name_model)
      
    } else {
      load(file_name_model)
    }

    #store mean and sd in df
    # calculate rowmean and it's SD of the retained 4k chains
    kds[[i]] <- rowMeans(degree$d.values,na.rm = TRUE) # calculate rowmean of netsize iterations: so the retained chains
    kdssd[[i]] <- matrixStats::rowSds(degree$d.values) # calculate sd of 4k estimates per row: sd for those values
    data_list[[i]] <- data.frame(cbind(kds[[i]],kdssd[[i]])) # combine and put in df
    
    # Save the data, new .txt for each iteration, if something goes wrong, we can always start at prior one and combine
    write.table(data_list[[i]], file = file_name_data, row.names = F) #store them in results.
}

toc() # time passed?
beep() # done!

---

Copyright © 2024 Jeroense Thijmen