1 Purpose

This document reproduces the descriptive analyses reported in the original firearm-storage publication and extends each analysis by comparing firearm owners with and without military experience.

The primary analyses are descriptive and survey weighted. They include:

  1. secure-container storage prevalence;
  2. secure-container storage across demographic and household characteristics;
  3. secure-container storage across firearm-carry frequency;
  4. current use of and willingness to adopt additional firearm-safety practices;
  5. motivations for changing firearm-storage behavior, stratified by intention to keep at least one firearm unlocked.

Inferential comparisons and adjusted models are placed in the final exploratory sections so that the central publication remains descriptive.

2 1. Packages

required_packages <- c(
  "foreign", "dplyr", "tidyr", "purrr", "survey",
  "broom", "ggplot2", "scales", "forcats", "gt",
  "patchwork", "knitr"
)

missing_packages <- required_packages[
  !required_packages %in% rownames(installed.packages())
]

if (length(missing_packages) > 0) {
  install.packages(missing_packages)
}

library(foreign)
library(dplyr)
library(tidyr)
library(purrr)
library(survey)
library(broom)
library(ggplot2)
library(scales)
library(forcats)
library(gt)
library(patchwork)
library(knitr)

3 2. Import and verify the dataset

The code will use an existing data frame named data when one is already loaded. Otherwise, it imports the SPSS file specified below.

data_file <- "KP_OMNI_2405_BMW_Client_File_03042024.sav"

data_exists <- exists(
  "data",
  envir = .GlobalEnv,
  inherits = FALSE
)

data_is_dataframe <- data_exists &&
  is.data.frame(get("data", envir = .GlobalEnv))

if (!data_is_dataframe) {

  if (!file.exists(data_file)) {
    stop(
      "The survey dataset was not found. Load a data frame named `data` ",
      "or update `data_file` to the correct location."
    )
  }

  data <- foreign::read.spss(
    data_file,
    to.data.frame = TRUE,
    use.value.labels = TRUE
  )
}

data <- as.data.frame(data)

required_columns <- c(
  "BMW1", "BMW3", "Status", "Weights",
  "ppage", "ppgender", "ppeduc5", "ppethm",
  "pphouse4", "ppinc7", "ppmarit5", "ppreg4",
  "ppemploy", "ppkid017", "xparty4", "xurbanicity",
  paste0("BMW4_", 1:7),
  paste0("BMW5_", 1:10),
  paste0("BMW6_", 1:10)
)

missing_columns <- setdiff(required_columns, names(data))

if (length(missing_columns) > 0) {
  stop(
    "The following required columns are missing: ",
    paste(missing_columns, collapse = ", ")
  )
}

4 3. Helper functions

# Convert checkbox fields to binary outcomes.
# Selected text = 1; unselected values such as 0 or blank = 0;
# skipped and missing values remain missing.
checkbox_binary <- function(x) {

  x_chr <- trimws(tolower(as.character(x)))

  dplyr::case_when(
    is.na(x) ~ NA_integer_,
    x_chr == "skipped" ~ NA_integer_,
    x_chr %in% c("", "0", "no", "not selected", "unchecked") ~ 0L,
    TRUE ~ 1L
  )
}

# Recode four-level agreement responses.
agreement_binary <- function(x) {

  x_chr <- trimws(tolower(as.character(x)))

  dplyr::case_when(
    x_chr %in% c("strongly agree", "somewhat agree") ~ 1L,
    x_chr %in% c("somewhat disagree", "strongly disagree") ~ 0L,
    TRUE ~ NA_integer_
  )
}

agreement_four <- function(x) {

  x_chr <- trimws(tolower(as.character(x)))

  dplyr::case_when(
    x_chr == "strongly disagree" ~ "Strongly disagree",
    x_chr == "somewhat disagree" ~ "Somewhat disagree",
    x_chr == "somewhat agree" ~ "Somewhat agree",
    x_chr == "strongly agree" ~ "Strongly agree",
    TRUE ~ NA_character_
  )
}

format_p <- function(x) {
  dplyr::case_when(
    is.na(x) ~ NA_character_,
    x < .001 ~ "<.001",
    TRUE ~ sprintf("%.3f", x)
  )
}

format_percent_ci <- function(est, lower, upper, digits = 1) {
  paste0(
    sprintf(paste0("%.", digits, "f"), 100 * est),
    "% (",
    sprintf(paste0("%.", digits, "f"), 100 * lower),
    "–",
    sprintf(paste0("%.", digits, "f"), 100 * upper),
    "%)"
  )
}

safe_svychisq <- function(formula, design) {
  tryCatch(
    survey::svychisq(formula, design, statistic = "F"),
    error = function(e) NULL
  )
}

extract_test_p <- function(test_object) {
  if (is.null(test_object)) return(NA_real_)
  unname(test_object$p.value)
}

5 4. Prepare the firearm-owner sample

Military experience is defined using the original Status item. Any response other than “None of the above” is classified as military experience, consistent with the prior analysis.

owners <- data %>%
  dplyr::filter(
    BMW1 == "Yes",
    !is.na(Status),
    !is.na(Weights)
  ) %>%
  dplyr::mutate(
    Respondent_ID = dplyr::row_number(),

    Military = dplyr::if_else(
      trimws(as.character(Status)) == "None of the above",
      "No Military Experience",
      "Military Experience"
    ),

    Military = factor(
      Military,
      levels = c(
        "No Military Experience",
        "Military Experience"
      )
    ),

    Secure_Container = checkbox_binary(BMW5_2),

    age.cat = cut(
      ppage,
      breaks = c(18, 25, 40, 60, Inf),
      labels = c("18–25", "26–40", "41–60", ">60"),
      right = FALSE
    ),

    Children = dplyr::case_when(
      is.na(ppkid017) ~ NA_character_,
      ppkid017 > 0 ~ "Children",
      ppkid017 == 0 ~ "No Children"
    ),

    Children = factor(
      Children,
      levels = c("No Children", "Children")
    ),

    Always_Unlocked = dplyr::case_when(
      BMW4_7 %in% c("Strongly agree", "Somewhat agree") ~
        "Endorses always-unlocked storage",
      BMW4_7 %in% c("Strongly disagree", "Somewhat disagree") ~
        "Rejects always-unlocked storage",
      TRUE ~ NA_character_
    ),

    Always_Unlocked = factor(
      Always_Unlocked,
      levels = c(
        "Rejects always-unlocked storage",
        "Endorses always-unlocked storage"
      )
    )
  )

for (j in 1:10) {
  owners[[paste0("Use_", j)]] <-
    checkbox_binary(owners[[paste0("BMW5_", j)]])

  owners[[paste0("Willing_", j)]] <-
    checkbox_binary(owners[[paste0("BMW6_", j)]])
}

for (j in 1:7) {
  owners[[paste0("Attitude_", j)]] <-
    agreement_binary(owners[[paste0("BMW4_", j)]])
}

design_owners <- survey::svydesign(
  ids = ~1,
  weights = ~Weights,
  data = owners
)

cat("Total firearm owners:", nrow(owners), "\n")
## Total firearm owners: 336
print(table(owners$Military, useNA = "ifany"))
## 
## No Military Experience    Military Experience 
##                    282                     54

6 5. Sample composition by military experience

sample_counts <- owners %>%
  dplyr::count(Military, name = "Unweighted_n")

weighted_group_share <- survey::svymean(
  ~Military,
  design_owners,
  na.rm = TRUE
)

sample_counts
weighted_group_share
##                                   mean     SE
## MilitaryNo Military Experience 0.83836 0.0213
## MilitaryMilitary Experience    0.16164 0.0213

7 6. Secure-container storage by military experience

This reproduces the primary prevalence outcome from the original paper and compares the weighted prevalence between military groups.

secure_by_military <- survey::svyby(
  ~Secure_Container,
  ~Military,
  design = design_owners,
  FUN = survey::svymean,
  vartype = c("se", "ci"),
  na.rm = TRUE,
  keep.names = FALSE
) %>%
  as.data.frame() %>%
  dplyr::mutate(
    Weighted_Percent = 100 * Secure_Container,
    CI_Lower_Percent = 100 * ci_l,
    CI_Upper_Percent = 100 * ci_u,
    Estimate = format_percent_ci(
      Secure_Container,
      ci_l,
      ci_u
    )
  )

secure_test <- safe_svychisq(
  ~Military + factor(Secure_Container),
  design_owners
)

secure_by_military %>%
  dplyr::select(
    Military,
    Weighted_Percent,
    CI_Lower_Percent,
    CI_Upper_Percent,
    Estimate
  ) %>%
  gt() %>%
  tab_header(
    title = "Secure-container storage by military experience",
    subtitle = paste0(
      "Survey-weighted prevalence; Rao–Scott p = ",
      format_p(extract_test_p(secure_test))
    )
  ) %>%
  cols_label(
    Military = "Military experience",
    Weighted_Percent = "Weighted %",
    CI_Lower_Percent = "95% CI lower",
    CI_Upper_Percent = "95% CI upper",
    Estimate = "Weighted % (95% CI)"
  ) %>%
  fmt_number(
    columns = c(
      Weighted_Percent,
      CI_Lower_Percent,
      CI_Upper_Percent
    ),
    decimals = 1
  )
Secure-container storage by military experience
Survey-weighted prevalence; Rao–Scott p = NA
Military experience Weighted % 95% CI lower 95% CI upper Weighted % (95% CI)
No Military Experience 57.2 51.1 63.4 57.2% (51.1–63.4%)
Military Experience 67.1 53.9 80.3 67.1% (53.9–80.3%)
 
ggplot(
  secure_by_military,
  aes(x = Military, y = Weighted_Percent)
) +
  geom_col(width = 0.62) +
  geom_errorbar(
    aes(
      ymin = CI_Lower_Percent,
      ymax = CI_Upper_Percent
    ),
    width = 0.12
  ) +
  geom_text(
    aes(label = sprintf("%.1f%%", Weighted_Percent)),
    vjust = -0.6,
    size = 4
  ) +
  scale_y_continuous(
    limits = c(0, 100),
    labels = function(x) paste0(x, "%"),
    expand = expansion(mult = c(0, .08))
  ) +
  labs(
    x = NULL,
    y = "Survey-weighted prevalence",
    title = "Secure-container storage by military experience"
  ) +
  theme_classic(base_size = 12)

8 7. Demographic and household characteristics

The original publication reported secure-container prevalence across demographic and household categories. This section reproduces those weighted prevalence estimates separately for military and nonmilitary firearm owners.

Because some military subgroups are small, these results should be treated as descriptive. Cells with fewer than five unweighted respondents are flagged.

demographic_variables <- c(
  "ppgender",
  "age.cat",
  "xparty4",
  "xurbanicity",
  "ppeduc5",
  "ppethm",
  "pphouse4",
  "ppinc7",
  "ppmarit5",
  "ppreg4",
  "ppemploy",
  "Children"
)

demographic_labels <- c(
  ppgender = "Gender",
  age.cat = "Age",
  xparty4 = "Political affiliation",
  xurbanicity = "Urbanicity",
  ppeduc5 = "Educational attainment",
  ppethm = "Race and ethnicity",
  pphouse4 = "Housing type",
  ppinc7 = "Household income",
  ppmarit5 = "Marital status",
  ppreg4 = "Geographic region",
  ppemploy = "Employment",
  Children = "Children in household"
)
weighted_secure_by_demographic <- function(var_name) {

  analysis_subset <- owners %>%
    dplyr::filter(
      !is.na(.data[[var_name]]),
      !is.na(Secure_Container),
      !is.na(Military),
      !is.na(Weights)
    ) %>%
    droplevels()

  design_subset <- survey::svydesign(
    ids = ~1,
    weights = ~Weights,
    data = analysis_subset
  )

  weighted <- survey::svyby(
    ~Secure_Container,
    as.formula(paste0("~Military + `", var_name, "`")),
    design = design_subset,
    FUN = survey::svymean,
    vartype = c("se", "ci"),
    na.rm = TRUE,
    keep.names = FALSE
  ) %>%
    as.data.frame()

  names(weighted)[names(weighted) == var_name] <- "Category"

  counts <- analysis_subset %>%
    dplyr::count(
      Military,
      Category = .data[[var_name]],
      name = "Unweighted_n"
    )

  weighted %>%
    dplyr::left_join(
      counts,
      by = c("Military", "Category")
    ) %>%
    dplyr::mutate(
      Characteristic = unname(demographic_labels[var_name]),
      Weighted_Percent = 100 * Secure_Container,
      CI_Lower_Percent = 100 * ci_l,
      CI_Upper_Percent = 100 * ci_u,
      Estimate = format_percent_ci(
        Secure_Container,
        ci_l,
        ci_u
      ),
      Sparse_Cell = Unweighted_n < 5
    ) %>%
    dplyr::select(
      Characteristic,
      Category,
      Military,
      Unweighted_n,
      Weighted_Percent,
      CI_Lower_Percent,
      CI_Upper_Percent,
      Estimate,
      Sparse_Cell
    )
}

demographic_secure_results <- purrr::map_dfr(
  demographic_variables,
  weighted_secure_by_demographic
)

demographic_secure_results
demographic_table <- demographic_secure_results %>%
  dplyr::mutate(
    Display = paste0(
      Estimate,
      dplyr::if_else(
        Sparse_Cell,
        "†",
        ""
      )
    )
  ) %>%
  dplyr::select(
    Characteristic,
    Category,
    Military,
    Display
  ) %>%
  tidyr::pivot_wider(
    names_from = Military,
    values_from = Display
  )

demographic_table %>%
  gt(groupname_col = "Characteristic") %>%
  tab_header(
    title = "Secure-container storage across demographic characteristics",
    subtitle = "Survey-weighted prevalence by military experience"
  ) %>%
  cols_label(
    Category = "Category",
    `No Military Experience` = "No military experience",
    `Military Experience` = "Military experience"
  ) %>%
  tab_source_note(
    "Values are weighted percentages with 95% confidence intervals. † Unweighted cell n < 5; interpret cautiously."
  )
Secure-container storage across demographic characteristics
Survey-weighted prevalence by military experience
Category No military experience Military experience
Gender
Male 56.5% (48.1–64.8%) 65.4% (51.8–79.0%)
Female 58.1% (49.1–67.1%) 100.0% (100.0–100.0%)†
Age
18–25 68.6% (49.2–88.0%) 100.0% (100.0–100.0%)†
26–40 65.9% (54.5–77.3%) 64.6% (8.3–120.9%)†
41–60 59.4% (49.7–69.1%) 74.7% (54.9–94.6%)
>60 45.0% (33.4–56.5%) 61.8% (43.2–80.4%)
Political affiliation
Republican 60.6% (51.2–70.0%) 58.7% (34.3–83.2%)
Democrat 53.1% (38.2–68.0%) 88.7% (67.0–110.4%)
Independent 54.1% (43.4–64.8%) 65.4% (45.8–85.0%)
Something else 60.9% (40.3–81.5%) 57.2% (9.0–105.5%)†
Urbanicity
Urban 54.3% (41.9–66.7%) 72.4% (50.2–94.7%)
Rural 54.6% (43.1–66.0%) 75.1% (51.9–98.2%)
Suburban 60.4% (51.5–69.3%) 59.6% (39.6–79.7%)
Educational attainment
No high school diploma or GED 29.5% (6.8–52.2%) 100.0% (100.0–100.0%)†
High school graduate (high school diploma or the equivalent GED) 57.9% (47.2–68.6%) 49.6% (15.0–84.1%)
Some college or Associate degree 50.2% (38.3–62.0%) 73.0% (52.9–93.2%)
Bachelor’s degree 67.4% (55.0–79.8%) 61.6% (34.6–88.7%)
Master’s degree or above 71.4% (56.8–86.0%) 76.5% (48.4–104.5%)
Race and ethnicity
White, Non-Hispanic 57.1% (50.4–63.8%) 64.2% (49.2–79.3%)
Black or African American, Non-Hispanic 46.0% (21.2–70.8%) 61.7% (19.1–104.4%)
Other, Non-Hispanic 86.2% (60.5–111.9%) 100.0% (100.0–100.0%)†
Hispanic 58.0% (37.3–78.7%) 82.4% (50.4–114.4%)
2+ races, Non-Hispanic 46.1% (18.3–73.9%) 0.0% (0.0–0.0%)†
Housing type
One-family house detached from any other house 59.8% (53.2–66.5%) 65.2% (51.1–79.4%)
One-family condo or townhouse attached to other units 53.7% (28.0–79.3%) 72.4% (23.6–121.2%)†
Building with 2 or more apartments 34.6% (13.2–56.0%) 100.0% (100.0–100.0%)†
Other (mobile home, boat, RV, van, etc.) 44.0% (12.6–75.3%) 100.0% (100.0–100.0%)†
Household income
Under $10,000 56.8% (19.3–94.3%) NA
$10,000 to $24,999 42.1% (21.4–62.8%) 0.0% (0.0–0.0%)†
$25,000 to $49,999 41.7% (23.9–59.4%) 65.2% (24.1–106.3%)
$50,000 to $74,999 48.2% (34.0–62.5%) 66.0% (34.5–97.6%)
$75,000 to $99,999 51.0% (34.3–67.7%) 60.9% (18.6–103.2%)
$100,000 to $149,999 64.1% (51.8–76.3%) 100.0% (100.0–100.0%)
$150,000 or more 76.3% (65.2–87.4%) 60.0% (38.9–81.2%)
Marital status
Now married 62.1% (54.3–69.9%) 69.0% (53.7–84.2%)
Widowed 17.5% (-5.8–40.8%) 63.3% (27.6–99.1%)
Divorced 56.0% (37.8–74.3%) 63.6% (20.0–107.2%)
Separated 67.9% (18.4–117.5%)† NA
Never married 52.3% (39.9–64.7%) 46.4% (-22.6–115.4%)†
Geographic region
Northeast 65.3% (48.9–81.7%) 68.8% (32.3–105.4%)
Midwest 67.7% (56.0–79.4%) 48.7% (17.7–79.7%)
South 48.5% (39.2–57.8%) 71.2% (54.2–88.1%)
West 58.9% (44.7–73.0%) 74.2% (42.0–106.4%)
Employment
Working full-time 61.9% (53.9–70.0%) 73.3% (53.4–93.2%)
Working part-time 71.1% (53.7–88.4%) 79.0% (42.5–115.5%)
Not working 46.8% (36.5–57.2%) 60.2% (41.0–79.4%)
Children in household
No Children 53.7% (46.4–61.1%) 58.2% (42.0–74.4%)
Children 65.8% (55.0–76.6%) 92.9% (79.3–106.4%)
Values are weighted percentages with 95% confidence intervals. † Unweighted cell n < 5; interpret cautiously.

9 8. Frequency of carrying a loaded firearm

This reproduces the carry-frequency analysis from the original publication and adds separate estimates for military and nonmilitary firearm owners.

carry_levels_original <- c(
  "I own a firearm but never carry it loaded",
  "Almost never",
  "At least once a year, but not every month",
  "At least once a month, but not every week",
  "At least once a week, but not every day",
  "Almost every day",
  "Daily"
)

carry_labels <- c(
  "I own a firearm but never carry it loaded" = "Never",
  "Almost never" = "Almost never",
  "At least once a year, but not every month" = "Annually",
  "At least once a month, but not every week" = "Monthly",
  "At least once a week, but not every day" = "Weekly",
  "Almost every day" = "Daily",
  "Daily" = "Daily"
)

carry_data <- owners %>%
  dplyr::filter(
    !is.na(BMW3),
    BMW3 != "Skipped",
    !is.na(Secure_Container),
    !is.na(Military)
  ) %>%
  dplyr::mutate(
    Carry = dplyr::recode(
      as.character(BMW3),
      !!!carry_labels,
      .default = NA_character_
    ),

    Carry = factor(
      Carry,
      levels = c(
        "Never",
        "Almost never",
        "Annually",
        "Monthly",
        "Weekly",
        "Daily"
      )
    )
  ) %>%
  dplyr::filter(!is.na(Carry)) %>%
  droplevels()

design_carry <- survey::svydesign(
  ids = ~1,
  weights = ~Weights,
  data = carry_data
)
carry_estimates <- survey::svyby(
  ~Secure_Container,
  ~Military + Carry,
  design = design_carry,
  FUN = survey::svymean,
  vartype = c("se", "ci"),
  na.rm = TRUE,
  keep.names = FALSE
) %>%
  as.data.frame() %>%
  dplyr::mutate(
    Percent = 100 * Secure_Container,
    CI_Lower = 100 * ci_l,
    CI_Upper = 100 * ci_u
  )

carry_counts <- carry_data %>%
  dplyr::count(Military, Carry, name = "Unweighted_n")

carry_estimates <- carry_estimates %>%
  dplyr::left_join(
    carry_counts,
    by = c("Military", "Carry")
  )

carry_estimates %>%
  dplyr::mutate(
    Estimate = sprintf(
      "%.1f%% (%.1f–%.1f)",
      Percent,
      CI_Lower,
      CI_Upper
    )
  ) %>%
  dplyr::select(
    Carry,
    Military,
    Unweighted_n,
    Estimate
  ) %>%
  tidyr::pivot_wider(
    names_from = Military,
    values_from = c(Unweighted_n, Estimate)
  ) %>%
  gt() %>%
  tab_header(
    title = "Secure-container storage by carry frequency and military experience"
  )
Secure-container storage by carry frequency and military experience
Carry Unweighted_n_No Military Experience Unweighted_n_Military Experience Estimate_No Military Experience Estimate_Military Experience
Never 134 20 54.2% (45.4–62.9) 57.4% (34.4–80.4)
Almost never 64 16 56.6% (43.7–69.5) 64.8% (40.5–89.1)
Annually 12 2 53.8% (25.1–82.6) 100.0% (100.0–100.0)
Monthly 13 4 68.8% (42.9–94.7) 76.0% (34.4–117.6)
Weekly 23 3 76.2% (57.5–94.8) 75.0% (30.0–120.1)
Daily 27 7 55.7% (36.3–75.0) 78.6% (50.0–107.1)
 
ggplot(
  carry_estimates,
  aes(
    x = Carry,
    y = Percent,
    group = Military,
    shape = Military,
    linetype = Military
  )
) +
  geom_line(
    position = position_dodge(width = 0.16)
  ) +
  geom_point(
    size = 3,
    position = position_dodge(width = 0.16)
  ) +
  geom_errorbar(
    aes(
      ymin = CI_Lower,
      ymax = CI_Upper
    ),
    width = 0.08,
    position = position_dodge(width = 0.16)
  ) +
  geom_text(
    aes(
      label = paste0(
        sprintf("%.0f%%", Percent),
        "\n(n=",
        Unweighted_n,
        ")"
      )
    ),
    size = 3,
    vjust = -1.0,
    position = position_dodge(width = 0.16),
    check_overlap = TRUE
  ) +
  scale_y_continuous(
    limits = c(0, 110),
    breaks = seq(0, 100, 25),
    labels = function(x) paste0(x, "%")
  ) +
  labs(
    x = "Frequency of carrying a loaded firearm",
    y = "Survey-weighted prevalence of secure-container storage",
    shape = "Military experience",
    linetype = "Military experience",
    title = "Secure-container storage by carry frequency"
  ) +
  theme_classic(base_size = 12) +
  theme(
    legend.position = "top",
    axis.text.x = element_text(
      angle = 25,
      hjust = 1
    )
  )

10 9. Current use of firearm-safety practices

method_labels <- c(
  Use_1 = "Keep all firearms unloaded",
  Use_2 = "Use a secure container",
  Use_3 = "Use a cable or trigger lock",
  Use_4 = "Use an access alarm or notification",
  Use_5 = "Use a firearm sensor",
  Use_6 = "Lock ammunition separately",
  Use_7 = "Disassemble firearms",
  Use_8 = "Entrust keys or parts to another",
  Use_9 = "Use another safety measure",
  Use_10 = "None of the above"
)

willing_labels <- c(
  Willing_1 = "Keep all firearms unloaded",
  Willing_2 = "Use a secure container",
  Willing_3 = "Use a cable or trigger lock",
  Willing_4 = "Use an access alarm or notification",
  Willing_5 = "Use a firearm sensor",
  Willing_6 = "Lock ammunition separately",
  Willing_7 = "Disassemble firearms",
  Willing_8 = "Entrust keys or parts to another",
  Willing_9 = "Use another safety measure",
  Willing_10 = "None of the above"
)
current_use_long <- owners %>%
  dplyr::select(
    Respondent_ID,
    Military,
    Secure_Container,
    Weights,
    dplyr::all_of(paste0("Use_", 1:10))
  ) %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("Use_"),
    names_to = "Outcome",
    values_to = "Endorsed"
  ) %>%
  dplyr::mutate(
    Method = unname(method_labels[Outcome]),
    Secure_Group = dplyr::if_else(
      Secure_Container == 1,
      "Uses secure container",
      "Does not use secure container"
    )
  )
estimate_weighted_binary <- function(df, group_variables) {

  design_temp <- survey::svydesign(
    ids = ~1,
    weights = ~Weights,
    data = df
  )

  grouping_formula <- as.formula(
    paste0(
      "~",
      paste(group_variables, collapse = " + ")
    )
  )

  survey::svyby(
    ~Endorsed,
    grouping_formula,
    design = design_temp,
    FUN = survey::svymean,
    vartype = c("se", "ci"),
    na.rm = TRUE,
    keep.names = FALSE
  ) %>%
    as.data.frame() %>%
    dplyr::mutate(
      Percent = 100 * Endorsed,
      CI_Lower = 100 * ci_l,
      CI_Upper = 100 * ci_u
    )
}

current_all <- estimate_weighted_binary(
  current_use_long,
  c("Military", "Outcome", "Method")
) %>%
  dplyr::mutate(Population = "All firearm owners")

current_nonsecure <- current_use_long %>%
  dplyr::filter(Secure_Container == 0) %>%
  estimate_weighted_binary(
    c("Military", "Outcome", "Method")
  ) %>%
  dplyr::mutate(
    Population = "Firearm owners not using a secure container"
  )

current_results <- dplyr::bind_rows(
  current_nonsecure,
  current_all
)

current_results

11 10. Willingness to adopt additional firearm-safety practices

The original survey used display logic: respondents who already used a method were not asked whether they would consider it. The primary willingness estimate below therefore uses respondents who did not currently use the corresponding method and provided a substantive BMW6 response.

A second “current or willing” estimate is also calculated to match the interpretation in the original publication.

willingness_long <- purrr::map_dfr(
  1:10,
  function(j) {

    use_var <- paste0("Use_", j)
    willing_var <- paste0("Willing_", j)

    owners %>%
      dplyr::transmute(
        Respondent_ID,
        Military,
        Secure_Container,
        Weights,
        Outcome = willing_var,
        Method = unname(willing_labels[willing_var]),
        Current = .data[[use_var]],
        Willing = .data[[willing_var]]
      )
  }
)
estimate_willingness <- function(df, population_label) {

  eligible <- df %>%
    dplyr::filter(
      Current == 0,
      !is.na(Willing)
    ) %>%
    dplyr::rename(Endorsed = Willing)

  estimates <- estimate_weighted_binary(
    eligible,
    c("Military", "Outcome", "Method")
  )

  counts <- eligible %>%
    dplyr::count(
      Military,
      Outcome,
      Method,
      name = "Unweighted_n"
    )

  estimates %>%
    dplyr::left_join(
      counts,
      by = c("Military", "Outcome", "Method")
    ) %>%
    dplyr::mutate(Population = population_label)
}

willing_all <- estimate_willingness(
  willingness_long,
  "All firearm owners not currently using each method"
)

willing_nonsecure <- willingness_long %>%
  dplyr::filter(Secure_Container == 0) %>%
  estimate_willingness(
    "Secure-container nonusers not currently using each method"
  )

willing_results <- dplyr::bind_rows(
  willing_nonsecure,
  willing_all
)

willing_results

12 11. Publication table: current use and willingness by military experience

This table mirrors Table 2 of the original paper but presents separate columns for military and nonmilitary respondents.

make_method_table <- function(
  current_df,
  willing_df,
  population_name
) {

  current_part <- current_df %>%
    dplyr::filter(Population == population_name) %>%
    dplyr::select(
      Method,
      Military,
      Current_Percent = Percent
    )

  willingness_name <- dplyr::case_when(
    population_name ==
      "All firearm owners" ~
      "All firearm owners not currently using each method",

    population_name ==
      "Firearm owners not using a secure container" ~
      "Secure-container nonusers not currently using each method"
  )

  willing_part <- willing_df %>%
    dplyr::filter(Population == willingness_name) %>%
    dplyr::select(
      Method,
      Military,
      Willing_Percent = Percent
    )

  current_part %>%
    dplyr::full_join(
      willing_part,
      by = c("Method", "Military")
    ) %>%
    tidyr::pivot_wider(
      names_from = Military,
      values_from = c(
        Current_Percent,
        Willing_Percent
      )
    ) %>%
    dplyr::mutate(Population = population_name)
}

table2_nonsecure <- make_method_table(
  current_results,
  willing_results,
  "Firearm owners not using a secure container"
)

table2_all <- make_method_table(
  current_results,
  willing_results,
  "All firearm owners"
)

table2_combined <- dplyr::bind_rows(
  table2_nonsecure,
  table2_all
)

table2_combined %>%
  gt(groupname_col = "Population") %>%
  tab_header(
    title = "Current use and willingness to adopt additional firearm-safety measures",
    subtitle = "Survey-weighted percentages by military experience"
  ) %>%
  cols_label(
    Method = "Safety measure",
    `Current_Percent_No Military Experience` =
      "Current: no military",
    `Current_Percent_Military Experience` =
      "Current: military",
    `Willing_Percent_No Military Experience` =
      "Willing: no military",
    `Willing_Percent_Military Experience` =
      "Willing: military"
  ) %>%
  fmt_number(
    columns = where(is.numeric),
    decimals = 1,
    pattern = "{x}%"
  ) %>%
  tab_source_note(
    "Willingness estimates include only respondents not currently using the corresponding method. Current-use and willingness percentages should therefore not be interpreted as estimates from identical denominators."
  )
Current use and willingness to adopt additional firearm-safety measures
Survey-weighted percentages by military experience
Safety measure Current: no military Current: military Willing: no military Willing: military
Firearm owners not using a secure container
Disassemble firearms 4.8% 11.0% 3.9% 0.0%
Entrust keys or parts to another 1.7% 0.0% 2.5% 0.0%
Keep all firearms unloaded 41.0% 59.5% 26.5% 22.6%
Lock ammunition separately 23.0% 40.3% 22.5% 6.9%
None of the above 38.3% 31.6% 20.6% 59.8%
Use a cable or trigger lock 15.5% 15.9% 15.5% 4.5%
Use a firearm sensor 1.0% 0.0% 6.1% 0.0%
Use a secure container 0.0% 0.0% 52.3% 28.5%
Use an access alarm or notification 0.0% 4.7% 13.1% 2.9%
Use another safety measure 4.0% 5.5% 5.3% 0.0%
All firearm owners
Disassemble firearms 5.8% 8.8% 5.8% 6.6%
Entrust keys or parts to another 5.5% 3.6% 4.1% 4.6%
Keep all firearms unloaded 50.5% 43.8% 24.6% 17.6%
Lock ammunition separately 36.4% 35.7% 24.4% 13.4%
None of the above 16.4% 10.4% 32.5% 44.9%
Use a cable or trigger lock 18.4% 26.5% 24.1% 24.9%
Use a firearm sensor 0.4% 2.3% 9.0% 3.4%
Use a secure container 57.2% 67.1% 52.3% 28.5%
Use an access alarm or notification 1.4% 8.0% 19.9% 9.5%
Use another safety measure 3.0% 1.8% 4.7% 1.8%
Willingness estimates include only respondents not currently using the corresponding method. Current-use and willingness percentages should therefore not be interpreted as estimates from identical denominators. 
plot_method_data <- table2_all %>%
  dplyr::select(-Population) %>%
  tidyr::pivot_longer(
    cols = -Method,
    names_to = c("Measure", "Military"),
    names_pattern =
      "(Current_Percent|Willing_Percent)_(.*)",
    values_to = "Percent"
  ) %>%
  dplyr::mutate(
    Measure = dplyr::recode(
      Measure,
      Current_Percent = "Currently uses",
      Willing_Percent = "Would consider"
    )
  )

ggplot(
  plot_method_data,
  aes(
    x = Percent,
    y = forcats::fct_rev(Method),
    shape = Military
  )
) +
  geom_point(
    size = 2.8,
    position = position_dodge(width = .55)
  ) +
  facet_wrap(
    ~Measure,
    ncol = 2
  ) +
  scale_x_continuous(
    limits = c(0, 100),
    labels = function(x) paste0(x, "%")
  ) +
  labs(
    x = "Survey-weighted percentage",
    y = NULL,
    shape = "Military experience",
    title = "Current use and willingness to adopt firearm-safety measures"
  ) +
  theme_classic(base_size = 12) +
  theme(
    legend.position = "top",
    strip.background = element_blank(),
    strip.text = element_text(face = "bold")
  )

13 12. Motivators for changing firearm-storage behavior

The original publication stratified motivators according to whether respondents endorsed keeping at least one firearm unlocked. The extension below further stratifies the estimates by military experience, producing four combinations:

  • no military experience / rejects always-unlocked storage;
  • no military experience / endorses always-unlocked storage;
  • military experience / rejects always-unlocked storage;
  • military experience / endorses always-unlocked storage.
motivator_labels <- c(
  BMW4_1 = "If I had children at home",
  BMW4_2 = "Household mental or physical health concerns",
  BMW4_3 = "A close friend or family member asked",
  BMW4_4 = "Prevent use, theft, or damage",
  BMW4_5 = "Prevent accidental injury",
  BMW4_6 = "Prevent suicide for me or others"
)

response_levels <- c(
  "Strongly disagree",
  "Somewhat disagree",
  "Somewhat agree",
  "Strongly agree"
)
motivator_long <- owners %>%
  dplyr::select(
    Respondent_ID,
    Military,
    Always_Unlocked,
    Weights,
    dplyr::all_of(paste0("BMW4_", 1:6))
  ) %>%
  tidyr::pivot_longer(
    cols = dplyr::all_of(paste0("BMW4_", 1:6)),
    names_to = "Question",
    values_to = "Response"
  ) %>%
  dplyr::mutate(
    Response = agreement_four(Response),
    Response = factor(
      Response,
      levels = response_levels
    ),
    Motivator = unname(motivator_labels[Question]),
    Top_Two = dplyr::if_else(
      Response %in% c(
        "Somewhat agree",
        "Strongly agree"
      ),
      1L,
      0L,
      missing = NA_integer_
    )
  ) %>%
  dplyr::filter(
    !is.na(Response),
    !is.na(Always_Unlocked),
    !is.na(Military)
  )
design_motivators <- survey::svydesign(
  ids = ~1,
  weights = ~Weights,
  data = motivator_long
)

motivator_top_two <- survey::svyby(
  ~Top_Two,
  ~Military + Always_Unlocked + Question + Motivator,
  design = design_motivators,
  FUN = survey::svymean,
  vartype = c("se", "ci"),
  na.rm = TRUE,
  keep.names = FALSE
) %>%
  as.data.frame() %>%
  dplyr::mutate(
    Percent = 100 * Top_Two,
    CI_Lower = 100 * ci_l,
    CI_Upper = 100 * ci_u
  )

motivator_counts <- motivator_long %>%
  dplyr::distinct(
    Respondent_ID,
    Military,
    Always_Unlocked
  ) %>%
  dplyr::count(
    Military,
    Always_Unlocked,
    name = "Unweighted_n"
  )

motivator_top_two <- motivator_top_two %>%
  dplyr::left_join(
    motivator_counts,
    by = c("Military", "Always_Unlocked")
  )

motivator_top_two
motivator_top_two %>%
  dplyr::mutate(
    Group = paste(
      Military,
      Always_Unlocked,
      sep = ": "
    ),
    Estimate = paste0(
      sprintf("%.1f%%", Percent),
      " (",
      sprintf("%.1f", CI_Lower),
      "–",
      sprintf("%.1f", CI_Upper),
      ")"
    )
  ) %>%
  dplyr::select(
    Motivator,
    Group,
    Estimate
  ) %>%
  tidyr::pivot_wider(
    names_from = Group,
    values_from = Estimate
  ) %>%
  gt() %>%
  tab_header(
    title = "Motivators for changing firearm-storage behavior",
    subtitle = "Top-two-box agreement by military experience and always-unlocked attitude"
  ) %>%
  tab_source_note(
    "Values are survey-weighted percentages with 95% confidence intervals."
  )
Motivators for changing firearm-storage behavior
Top-two-box agreement by military experience and always-unlocked attitude
Motivator No Military Experience: Rejects always-unlocked storage Military Experience: Rejects always-unlocked storage No Military Experience: Endorses always-unlocked storage Military Experience: Endorses always-unlocked storage
A close friend or family member asked 92.5% (87.7–97.4) 72.4% (53.5–91.4) 68.3% (60.5–76.1) 38.9% (21.1–56.7)
Household mental or physical health concerns 98.0% (95.1–100.8) 91.6% (80.3–103.0) 93.5% (89.3–97.7) 92.1% (82.8–101.4)
If I had children at home 98.7% (96.9–100.5) 87.7% (74.4–101.1) 89.0% (83.8–94.2) 88.2% (76.4–100.1)
Prevent accidental injury 96.8% (94.0–99.6) 85.6% (71.8–99.4) 80.0% (73.5–86.6) 58.0% (39.8–76.1)
Prevent suicide for me or others 94.6% (90.9–98.3) 89.4% (77.7–101.1) 77.7% (70.7–84.7) 72.3% (56.8–87.9)
Prevent use, theft, or damage 97.2% (94.4–100.0) 88.5% (76.0–101.0) 80.9% (74.3–87.4) 69.1% (51.5–86.7)
Values are survey-weighted percentages with 95% confidence intervals. 
ggplot(
  motivator_top_two,
  aes(
    x = Percent,
    y = forcats::fct_reorder(Motivator, Percent),
    shape = Military
  )
) +
  geom_errorbar(
    aes(
      xmin = CI_Lower,
      xmax = CI_Upper
    ),
    width = .14,
    position = position_dodge(width = .55)
  ) +
  geom_point(
    size = 2.8,
    position = position_dodge(width = .55)
  ) +
  facet_wrap(
    ~Always_Unlocked,
    ncol = 1
  ) +
  scale_x_continuous(
    limits = c(0, 100),
    labels = function(x) paste0(x, "%")
  ) +
  labs(
    x = "Survey-weighted top-two-box agreement",
    y = NULL,
    shape = "Military experience",
    title = "Motivations for changing firearm-storage behavior"
  ) +
  theme_classic(base_size = 12) +
  theme(
    legend.position = "top",
    strip.background = element_blank(),
    strip.text = element_text(face = "bold")
  )

13.1 Response-intensity distributions

The following plot reproduces the original distinction between “somewhat agree” and “strongly agree,” while adding military experience.

motivator_stacked <- motivator_long %>%
  dplyr::filter(
    Response %in% c(
      "Somewhat agree",
      "Strongly agree"
    )
  ) %>%
  dplyr::group_by(
    Military,
    Always_Unlocked,
    Question,
    Motivator,
    Response
  ) %>%
  dplyr::summarise(
    Weighted_Count = sum(Weights, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  dplyr::group_by(
    Military,
    Always_Unlocked,
    Question,
    Motivator
  ) %>%
  dplyr::mutate(
    Agree_Weighted_Total =
      sum(Weighted_Count, na.rm = TRUE),
    Within_Agreement_Percent =
      100 * Weighted_Count / Agree_Weighted_Total
  ) %>%
  dplyr::ungroup() %>%
  dplyr::left_join(
    motivator_top_two %>%
      dplyr::select(
        Military,
        Always_Unlocked,
        Question,
        Percent
      ),
    by = c(
      "Military",
      "Always_Unlocked",
      "Question"
    )
  ) %>%
  dplyr::mutate(
    Segment_Percent =
      Percent * Within_Agreement_Percent / 100
  )
ggplot(
  motivator_stacked,
  aes(
    x = Segment_Percent,
    y = forcats::fct_reorder(Motivator, Percent),
    fill = Response
  )
) +
  geom_col() +
  facet_grid(
    Military ~ Always_Unlocked,
    scales = "free_y"
  ) +
  scale_x_continuous(
    limits = c(0, 100),
    labels = function(x) paste0(x, "%")
  ) +
  labs(
    x = "Survey-weighted agreement",
    y = NULL,
    fill = NULL,
    title = "Motivational endorsement by military experience",
    subtitle = "Bars show somewhat agree and strongly agree"
  ) +
  theme_classic(base_size = 11) +
  theme(
    legend.position = "top",
    strip.background = element_blank(),
    strip.text = element_text(face = "bold")
  )

14 13. Descriptive results summary

secure_no <- secure_by_military %>%
  dplyr::filter(
    Military == "No Military Experience"
  ) %>%
  dplyr::pull(Weighted_Percent)

secure_military <- secure_by_military %>%
  dplyr::filter(
    Military == "Military Experience"
  ) %>%
  dplyr::pull(Weighted_Percent)

cat("## Main descriptive findings\n\n")

14.1 Main descriptive findings

cat(
  "The firearm-owner sample included ",
  nrow(owners),
  " respondents: ",
  sum(owners$Military == "No Military Experience"),
  " without military experience and ",
  sum(owners$Military == "Military Experience"),
  " with military experience. ",
  "Survey-weighted secure-container storage prevalence was ",
  sprintf("%.1f%%", secure_no),
  " among respondents without military experience and ",
  sprintf("%.1f%%", secure_military),
  " among respondents with military experience.\n\n",
  sep = ""
)

The firearm-owner sample included 336 respondents: 282 without military experience and 54 with military experience. Survey-weighted secure-container storage prevalence was 57.2% among respondents without military experience and 67.1% among respondents with military experience.

cat(
  "The tables and figures above reproduce the original publication's analyses ",
  "while displaying every estimate separately by military experience. ",
  "Particular attention should be given to differences in secure-container use, ",
  "carry-frequency patterns, alternative safety practices, willingness to adopt ",
  "unused practices, intention to keep a firearm unlocked, and the circumstances ",
  "that respondents identified as potential motivators for changing storage behavior.\n"
)

The tables and figures above reproduce the original publication’s analyses while displaying every estimate separately by military experience. Particular attention should be given to differences in secure-container use, carry-frequency patterns, alternative safety practices, willingness to adopt unused practices, intention to keep a firearm unlocked, and the circumstances that respondents identified as potential motivators for changing storage behavior.

15 14. Exploratory extensions for possible inclusion

The sections below go beyond the purely descriptive replication. They are separated from the main analyses so they can be omitted from a descriptive Brief Report or moved to supplemental material.

15.1 14.1 Adjusted military association with secure-container storage

adjusted_data <- owners %>%
  dplyr::filter(
    !is.na(Secure_Container),
    !is.na(Military),
    !is.na(age.cat),
    !is.na(ppgender),
    !is.na(ppeduc5),
    !is.na(ppmarit5),
    !is.na(Weights)
  ) %>%
  droplevels()

design_adjusted <- survey::svydesign(
  ids = ~1,
  weights = ~Weights,
  data = adjusted_data
)

adjusted_secure_model <- survey::svyglm(
  Secure_Container ~
    Military +
    age.cat +
    ppgender +
    ppeduc5 +
    ppmarit5,
  design = design_adjusted,
  family = quasibinomial()
)

adjusted_secure_result <- broom::tidy(
  adjusted_secure_model,
  exponentiate = TRUE,
  conf.int = TRUE
) %>%
  dplyr::filter(
    term == "MilitaryMilitary Experience"
  )

adjusted_secure_result

15.2 14.2 Does the carry-frequency pattern differ by military experience?

The full six-category interaction may be unstable because some military cells are small. The first model uses the original categories. The second collapses carry into never, occasional, and frequent categories as a sensitivity analysis.

carry_interaction_model <- survey::svyglm(
  Secure_Container ~
    Military * Carry +
    age.cat +
    ppgender +
    ppeduc5 +
    ppmarit5,
  design = design_carry,
  family = quasibinomial()
)

carry_interaction_test <- tryCatch(
  survey::regTermTest(
    carry_interaction_model,
    ~Military:Carry
  ),
  error = function(e) NULL
)

carry_interaction_test
## Wald test for Military:Carry
##  in svyglm(formula = Secure_Container ~ Military * Carry + age.cat + 
##     ppgender + ppeduc5 + ppmarit5, design = design_carry, family = quasibinomial())
## F =  37.98631  on  5  and  301  df: p= < 2.22e-16
carry_data3 <- carry_data %>%
  dplyr::mutate(
    Carry3 = dplyr::case_when(
      Carry == "Never" ~ "Never",
      Carry %in% c(
        "Almost never",
        "Annually",
        "Monthly"
      ) ~ "Occasional",
      Carry %in% c(
        "Weekly",
        "Daily"
      ) ~ "Frequent",
      TRUE ~ NA_character_
    ),

    Carry3 = factor(
      Carry3,
      levels = c(
        "Never",
        "Occasional",
        "Frequent"
      )
    )
  ) %>%
  dplyr::filter(!is.na(Carry3)) %>%
  droplevels()

design_carry3 <- survey::svydesign(
  ids = ~1,
  weights = ~Weights,
  data = carry_data3
)

carry3_model <- survey::svyglm(
  Secure_Container ~
    Military * Carry3 +
    age.cat +
    ppgender +
    ppeduc5 +
    ppmarit5,
  design = design_carry3,
  family = quasibinomial()
)

carry3_interaction_test <- tryCatch(
  survey::regTermTest(
    carry3_model,
    ~Military:Carry3
  ),
  error = function(e) NULL
)

carry3_interaction_test
## Wald test for Military:Carry3
##  in svyglm(formula = Secure_Container ~ Military * Carry3 + age.cat + 
##     ppgender + ppeduc5 + ppmarit5, design = design_carry3, family = quasibinomial())
## F =  0.3576757  on  2  and  307  df: p= 0.69959

15.3 14.3 Children in the household and military experience

children_data <- owners %>%
  dplyr::filter(
    !is.na(Secure_Container),
    !is.na(Military),
    !is.na(Children),
    !is.na(age.cat),
    !is.na(ppgender),
    !is.na(ppeduc5),
    !is.na(ppmarit5),
    !is.na(Weights)
  ) %>%
  droplevels()

design_children <- survey::svydesign(
  ids = ~1,
  weights = ~Weights,
  data = children_data
)

children_interaction_model <- survey::svyglm(
  Secure_Container ~
    Military * Children +
    age.cat +
    ppgender +
    ppeduc5 +
    ppmarit5,
  design = design_children,
  family = quasibinomial()
)

children_interaction_test <- survey::regTermTest(
  children_interaction_model,
  ~Military:Children
)

children_interaction_test
## Wald test for Military:Children
##  in svyglm(formula = Secure_Container ~ Military * Children + age.cat + 
##     ppgender + ppeduc5 + ppmarit5, design = design_children, 
##     family = quasibinomial())
## F =  1.565524  on  1  and  312  df: p= 0.2118

15.4 14.4 Method-specific military comparisons with FDR correction

These analyses are exploratory and are best placed in supplemental material. Extremely sparse methods are flagged.

fit_method_model <- function(j) {

  outcome <- paste0("Use_", j)

  model_data <- owners %>%
    dplyr::filter(
      !is.na(.data[[outcome]]),
      !is.na(Military),
      !is.na(age.cat),
      !is.na(ppgender),
      !is.na(ppeduc5),
      !is.na(ppmarit5),
      !is.na(Weights)
    ) %>%
    droplevels()

  military_events <- sum(
    model_data[[outcome]] == 1 &
      model_data$Military == "Military Experience",
    na.rm = TRUE
  )

  nonmilitary_events <- sum(
    model_data[[outcome]] == 1 &
      model_data$Military == "No Military Experience",
    na.rm = TRUE
  )

  design_method <- survey::svydesign(
    ids = ~1,
    weights = ~Weights,
    data = model_data
  )

  formula_method <- as.formula(
    paste(
      outcome,
      "~ Military + age.cat + ppgender + ppeduc5 + ppmarit5"
    )
  )

  fit <- tryCatch(
    survey::svyglm(
      formula_method,
      design = design_method,
      family = quasibinomial()
    ),
    error = function(e) NULL
  )

  if (is.null(fit)) {
    return(
      tibble::tibble(
        Method = unname(method_labels[outcome]),
        Military_events = military_events,
        Nonmilitary_events = nonmilitary_events,
        OR = NA_real_,
        CI_Lower = NA_real_,
        CI_Upper = NA_real_,
        p_value = NA_real_,
        Status = "Model failed"
      )
    )
  }

  result <- broom::tidy(
    fit,
    exponentiate = TRUE,
    conf.int = TRUE
  ) %>%
    dplyr::filter(
      term == "MilitaryMilitary Experience"
    )

  if (nrow(result) == 0) {
    return(
      tibble::tibble(
        Method = unname(method_labels[outcome]),
        Military_events = military_events,
        Nonmilitary_events = nonmilitary_events,
        OR = NA_real_,
        CI_Lower = NA_real_,
        CI_Upper = NA_real_,
        p_value = NA_real_,
        Status = "Military estimate unavailable"
      )
    )
  }

  result %>%
    dplyr::transmute(
      Method = unname(method_labels[outcome]),
      Military_events = military_events,
      Nonmilitary_events = nonmilitary_events,
      OR = estimate,
      CI_Lower = conf.low,
      CI_Upper = conf.high,
      p_value = p.value,
      Status = dplyr::if_else(
        military_events < 5 |
          nonmilitary_events < 5,
        "Sparse; interpret cautiously",
        "Estimated"
      )
    )
}

method_models <- purrr::map_dfr(
  1:10,
  fit_method_model
) %>%
  dplyr::mutate(
    p_FDR = p.adjust(
      p_value,
      method = "BH"
    )
  )

method_models %>%
  dplyr::mutate(
    `Adjusted OR (95% CI)` = dplyr::if_else(
      is.na(OR),
      "Not estimable",
      sprintf(
        "%.2f (%.2f–%.2f)",
        OR,
        CI_Lower,
        CI_Upper
      )
    ),
    `p value` = format_p(p_value),
    `FDR p` = format_p(p_FDR)
  ) %>%
  dplyr::select(
    Method,
    Military_events,
    Nonmilitary_events,
    `Adjusted OR (95% CI)`,
    `p value`,
    `FDR p`,
    Status
  ) %>%
  gt() %>%
  tab_header(
    title = "Exploratory adjusted comparisons of individual storage methods"
  ) %>%
  tab_source_note(
    "Models adjust for age, gender, education, and marital status. Benjamini–Hochberg correction is applied across storage-method outcomes."
  )
Exploratory adjusted comparisons of individual storage methods
Method Military_events Nonmilitary_events Adjusted OR (95% CI) p value FDR p Status
Keep all firearms unloaded 24 138 0.75 (0.36–1.58) 0.452 0.646 Estimated
Use a secure container 35 158 2.09 (0.95–4.62) 0.068 0.171 Estimated
Use a cable or trigger lock 16 48 1.59 (0.72–3.53) 0.255 0.425 Estimated
Use an access alarm or notification 4 4 9.28 (1.03–83.50) 0.047 0.156 Sparse; interpret cautiously
Use a firearm sensor 1 1 0.00 (0.00–0.00) <.001 <.001 Sparse; interpret cautiously
Lock ammunition separately 20 100 1.18 (0.57–2.41) 0.657 0.657 Estimated
Disassemble firearms 5 16 5.74 (1.29–25.61) 0.022 0.111 Estimated
Entrust keys or parts to another 2 15 0.69 (0.17–2.78) 0.600 0.657 Sparse; interpret cautiously
Use another safety measure 1 8 0.50 (0.05–5.31) 0.568 0.657 Sparse; interpret cautiously
None of the above 6 47 0.51 (0.17–1.52) 0.227 0.425 Estimated
Models adjust for age, gender, education, and marital status. Benjamini–Hochberg correction is applied across storage-method outcomes.

16 15. Save key outputs

The document prints all results directly. This section also saves the central tables and figures for manuscript preparation.

output_directory <- "military_descriptive_outputs"

dir.create(
  output_directory,
  showWarnings = FALSE
)

write.csv(
  secure_by_military,
  file.path(
    output_directory,
    "Secure_Container_by_Military.csv"
  ),
  row.names = FALSE
)

write.csv(
  demographic_secure_results,
  file.path(
    output_directory,
    "Demographic_Secure_Storage_by_Military.csv"
  ),
  row.names = FALSE
)

write.csv(
  carry_estimates,
  file.path(
    output_directory,
    "Carry_Frequency_by_Military.csv"
  ),
  row.names = FALSE
)

write.csv(
  table2_combined,
  file.path(
    output_directory,
    "Current_and_Willingness_by_Military.csv"
  ),
  row.names = FALSE
)

write.csv(
  motivator_top_two,
  file.path(
    output_directory,
    "Motivators_by_Military_and_Unlocked_Attitude.csv"
  ),
  row.names = FALSE
)

write.csv(
  method_models,
  file.path(
    output_directory,
    "Exploratory_Method_Models.csv"
  ),
  row.names = FALSE
)

18 17. Reproducibility information

sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.7.4
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] knitr_1.51      patchwork_1.3.2 gt_1.3.0        forcats_1.0.1  
##  [5] scales_1.4.0    ggplot2_4.0.2   broom_1.0.12    survey_4.5     
##  [9] survival_3.8-3  Matrix_1.7-4    purrr_1.2.2     tidyr_1.3.2    
## [13] dplyr_1.2.1     foreign_0.8-90 
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     xml2_1.5.2         lattice_0.22-7    
##  [5] digest_0.6.39      magrittr_2.0.4     evaluate_1.0.5     RColorBrewer_1.1-3
##  [9] fastmap_1.2.0      jsonlite_2.0.0     backports_1.5.0    DBI_1.3.0         
## [13] jquerylib_0.1.4    cli_3.6.5          mitools_2.4        rlang_1.3.0       
## [17] splines_4.5.2      withr_3.0.2        cachem_1.1.0       yaml_2.3.12       
## [21] otel_0.2.0         tools_4.5.2        vctrs_0.7.1        R6_2.6.1          
## [25] lifecycle_1.0.5    fs_1.6.7           pkgconfig_2.0.3    pillar_1.11.1     
## [29] bslib_0.10.0       gtable_0.3.6       glue_1.8.1         Rcpp_1.1.2        
## [33] xfun_0.56          tibble_3.3.1       tidyselect_1.2.1   rstudioapi_0.18.0 
## [37] farver_2.1.2       htmltools_0.5.9    labeling_0.4.3     rmarkdown_2.30    
## [41] compiler_4.5.2     S7_0.2.1