## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(collapse = TRUE, comment = "#>") # FALSE: build the document without running the calculations. run_examples <- FALSE # TRUE: use reduced calculations to check that the code runs. # FALSE: use the full computational settings. quick_check <- TRUE nboot <- if (quick_check) 10L else 1000L ntop_finite <- if (quick_check) 2L else 60L ntop_unrestricted <- if (quick_check) 2L else Inf ## ----reproduction-setup, eval=run_examples------------------------------------ # library(MultipleSystemsEstimation) # # vary_ntop_bca <- getFromNamespace("vary_ntop_bca", "MultipleSystemsEstimation") # downhill_fit <- getFromNamespace("downhill_fit", "MultipleSystemsEstimation") # downhill_bootstrapcal <- getFromNamespace("downhill_bootstrapcal", "MultipleSystemsEstimation") # downhill_jackknifecal <- getFromNamespace("downhill_jackknifecal", "MultipleSystemsEstimation") # assemble_bic <- getFromNamespace("assemble_bic", "MultipleSystemsEstimation") # bootstrapcal <- getFromNamespace("bootstrapcal", "MultipleSystemsEstimation") # degree2_order <- getFromNamespace(".degree2_order", "MultipleSystemsEstimation") # bcaconfvalues <- getFromNamespace("bcaconfvalues", "MultipleSystemsEstimation") # # # Load all datasets used in the calculations. # data(Kosovo) # data(Korea) # data(UKdat_5) # data(NewOrl_5) # # # Check that the corrected New Orleans data are being used. # if (sum(NewOrl_5[, ncol(NewOrl_5)]) != 185) { # stop("This vignette requires the corrected 185-case New Orleans data.") # } # # # This check applies to the complete candidate space, before any ntop restriction. # checked_assemble <- function(zdat, maxorder, expected_models = NULL) { # z <- assemble_bic(zdat, maxorder = maxorder) # if (!isTRUE(z$enumeration_maxorder == z$requested_maxorder)) { # stop("Enumeration order was reduced: requested ", z$requested_maxorder, # "; actual ", z$enumeration_maxorder, # ". Increase maxmodel before treating this search as exhaustive.") # } # if (!nrow(z$res)) stop("No candidate model has a valid original-data fit.") # if (!is.null(expected_models) && nrow(z$res) != expected_models) { # stop("Expected ", expected_models, " valid models, but found ", nrow(z$res), # ". Check the package version and data before comparing with the paper.") # } # z # } # # # Check the requested searches and the published Kosovo and Korea model counts. # model_spaces <- list( # Kosovo = list(data = Kosovo, order = ncol(Kosovo) - 2, expected = 113L), # Korea = list(data = Korea, order = 2, expected = 6L), # UK2 = list(data = UKdat_5, order = 2, expected = NULL), # UK4 = list(data = UKdat_5, order = 4, expected = NULL), # NewOrleans2 = list(data = NewOrl_5, order = 2, expected = NULL), # NewOrleans4 = list(data = NewOrl_5, order = 4, expected = NULL) # ) # for (spec in model_spaces) { # invisible(checked_assemble(spec$data, spec$order, spec$expected)) # } # rm(model_spaces, spec) ## ----table2-calculation, eval=run_examples------------------------------------ # data(Kosovo) # # kos <- vary_ntop_bca( # Kosovo, # maxorder = ncol(Kosovo) - 2, # ntopmax = ntop_unrestricted, # degree = 1, # nboot = nboot, # iseed = 1234, # alpha = c(0.025, 0.1, 0.9, 0.975) # ) ## ----table2-extract, eval=run_examples---------------------------------------- # kos$inference[kos$inference$ntop %in% c(1, 5, 10, max(kos$inference$ntop)), ] ## ----table3-calculation, eval=run_examples------------------------------------ # data(Korea) # # kor <- vary_ntop_bca( # Korea, # maxorder = 2, # ntopmax = ntop_unrestricted, # degree = 1, # nboot = nboot, # iseed = 1234, # alpha = c(0.025, 0.1, 0.9, 0.975) # ) # ## ----table3-extract, eval=run_examples---------------------------------------- # kor$inference[ # kor$inference$ntop %in% c(1, 2, max(kor$inference$ntop)), # , # drop = FALSE # ] ## ----fig1-calc, eval=run_examples--------------------------------------------- # data(UKdat_5) # # uk2 <- vary_ntop_bca( # UKdat_5, # maxorder = 2, # ntopmax = ntop_finite, # degree = 1, # nboot = nboot, # iseed = 1234, # alpha = c(0.025, 0.1, 0.9, 0.975) # ) # # uk4 <- vary_ntop_bca( # UKdat_5, # maxorder = 4, # ntopmax = ntop_finite, # degree = 1, # nboot = nboot, # iseed = 1234, # alpha = c(0.025, 0.1, 0.9, 0.975) # ) ## ----fig2-calc, eval=run_examples--------------------------------------------- # data(NewOrl_5) # # no2 <- vary_ntop_bca( # NewOrl_5, # maxorder = 2, # ntopmax = ntop_finite, # degree = 1, # nboot = nboot, # iseed = 1234, # alpha = c(0.025, 0.1, 0.9, 0.975) # ) # # no4 <- vary_ntop_bca( # NewOrl_5, # maxorder = 4, # ntopmax = ntop_finite, # degree = 1, # nboot = nboot, # iseed = 1234, # alpha = c(0.025, 0.1, 0.9, 0.975) # ) ## ----table4-calc, eval=run_examples------------------------------------------- # bootstrap_optimal_ranks <- function(zdat, maxorder, nboot = 1000, iseed = 1234) { # # Exercise the exhaustive-search code on a small dataset in quick mode. # if (quick_check) { # zdat <- Korea # maxorder <- 2 # } # z <- checked_assemble( # zdat, maxorder = maxorder # ) # # z <- bootstrapcal( # z, nboot = nboot, iseed = iseed # ) # # original_order <- order(z$res[, "BIC"]) # # vapply(seq_len(nboot), function(b) { # v <- z$bootbic[, b] # ok <- which(is.finite(v)) # if (!length(ok)) return(NA_integer_) # # best <- ok[which.min(v[ok])] # match(best, original_order) # }, integer(1)) # } # # rank_counts <- function(r, ntop = c(1, 5, 10, 50, 100)) { # vapply( # ntop, # function(k) sum(r <= k, na.rm = TRUE), # integer(1) # ) # } # # ranks4 <- list( # Kosovo = bootstrap_optimal_ranks( # Kosovo, maxorder = ncol(Kosovo) - 2, nboot = nboot # ), # UK2 = bootstrap_optimal_ranks(UKdat_5, maxorder = 2, nboot = nboot), # UK4 = bootstrap_optimal_ranks(UKdat_5, maxorder = 4, nboot = nboot), # NewOrleans2 = bootstrap_optimal_ranks(NewOrl_5, maxorder = 2, nboot = nboot), # NewOrleans4 = bootstrap_optimal_ranks(NewOrl_5, maxorder = 4, nboot = nboot) # ) # # sapply(ranks4, rank_counts) # # # Counts above exclude failed replications; disclose their denominator. # rank_diagnostics4 <- t(vapply(ranks4, function(r) { # c(attempted = length(r), usable = sum(is.finite(r)), # failed = sum(!is.finite(r))) # }, numeric(3))) # rank_diagnostics4 ## ----table5-kosovo, eval=run_examples----------------------------------------- # kos_d2 <- vary_ntop_bca( # Kosovo, # maxorder = ncol(Kosovo) - 2, # ntopmax = ntop_unrestricted, # degree = 2, # nboot = nboot, # iseed = 1234, # alpha = c(0.025, 0.975) # ) # # kos$inference[kos$inference$ntop %in% c(1, 5, 10, max(kos$inference$ntop)), ] # kos_d2$inference[kos_d2$inference$ntop %in% c(1, 5, 10, max(kos_d2$inference$ntop)), ] ## ----fig3-reproduction-code, eval=run_examples-------------------------------- # uk_d1 <- vary_ntop_bca( # UKdat_5, maxorder = 4, ntopmax = ntop_unrestricted, # degree = 1, nboot = nboot, iseed = 1234, # alpha = c(0.025, 0.975) # ) # # uk_d2 <- vary_ntop_bca( # UKdat_5, maxorder = 4, ntopmax = ntop_unrestricted, # degree = 2, nboot = nboot, iseed = 1234, # alpha = c(0.025, 0.975) # ) # # no_d1 <- vary_ntop_bca( # NewOrl_5, maxorder = 4, ntopmax = ntop_unrestricted, # degree = 1, nboot = nboot, iseed = 1234, # alpha = c(0.025, 0.975) # ) # # no_d2 <- vary_ntop_bca( # NewOrl_5, maxorder = 4, ntopmax = ntop_unrestricted, # degree = 2, nboot = nboot, iseed = 1234, # alpha = c(0.025, 0.975) # ) ## ----fig4-calc, eval=run_examples--------------------------------------------- # bootstrap_rank_pair <- function(zdat, maxorder, # nboot = 1000, iseed = 1234) { # # Exercise the exhaustive-search code on a small dataset in quick mode. # if (quick_check) { # zdat <- Korea # maxorder <- 2 # } # z <- checked_assemble( # zdat, maxorder = maxorder # ) # # z <- bootstrapcal( # z, nboot = nboot, iseed = iseed # ) # # r1_order <- order(z$res[, "BIC"]) # # r2_order <- degree2_order( # rownames(z$res), # nlists = ncol(z$xdata) - 1, # maxorder = z$maxorder # ) # # r1 <- rep(NA_integer_, nboot) # r2 <- rep(NA_integer_, nboot) # # for (b in seq_len(nboot)) { # v <- z$bootbic[, b] # ok <- which(is.finite(v)) # if (!length(ok)) next # # best <- ok[which.min(v[ok])] # r1[b] <- match(best, r1_order) # r2[b] <- match(best, r2_order) # } # # usable <- is.finite(r1) & is.finite(r2) # keep <- usable & r1 > 1 # # pairs <- data.frame( # ordinary = r1[keep], # degree2 = r2[keep] # ) # attr(pairs, "replication_counts") <- c( # attempted = nboot, usable = sum(usable), failed = sum(!usable), # rank_one_omitted = sum(usable & r1 == 1), retained_pairs = sum(keep) # ) # pairs # } # # kos_rank <- bootstrap_rank_pair( # Kosovo, # maxorder = ncol(Kosovo) - 2, nboot = nboot # ) # # no2_rank <- bootstrap_rank_pair(NewOrl_5, maxorder = 2, nboot = nboot) # no4_rank <- bootstrap_rank_pair(NewOrl_5, maxorder = 4, nboot = nboot) # uk2_rank <- bootstrap_rank_pair(UKdat_5, maxorder = 2, nboot = nboot) # uk4_rank <- bootstrap_rank_pair(UKdat_5, maxorder = 4, nboot = nboot) # # rank_sets <- list( # Kos = kos_rank, # NO2 = no2_rank, # NO4 = no4_rank, # UK2 = uk2_rank, # UK4 = uk4_rank # ) # # # Numerical rank pairs underlying Figure 4. # rank_sets # rank_diagnostics_fig4 <- t(vapply(rank_sets, function(x) { # attr(x, "replication_counts") # }, numeric(5))) # rank_diagnostics_fig4 ## ----table6-korea, eval=run_examples------------------------------------------ # downhill_korea <- list( # estimate = downhill_fit( # counts = Korea[, ncol(Korea)], # desmat = Korea[, -ncol(Korea), drop = FALSE], # maxorder = 2 # ), # bootreps = downhill_bootstrapcal( # Korea, # nboot = nboot, # iseed = 1234, # maxorder = 2 # ), # ahat = downhill_jackknifecal( # Korea, # maxorder = 2 # ) # ) # # c( # estimate = downhill_korea$estimate, # acceleration = downhill_korea$ahat # ) # # usable <- is.finite(downhill_korea$bootreps) # downhill_korea$replication_counts <- c( # attempted = length(usable), usable = sum(usable), failed = sum(!usable) # ) # downhill_korea$replication_counts # # downhill_korea$confidence_limits <- NULL # if (quick_check) { # message("Korea BCa limits omitted: ten replications are insufficient for meaningful tail estimates.") # } else if (!all(usable) || !is.finite(downhill_korea$estimate) || # !is.finite(downhill_korea$ahat)) { # warning("Korea BCa limits omitted: investigate failed fits or unavailable acceleration.") # } else { # # Extreme bias correction can make BCa probabilities undefined. # p_below <- mean(downhill_korea$bootreps < downhill_korea$estimate) # if (p_below == 0 || p_below == 1) { # warning("Korea BCa limits omitted: all bootstrap estimates lie on one side of the point estimate.") # } else { # downhill_korea$confidence_limits <- bcaconfvalues( # bootreps = downhill_korea$bootreps, # popest = downhill_korea$estimate, # ahat = downhill_korea$ahat, # alpha = c(0.025, 0.975) # ) # downhill_korea$confidence_limits # } # }