Test sim kwan em edits

DESCRIPTION

@@ -1,12 +1,13 @@
 Type: Package
 Package: serocalculator
 Title: Estimating Infection Rates from Serological Data
-Version: 1.3.0.9047
+Version: 1.3.0.9048
 Authors@R: c(
     person("Peter", "Teunis", , "p.teunis@emory.edu", role = c("aut", "cph"),
            comment = "Author of the method and original code."),
     person("Kristina", "Lai", , "kwlai@ucdavis.edu", role = c("aut", "cre")),
     person("Chris", "Orwa", role = "aut"),
+    person("Kwan Ho", "Lee", , "ksjlee@ucdavis.edu", role = "ctb"),
     person("Kristen", "Aiemjoy", , "kaiemjoy@ucdavis.edu", role = "aut"),
     person("Douglas Ezra", "Morrison", , "demorrison@ucdavis.edu", role = "aut")
   )

NAMESPACE

@@ -5,6 +5,7 @@ S3method(autoplot,curve_params)
 S3method(autoplot,pop_data)
 S3method(autoplot,seroincidence)
 S3method(autoplot,seroincidence.by)
+S3method(autoplot,sim_results)
 S3method(autoplot,summary.seroincidence.by)
 S3method(print,seroincidence)
 S3method(print,seroincidence.by)
@@ -14,6 +15,7 @@ S3method(strata,default)
 S3method(summary,pop_data)
 S3method(summary,seroincidence)
 S3method(summary,seroincidence.by)
+export(analyze_sims)
 export(as_curve_params)
 export(as_noise_params)
 export(as_pop_data)

NEWS.md

@@ -2,6 +2,8 @@
 
 ## New features
 
+* Extended `sim_pop_data_multi()` to loop over multiple sample sizes (#444)
+* Added new functions `analyze_sims()` and `autoplot.sim_results()` (#444)
 * Rename `estimate_scr()` to `est_seroincidence_by()` (#439)
 * Rename `estimate_scr()` to `est_seroincidence()` (#432)
 * Rename argument `curve_params` to `sr_params` (#424)

R/analyze_sims.R

@@ -0,0 +1,69 @@
+#' Analyze simulation results
+#'
+#' @param data a [tibble::tbl_df] with columns:
+#' * `lambda.sim`,
+#' * `incidence.rate`,
+#' * `SE`,
+#' * `CI.lwr`,
+#' * `CI.upr`
+#' for example, as produced by [summary.seroincidence.by()] with
+#' `lambda.sim` as a stratifying variable
+#'
+#' @returns a `sim_results` object (extends [tibble::tbl_df])
+#' @export
+#'
+#' @example inst/examples/exm-analyze_sims.R
+#'
+analyze_sims <- function(
+    data) {
+
+  to_return <-
+    data |>
+    split(
+      f = ~  sample_size + lambda.sim
+    ) |>
+    lapply(FUN = analyze_sims_one_stratum) |>
+    bind_rows()
+
+  class(to_return) <- union("sim_results", class(to_return))
+
+  return(to_return)
+}
+
+analyze_sims_one_stratum <- function(
+    data,
+    true_lambda = data$lambda.sim,
+    sample_size = data$sample_size) {
+
+  # Filter out rows where CI.lwr or CI.upr is Inf or NaN
+  data <- data |>
+    filter(is.finite(.data$CI.lwr) & is.finite(.data$CI.upr))
+
+  # Compute Bias
+  bias <- mean(data$incidence.rate - true_lambda, na.rm = TRUE)
+
+  # Standard Error (Mean of reported standard errors)
+  standard_error <- mean(data$SE, na.rm = TRUE)
+
+  # RMSE (Root Mean Square Error)
+  rmse <- mean((data$incidence.rate - true_lambda)^2, na.rm = TRUE) |> sqrt()
+
+  # Confidence Interval Width (Mean of Upper - Lower bounds, without Inf values)
+  ci_width <- mean(data$CI.upr - data$CI.lwr, na.rm = TRUE)
+
+  coverage_prop <-
+    mean(data$CI.lwr <= true_lambda & data$CI.upr >= true_lambda, na.rm = TRUE)
+
+  to_return <- tibble(
+    lambda.sim = mean(true_lambda),
+    sample_size = mean(sample_size),
+    Bias = bias,
+    Mean_Est_SE = standard_error,
+    Empirical_SE = stats::sd(data$incidence.rate, na.rm = TRUE),
+    RMSE = rmse,
+    Mean_CI_Width = ci_width,
+    CI_Coverage = coverage_prop
+  )
+  # Return computed statistics as a list
+  return(to_return)
+}

R/autoplot.sim_results.R

@@ -0,0 +1,27 @@
+#' Plot simulation results
+#' `autoplot()` method for `sim_results` objects
+#'
+#' @param object a `sim_results` object (from [analyze_sims()])
+#' @param statistic which column of `object` should be the y-axis?
+#' @param ... unused
+#' @returns a [ggplot2::ggplot]
+#' @export
+#'
+#' @example inst/examples/exm-autoplot.sim_results.R
+autoplot.sim_results <- function(
+    object,
+    statistic = "Empirical_SE",
+    ...) {
+  object |>
+    dplyr::mutate(lambda.sim = factor(.data$lambda.sim)) |>
+    ggplot2::ggplot() +
+    ggplot2::aes(
+      x = .data$sample_size,
+      group = .data$lambda.sim,
+      col = .data$lambda.sim,
+      y = .data[[statistic]]
+    ) +
+    ggplot2::geom_point() +
+    ggplot2::geom_line() +
+    ggplot2::theme(legend.position = "bottom")
+}

R/sim_pop_data.R

@@ -106,7 +106,7 @@ sim_pop_data <- function(
 
   predpar <-
     curve_params |>
-    filter(.data$antigen_iso %in% antigen_isos) |>
+    dplyr::filter(.data$antigen_iso %in% antigen_isos) |>
     droplevels() |>
     prep_curve_params_for_array() |>
     df_to_array(dim_var_names = c("antigen_iso", "parameter"))
@@ -196,6 +196,7 @@ sim_pop_data <- function(
 #' consistent API.
 #' @keywords internal
 #' @export
+
 sim.cs <- function( # nolint: object_name_linter
   n.smpl, # nolint: object_name_linter
   age.rng, # nolint: object_name_linter

R/sim_pop_data_multi.R

@@ -3,53 +3,18 @@
 #' @param nclus  number of clusters
 #' @param rng_seed starting seed for random number generator,
 #' passed to [rngtools::RNGseq()]
-#' @param lambdas #incidence rate, in events/person*year
+#' @param sample_sizes sample sizes to simulate
+#' @param lambdas incidence rate, in events/person*year
 #' @param num_cores number of cores to use for parallel computations
 #' @param verbose whether to report verbose information
 #' @param ... arguments passed to [sim.cs()]
 #' @inheritDotParams sim_pop_data
 #' @return a [tibble::tibble()]
 #' @export
-#' @examples
-#' # Load curve parameters
-#' dmcmc <- typhoid_curves_nostrat_100
-#'
-#' # Specify the antibody-isotype responses to include in analyses
-#' antibodies <- c("HlyE_IgA", "HlyE_IgG")
-#'
-#' # Set seed to reproduce results
-#' set.seed(54321)
-#'
-#' # Simulated incidence rate per person-year
-#' lambdas = c(.05, .1, .15, .2, .3)
-#'
-#' # Range covered in simulations
-#' lifespan <- c(0, 10);
-#'
-#' # Cross-sectional sample size
-#' nrep <- 100
-#'
-#' # Biologic noise distribution
-#' dlims <- rbind(
-#'   "HlyE_IgA" = c(min = 0, max = 0.5),
-#'   "HlyE_IgG" = c(min = 0, max = 0.5)
-#' )
-#'
-#' sim_pop_data_multi(
-#'   curve_params = dmcmc,
-#'   lambdas = lambdas,
-#'   n_samples = nrep,
-#'   age_range = lifespan,
-#'   antigen_isos = antibodies,
-#'   n_mcmc_samples = 0,
-#'   renew_params = TRUE,
-#'   add_noise = TRUE,
-#'   noise_limits = dlims,
-#'   format = "long",
-#'   nclus = 10)
-#'
+#' @example inst/examples/exm-sim_pop_data_multi.R
 sim_pop_data_multi <- function(
     nclus = 10,
+    sample_sizes = 100,
     lambdas = c(.05, .1, .15, .2, .3),
     num_cores = max(1, parallel::detectCores() - 1),
     rng_seed = 1234,
@@ -81,31 +46,61 @@
   }
 
   doParallel::registerDoParallel(cores = num_cores)
-
+  n_sample_sizes <- length(sample_sizes)
   n_lambda <- length(lambdas)
 
   # trying to reproduce results using parallel
-  rng <- rngtools::RNGseq(n_lambda * nclus, rng_seed)
+  rng <- rngtools::RNGseq(n_sample_sizes * n_lambda * nclus, rng_seed)
+
+  dimnames1 <-
+    list(
+      "iteration" = 1:nclus,
+      "lambda" = lambdas,
+      "sample size" = sample_sizes
+    )
+
+  dims1 <-
+    sapply(FUN = length, dimnames1)
+
+  rng <- rng |>
+    array(
+      dimnames = dimnames1,
+      dim = dims1
+    )
   i <- NA
+  j <- 1
   r <- NA
 
   sim_df <-
+    foreach::foreach(
+      .combine = bind_rows,
+      j = seq_along(sample_sizes)
+
+    ) %:%
     foreach::foreach(
       .combine = bind_rows,
       i = seq_along(lambdas)
+
     ) %:%
     foreach::foreach(
       .combine = bind_rows,
       n = 1:nclus,
-      r = rng[(i - 1) * nclus + 1:nclus]
+      r = rng[1:nclus, i, j]
     ) %dopar% {
       l <- lambdas[i]
+      ns <- sample_sizes[j]
       rngtools::setRNG(r)
       sim_pop_data(
         lambda = l,
+        n_samples = ns,
         ...
       ) |>
-        mutate(lambda.sim = l, cluster = n)
+        mutate(
+          lambda.sim = l,
+          sample_size = ns,
+          cluster = n
+        ) |>
+        structure(r = r)
     }
   doParallel::stopImplicitCluster()
   sim_df <- sim_df |> set_biomarker_var(biomarker = "antigen_iso")

data-raw/.gitignore

@@ -1 +1,2 @@
-test_sim_020425_2.pdf
+*.pdf
+*.rmarkdown

data-raw/coverage-CI.R

@@ -0,0 +1,22 @@
+# Compute Coverage and its Confidence Interval
+compute_coverage_ci <- function(coverage_count, total_count) {
+  test_result <- stats::binom.test(
+    coverage_count,
+    total_count,
+    conf.level = 0.95
+  )
+  # 95% CI
+  coverage_proportion <- coverage_count / total_count
+  ci_lower <- test_result$conf.int[1]
+  ci_upper <- test_result$conf.int[2]
+
+  return(
+    list(
+      coverage = coverage_proportion,
+      ci_lower = ci_lower,
+      ci_upper = ci_upper
+    )
+  )
+}
+
+# coverage_result <- compute_coverage_ci(coverage_count, nrow(data)) # nolint: object_usage_linter

data-raw/generate_sim_table.R

@@ -0,0 +1,42 @@
+#' Generate table of simulation results
+#'
+#' @param results_list output from [simulate_seroincidence()]
+#' @param sample_size sample size of simulated data sets
+#' @noRd
+#' @returns a [tibble::tbl_df]
+generate_sim_table <- function(
+    results_list,
+    sample_size = results_list |> attr("sample_size")) {
+  # Initialize an empty list to store the results
+  summary_results <- list()
+
+  # Loop through each of the 100 results and extract the required columns
+  for (i in 1:300) {
+    # Extract the summary for each result
+    result_summary <- summary(results_list[[i]]$est1)
+
+    # Select the required columns
+    extracted_columns <- result_summary %>%
+      select(incidence.rate, SE, CI.lwr, CI.upr)
+
+    # Add a column for the index (optional, for tracking)
+    extracted_columns <- extracted_columns %>%
+      mutate(index = i)
+
+    # Append to the list
+    summary_results[[i]] <- extracted_columns
+  }
+
+  # Combine all results into a single data frame
+  final_table <- bind_rows(summary_results) %>%
+    mutate(sample_size = sample_size) # Add sample size column for clarity
+
+  final_table <-
+    final_table |>
+    structure(
+      sample_size = sample_size,
+      true_lambda = results_list |> attr("lambda_true")
+    )
+
+  return(final_table)
+}