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:
- secure-container storage prevalence;
- secure-container storage across demographic and household
characteristics;
- secure-container storage across firearm-carry frequency;
- current use of and willingness to adopt additional firearm-safety
practices;
- 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.
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)
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 = ", ")
)
}
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)
}
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
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
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)
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. |
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
)
)
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
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
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")
)
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")
)
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")
)
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")
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.
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.
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
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
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
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. |
14.5 Method-specific military comparisons of willingness with FDR
correction
fit_willing_model <- function(j) {
use_outcome <- paste0("Use_", j)
willing_outcome <- paste0("Willing_", j)
# Willingness was only asked of respondents who were NOT
# currently using the corresponding method.
model_data <- owners %>%
dplyr::filter(
.data[[use_outcome]] == 0,
!is.na(.data[[willing_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[[willing_outcome]] == 1 &
model_data$Military == "Military Experience",
na.rm = TRUE
)
nonmilitary_events <- sum(
model_data[[willing_outcome]] == 1 &
model_data$Military == "No Military Experience",
na.rm = TRUE
)
design_willing <- survey::svydesign(
ids = ~1,
weights = ~Weights,
data = model_data
)
formula_willing <- as.formula(
paste(
willing_outcome,
"~ Military + age.cat + ppgender + ppeduc5 + ppmarit5"
)
)
fit <- tryCatch(
survey::svyglm(
formula_willing,
design = design_willing,
family = quasibinomial()
),
error = function(e) NULL
)
if (is.null(fit)) {
return(
tibble::tibble(
Method = unname(willing_labels[willing_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(willing_labels[willing_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(willing_labels[willing_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"
)
)
}
willing_models <- purrr::map_dfr(
1:9,
fit_willing_model
) %>%
dplyr::mutate(
p_FDR = p.adjust(
p_value,
method = "BH"
)
)
willing_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 = "Adjusted comparisons of willingness to adopt firearm-safety methods"
) %>%
tab_source_note(
"Models include respondents not currently using each method and adjust for age, gender, education, and marital status. Benjamini–Hochberg correction is applied across outcomes."
)
| Adjusted comparisons of willingness to adopt firearm-safety methods |
| Method |
Military_events |
Nonmilitary_events |
Adjusted OR (95% CI) |
p value |
FDR p |
Status |
| Keep all firearms unloaded |
6 |
33 |
0.50 (0.15–1.65) |
0.252 |
0.713 |
Estimated |
| Use a secure container |
5 |
60 |
0.49 (0.12–2.08) |
0.330 |
0.713 |
Estimated |
| Use a cable or trigger lock |
8 |
56 |
1.00 (0.39–2.60) |
0.997 |
0.997 |
Estimated |
| Use an access alarm or notification |
5 |
54 |
0.66 (0.21–2.09) |
0.481 |
0.713 |
Estimated |
| Use a firearm sensor |
2 |
25 |
0.62 (0.12–3.08) |
0.554 |
0.713 |
Sparse; interpret cautiously |
| Lock ammunition separately |
5 |
43 |
0.42 (0.14–1.28) |
0.127 |
0.713 |
Estimated |
| Disassemble firearms |
4 |
17 |
1.69 (0.45–6.33) |
0.439 |
0.713 |
Sparse; interpret cautiously |
| Entrust keys or parts to another |
3 |
13 |
1.33 (0.32–5.53) |
0.695 |
0.782 |
Sparse; interpret cautiously |
| Use another safety measure |
1 |
11 |
0.49 (0.06–4.16) |
0.516 |
0.713 |
Sparse; interpret cautiously |
| Models include respondents not currently using each method and adjust for age, gender, education, and marital status. Benjamini–Hochberg correction is applied across outcomes. |
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
)
## ============================================================
## Adjusted military comparisons of motivations for changing
## firearm-storage behavior
##
## Outcome:
## 1 = Somewhat agree OR Strongly agree
## 0 = Somewhat disagree OR Strongly disagree
##
## Predictor of interest:
## Military experience
##
## Covariates:
## Age, gender, education, marital status
## ============================================================
# Labels for the six motivation items
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"
)
# ------------------------------------------------------------
# 1. Build long dataset
# ------------------------------------------------------------
motivator_long_adjusted <- owners %>%
dplyr::select(
Respondent_ID,
Military,
Weights,
age.cat,
ppgender,
ppeduc5,
ppmarit5,
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(
# Use your existing function to standardize the
# four agreement response categories
Response = agreement_four(Response),
Response = factor(
Response,
levels = c(
"Strongly disagree",
"Somewhat disagree",
"Somewhat agree",
"Strongly agree"
)
),
Motivator = unname(motivator_labels[Question]),
# Binary outcome requested:
# Somewhat/Strongly Agree vs. Disagree
Top_Two = dplyr::case_when(
Response %in% c(
"Somewhat agree",
"Strongly agree"
) ~ 1L,
Response %in% c(
"Somewhat disagree",
"Strongly disagree"
) ~ 0L,
TRUE ~ NA_integer_
)
) %>%
dplyr::filter(
!is.na(Top_Two),
!is.na(Military),
!is.na(Weights)
)
# Make sure non-military is the reference group
motivator_long_adjusted <- motivator_long_adjusted %>%
dplyr::mutate(
Military = factor(
Military,
levels = c(
"No Military Experience",
"Military Experience"
)
)
)
# ------------------------------------------------------------
# 2. Function to fit one adjusted model per motivation
# ------------------------------------------------------------
fit_motivator_model <- function(question_name) {
model_data <- motivator_long_adjusted %>%
dplyr::filter(
Question == question_name,
!is.na(age.cat),
!is.na(ppgender),
!is.na(ppeduc5),
!is.na(ppmarit5)
) %>%
droplevels()
# Survey-weighted design
design_temp <- survey::svydesign(
ids = ~1,
weights = ~Weights,
data = model_data
)
# Quasibinomial regression
fit <- survey::svyglm(
Top_Two ~
Military +
age.cat +
ppgender +
ppeduc5 +
ppmarit5,
design = design_temp,
family = quasibinomial()
)
# Pull military effect only
broom::tidy(
fit,
conf.int = TRUE,
exponentiate = TRUE
) %>%
dplyr::filter(
term == "MilitaryMilitary Experience"
) %>%
dplyr::transmute(
Question = question_name,
Motivator = unname(motivator_labels[question_name]),
OR = estimate,
CI_Lower = conf.low,
CI_Upper = conf.high,
p_value = p.value,
N = nrow(model_data)
)
}
# ------------------------------------------------------------
# 3. Run models for all six motivations
# ------------------------------------------------------------
motivator_models <- purrr::map_dfr(
paste0("BMW4_", 1:6),
fit_motivator_model
) %>%
dplyr::mutate(
# Benjamini-Hochberg FDR correction across the six tests
p_FDR = p.adjust(
p_value,
method = "BH"
)
)
motivator_models
## ------------------------------------------------------------
## 4. Weighted top-two percentages by military status
## ------------------------------------------------------------
design_motivator_adjusted <- survey::svydesign(
ids = ~1,
weights = ~Weights,
data = motivator_long_adjusted
)
motivator_percentages <- survey::svyby(
~Top_Two,
~Military + Question + Motivator,
design = design_motivator_adjusted,
FUN = survey::svymean,
vartype = c("ci"),
na.rm = TRUE,
keep.names = FALSE
) %>%
as.data.frame() %>%
dplyr::mutate(
Percent = 100 * Top_Two,
CI_Lower_Percent = 100 * ci_l,
CI_Upper_Percent = 100 * ci_u
)
motivator_percentages
## ------------------------------------------------------------
## 5. Combine weighted percentages and adjusted odds ratios
## ------------------------------------------------------------
motivator_percent_wide <- motivator_percentages %>%
dplyr::select(
Question,
Motivator,
Military,
Percent
) %>%
tidyr::pivot_wider(
names_from = Military,
values_from = Percent
) %>%
dplyr::rename(
NonMilitary_Percent = `No Military Experience`,
Military_Percent = `Military Experience`
)
motivator_final_table <- motivator_percent_wide %>%
dplyr::left_join(
motivator_models,
by = c("Question", "Motivator")
) %>%
dplyr::mutate(
NonMilitary_Percent = sprintf(
"%.1f%%",
NonMilitary_Percent
),
Military_Percent = sprintf(
"%.1f%%",
Military_Percent
),
Adjusted_OR = sprintf(
"%.2f (%.2f–%.2f)",
OR,
CI_Lower,
CI_Upper
),
p_value_display = dplyr::case_when(
p_value < .001 ~ "<.001",
TRUE ~ sprintf("%.3f", p_value)
),
p_FDR_display = dplyr::case_when(
p_FDR < .001 ~ "<.001",
TRUE ~ sprintf("%.3f", p_FDR)
)
)
motivator_final_table
motivator_final_table %>%
dplyr::select(
Motivator,
NonMilitary_Percent,
Military_Percent,
Adjusted_OR,
p_value_display,
p_FDR_display
) %>%
dplyr::rename(
`Motivation for changing storage` = Motivator,
`No military experience` = NonMilitary_Percent,
`Military experience` = Military_Percent,
`Adjusted OR (95% CI)` = Adjusted_OR,
`p` = p_value_display,
`FDR-adjusted p` = p_FDR_display
) %>%
gt::gt() %>%
gt::tab_header(
title = "Motivations for Changing Firearm-Storage Practices by Military Experience"
) %>%
gt::tab_source_note(
gt::md(
"Percentages represent the survey-weighted proportion who **somewhat or strongly agreed**. Adjusted odds ratios compare respondents with military experience to those without military experience and adjust for age, gender, education, and marital status."
)
)
| Motivations for Changing Firearm-Storage Practices by Military Experience |
| Motivation for changing storage |
No military experience |
Military experience |
Adjusted OR (95% CI) |
p |
FDR-adjusted p |
| A close friend or family member asked |
79.9% |
54.2% |
0.49 (0.24–1.00) |
0.051 |
0.306 |
| Household mental or physical health concerns |
95.6% |
92.0% |
1.09 (0.28–4.17) |
0.900 |
0.922 |
| If I had children at home |
93.7% |
88.2% |
1.19 (0.39–3.65) |
0.758 |
0.922 |
| Prevent accidental injury |
88.1% |
70.4% |
0.71 (0.31–1.66) |
0.430 |
0.922 |
| Prevent suicide for me or others |
85.9% |
80.1% |
1.04 (0.43–2.52) |
0.922 |
0.922 |
| Prevent use, theft, or damage |
88.7% |
77.9% |
0.86 (0.36–2.09) |
0.740 |
0.922 |
| Percentages represent the survey-weighted proportion who somewhat or strongly agreed. Adjusted odds ratios compare respondents with military experience to those without military experience and adjust for age, gender, education, and marital status. |
Formal difference in keeping one unlocked
# ============================================================
# Military experience and endorsement of keeping at least
# one firearm unlocked
# Standalone version starting from `data`
# ============================================================
library(dplyr)
library(survey)
library(broom)
library(emmeans)
# ------------------------------------------------------------
# 1. Build analysis dataset directly from data
# ------------------------------------------------------------
data <- foreign::read.spss(
"KP_OMNI_2405_BMW_Client_File_03042024.sav",
to.data.frame = TRUE,
use.value.labels = TRUE
)
data <- as.data.frame(data)
unlocked_data <- data %>%
dplyr::filter(
BMW1 == "Yes",
!is.na(Status),
!is.na(BMW4_7),
!is.na(Weights)
) %>%
dplyr::mutate(
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"
)
),
Unlocked_Binary = dplyr::case_when(
BMW4_7 %in% c(
"Strongly agree",
"Somewhat agree"
) ~ 1L,
BMW4_7 %in% c(
"Somewhat disagree",
"Strongly disagree"
) ~ 0L,
TRUE ~ NA_integer_
),
age.cat = cut(
ppage,
breaks = c(18, 25, 40, 60, Inf),
labels = c("18-25", "26-40", "41-60", ">60"),
right = FALSE
),
ppgender = factor(ppgender),
ppeduc5 = factor(ppeduc5),
ppmarit5 = factor(ppmarit5)
) %>%
dplyr::filter(
!is.na(Unlocked_Binary)
) %>%
droplevels()
# ------------------------------------------------------------
# 2. Unweighted counts
# ------------------------------------------------------------
table(
unlocked_data$Military,
unlocked_data$Unlocked_Binary
)
##
## 0 1
## No Military Experience 133 146
## Military Experience 22 31
unlocked_data %>%
dplyr::group_by(Military) %>%
dplyr::summarise(
n = dplyr::n(),
n_unlocked = sum(Unlocked_Binary == 1),
percent_unweighted = 100 * mean(Unlocked_Binary == 1),
.groups = "drop"
)
# ------------------------------------------------------------
# 3. Survey design
# ------------------------------------------------------------
design_unlocked <- survey::svydesign(
ids = ~1,
weights = ~Weights,
data = unlocked_data
)
# ------------------------------------------------------------
# 4. Weighted prevalence by military experience
# ------------------------------------------------------------
unlocked_weighted <- survey::svyby(
~Unlocked_Binary,
~Military,
design = design_unlocked,
FUN = survey::svymean,
vartype = c("se", "ci"),
na.rm = TRUE,
keep.names = FALSE
) %>%
as.data.frame() %>%
dplyr::mutate(
Percent = 100 * Unlocked_Binary,
CI_Lower = 100 * ci_l,
CI_Upper = 100 * ci_u
)
unlocked_weighted
# ------------------------------------------------------------
# 5. Unadjusted survey-weighted comparison
# ------------------------------------------------------------
unlocked_data <- unlocked_data %>%
dplyr::mutate(
Unlocked_Factor = factor(
Unlocked_Binary,
levels = c(0, 1),
labels = c(
"Does not endorse keeping one unlocked",
"Endorses keeping one unlocked"
)
)
)
design_unlocked <- survey::svydesign(
ids = ~1,
weights = ~Weights,
data = unlocked_data
)
unlocked_test <- survey::svychisq(
~Military + Unlocked_Factor,
design = design_unlocked,
statistic = "F"
)
unlocked_test
##
## Pearson's X^2: Rao & Scott adjustment
##
## data: survey::svychisq(~Military + Unlocked_Factor, design = design_unlocked, statistic = "F")
## F = 0.33309, ndf = 1, ddf = 331, p-value = 0.5642
# ------------------------------------------------------------
# 6. Adjusted model
# ------------------------------------------------------------
unlocked_adjusted_data <- unlocked_data %>%
dplyr::filter(
!is.na(age.cat),
!is.na(ppgender),
!is.na(ppeduc5),
!is.na(ppmarit5)
) %>%
droplevels()
design_unlocked_adjusted <- survey::svydesign(
ids = ~1,
weights = ~Weights,
data = unlocked_adjusted_data
)
unlocked_model <- survey::svyglm(
Unlocked_Binary ~
Military +
age.cat +
ppgender +
ppeduc5 +
ppmarit5,
design = design_unlocked_adjusted,
family = quasibinomial()
)
summary(unlocked_model)
##
## Call:
## svyglm(formula = Unlocked_Binary ~ Military + age.cat + ppgender +
## ppeduc5 + ppmarit5, design = design_unlocked_adjusted, family = quasibinomial())
##
## Survey design:
## survey::svydesign(ids = ~1, weights = ~Weights, data = unlocked_adjusted_data)
##
## Coefficients:
## Estimate
## (Intercept) 0.71465
## MilitaryMilitary Experience -0.08155
## age.cat26-40 -0.58235
## age.cat41-60 -0.11256
## age.cat>60 -0.12118
## ppgenderFemale -0.67474
## ppeduc5High school graduate (high school diploma or the equivalent GED) 0.27702
## ppeduc5Some college or Associate degree 0.05477
## ppeduc5Bachelor’s degree -0.53128
## ppeduc5Master’s degree or above -0.97311
## ppmarit5Widowed 0.13774
## ppmarit5Divorced -0.04964
## ppmarit5Separated -0.25173
## ppmarit5Never married 0.15124
## Std. Error
## (Intercept) 0.80292
## MilitaryMilitary Experience 0.38145
## age.cat26-40 0.54340
## age.cat41-60 0.56839
## age.cat>60 0.60826
## ppgenderFemale 0.26578
## ppeduc5High school graduate (high school diploma or the equivalent GED) 0.58133
## ppeduc5Some college or Associate degree 0.58041
## ppeduc5Bachelor’s degree 0.59566
## ppeduc5Master’s degree or above 0.62287
## ppmarit5Widowed 0.57892
## ppmarit5Divorced 0.40648
## ppmarit5Separated 1.26523
## ppmarit5Never married 0.38199
## t value
## (Intercept) 0.890
## MilitaryMilitary Experience -0.214
## age.cat26-40 -1.072
## age.cat41-60 -0.198
## age.cat>60 -0.199
## ppgenderFemale -2.539
## ppeduc5High school graduate (high school diploma or the equivalent GED) 0.477
## ppeduc5Some college or Associate degree 0.094
## ppeduc5Bachelor’s degree -0.892
## ppeduc5Master’s degree or above -1.562
## ppmarit5Widowed 0.238
## ppmarit5Divorced -0.122
## ppmarit5Separated -0.199
## ppmarit5Never married 0.396
## Pr(>|t|)
## (Intercept) 0.3741
## MilitaryMilitary Experience 0.8309
## age.cat26-40 0.2847
## age.cat41-60 0.8431
## age.cat>60 0.8422
## ppgenderFemale 0.0116
## ppeduc5High school graduate (high school diploma or the equivalent GED) 0.6340
## ppeduc5Some college or Associate degree 0.9249
## ppeduc5Bachelor’s degree 0.3731
## ppeduc5Master’s degree or above 0.1192
## ppmarit5Widowed 0.8121
## ppmarit5Divorced 0.9029
## ppmarit5Separated 0.8424
## ppmarit5Never married 0.6924
##
## (Intercept)
## MilitaryMilitary Experience
## age.cat26-40
## age.cat41-60
## age.cat>60
## ppgenderFemale *
## ppeduc5High school graduate (high school diploma or the equivalent GED)
## ppeduc5Some college or Associate degree
## ppeduc5Bachelor’s degree
## ppeduc5Master’s degree or above
## ppmarit5Widowed
## ppmarit5Divorced
## ppmarit5Separated
## ppmarit5Never married
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 1.000524)
##
## Number of Fisher Scoring iterations: 4
unlocked_OR <- broom::tidy(
unlocked_model,
exponentiate = TRUE,
conf.int = TRUE
) %>%
dplyr::filter(
term == "MilitaryMilitary Experience"
)
unlocked_OR
# ------------------------------------------------------------
# 7. Adjusted probabilities
# ------------------------------------------------------------
unlocked_emm <- emmeans::emmeans(
unlocked_model,
~Military,
type = "response"
)
summary(unlocked_emm)
# ============================================================
# Friend/family motivation adjusted for unlocked attitude
# ============================================================
friend_family_data <- data %>%
dplyr::filter(
BMW1 == "Yes",
!is.na(Status),
!is.na(BMW4_3),
!is.na(BMW4_7),
!is.na(Weights)
) %>%
dplyr::mutate(
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"
)
),
Friend_Family = dplyr::case_when(
BMW4_3 %in% c(
"Strongly agree",
"Somewhat agree"
) ~ 1L,
BMW4_3 %in% c(
"Somewhat disagree",
"Strongly disagree"
) ~ 0L,
TRUE ~ NA_integer_
),
Unlocked_Binary = dplyr::case_when(
BMW4_7 %in% c(
"Strongly agree",
"Somewhat agree"
) ~ 1L,
BMW4_7 %in% c(
"Somewhat disagree",
"Strongly disagree"
) ~ 0L,
TRUE ~ NA_integer_
),
age.cat = cut(
ppage,
breaks = c(18, 25, 40, 60, Inf),
labels = c("18-25", "26-40", "41-60", ">60"),
right = FALSE
),
ppgender = factor(ppgender),
ppeduc5 = factor(ppeduc5),
ppmarit5 = factor(ppmarit5)
) %>%
dplyr::filter(
!is.na(Friend_Family),
!is.na(Unlocked_Binary),
!is.na(age.cat),
!is.na(ppgender),
!is.na(ppeduc5),
!is.na(ppmarit5)
) %>%
droplevels()
design_friend_family <- survey::svydesign(
ids = ~1,
weights = ~Weights,
data = friend_family_data
)
# Main model: does military experience predict friend/family motivation,
# controlling for unlocked preference and demographics?
friend_family_model <- survey::svyglm(
Friend_Family ~
Military +
Unlocked_Binary +
age.cat +
ppgender +
ppeduc5 +
ppmarit5,
design = design_friend_family,
family = quasibinomial()
)
broom::tidy(
friend_family_model,
exponentiate = TRUE,
conf.int = TRUE
)
# Interaction model:
# does the military difference depend on unlocked preference?
friend_family_interaction <- survey::svyglm(
Friend_Family ~
Military * Unlocked_Binary +
age.cat +
ppgender +
ppeduc5 +
ppmarit5,
design = design_friend_family,
family = quasibinomial()
)
summary(friend_family_interaction)
##
## Call:
## svyglm(formula = Friend_Family ~ Military * Unlocked_Binary +
## age.cat + ppgender + ppeduc5 + ppmarit5, design = design_friend_family,
## family = quasibinomial())
##
## Survey design:
## survey::svydesign(ids = ~1, weights = ~Weights, data = friend_family_data)
##
## Coefficients:
## Estimate
## (Intercept) 1.16466
## MilitaryMilitary Experience -1.27336
## Unlocked_Binary -1.89101
## age.cat26-40 0.67360
## age.cat41-60 0.49948
## age.cat>60 0.86242
## ppgenderFemale 0.65930
## ppeduc5High school graduate (high school diploma or the equivalent GED) 0.54714
## ppeduc5Some college or Associate degree -0.09859
## ppeduc5Bachelor’s degree -0.17811
## ppeduc5Master’s degree or above 0.36716
## ppmarit5Widowed 0.49765
## ppmarit5Divorced 0.75374
## ppmarit5Separated 15.94217
## ppmarit5Never married 1.29121
## MilitaryMilitary Experience:Unlocked_Binary 0.56500
## Std. Error
## (Intercept) 1.36103
## MilitaryMilitary Experience 0.61866
## Unlocked_Binary 0.41977
## age.cat26-40 0.72303
## age.cat41-60 0.72602
## age.cat>60 0.76035
## ppgenderFemale 0.35979
## ppeduc5High school graduate (high school diploma or the equivalent GED) 1.01220
## ppeduc5Some college or Associate degree 1.01029
## ppeduc5Bachelor’s degree 1.03940
## ppeduc5Master’s degree or above 1.08963
## ppmarit5Widowed 0.60791
## ppmarit5Divorced 0.64768
## ppmarit5Separated 0.83400
## ppmarit5Never married 0.50888
## MilitaryMilitary Experience:Unlocked_Binary 0.72003
## t value
## (Intercept) 0.856
## MilitaryMilitary Experience -2.058
## Unlocked_Binary -4.505
## age.cat26-40 0.932
## age.cat41-60 0.688
## age.cat>60 1.134
## ppgenderFemale 1.832
## ppeduc5High school graduate (high school diploma or the equivalent GED) 0.541
## ppeduc5Some college or Associate degree -0.098
## ppeduc5Bachelor’s degree -0.171
## ppeduc5Master’s degree or above 0.337
## ppmarit5Widowed 0.819
## ppmarit5Divorced 1.164
## ppmarit5Separated 19.115
## ppmarit5Never married 2.537
## MilitaryMilitary Experience:Unlocked_Binary 0.785
## Pr(>|t|)
## (Intercept) 0.3928
## MilitaryMilitary Experience 0.0404
## Unlocked_Binary 9.39e-06
## age.cat26-40 0.3522
## age.cat41-60 0.4920
## age.cat>60 0.2576
## ppgenderFemale 0.0678
## ppeduc5High school graduate (high school diploma or the equivalent GED) 0.5892
## ppeduc5Some college or Associate degree 0.9223
## ppeduc5Bachelor’s degree 0.8641
## ppeduc5Master’s degree or above 0.7364
## ppmarit5Widowed 0.4136
## ppmarit5Divorced 0.2454
## ppmarit5Separated < 2e-16
## ppmarit5Never married 0.0117
## MilitaryMilitary Experience:Unlocked_Binary 0.4332
##
## (Intercept)
## MilitaryMilitary Experience *
## Unlocked_Binary ***
## age.cat26-40
## age.cat41-60
## age.cat>60
## ppgenderFemale .
## ppeduc5High school graduate (high school diploma or the equivalent GED)
## ppeduc5Some college or Associate degree
## ppeduc5Bachelor’s degree
## ppeduc5Master’s degree or above
## ppmarit5Widowed
## ppmarit5Divorced
## ppmarit5Separated ***
## ppmarit5Never married *
## MilitaryMilitary Experience:Unlocked_Binary
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.9568434)
##
## Number of Fisher Scoring iterations: 15
survey::regTermTest(
friend_family_interaction,
~Military:Unlocked_Binary
)
## Wald test for Military:Unlocked_Binary
## in svyglm(formula = Friend_Family ~ Military * Unlocked_Binary +
## age.cat + ppgender + ppeduc5 + ppmarit5, design = design_friend_family,
## family = quasibinomial())
## F = 0.6157466 on 1 and 313 df: p= 0.43323
# Adjusted probabilities for the four groups
friend_family_emmeans <- emmeans::emmeans(
friend_family_interaction,
~Military | Unlocked_Binary,
type = "response"
)
summary(friend_family_emmeans)
pairs(friend_family_emmeans)
## Unlocked_Binary = 0:
## contrast odds.ratio SE df null z.ratio
## No Military Experience / Military Experience 3.57 2.21 Inf 1 2.058
## p.value
## 0.0396
##
## Unlocked_Binary = 1:
## contrast odds.ratio SE df null z.ratio
## No Military Experience / Military Experience 2.03 0.93 Inf 1 1.547
## p.value
## 0.1219
##
## Results are averaged over the levels of: age.cat, ppgender, ppeduc5, ppmarit5
## Tests are performed on the log odds ratio scale
Analysis for willingness
# ============================================================
# THREE-LEVEL STORAGE STATUS:
# CURRENTLY USES / WILLING TO USE / NOT WILLING TO USE
# ============================================================
library(dplyr)
library(tidyr)
library(purrr)
library(survey)
library(broom)
library(gt)
options(survey.lonely.psu = "adjust")
# ------------------------------------------------------------
# 1. Start from firearm owners
# ------------------------------------------------------------
storage3_data <- data %>%
dplyr::filter(
BMW1 == "Yes",
!is.na(Status),
!is.na(Weights)
) %>%
dplyr::mutate(
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"
)
),
age.cat = cut(
ppage,
breaks = c(18, 25, 40, 60, Inf),
labels = c("18-25", "26-40", "41-60", ">60"),
right = FALSE
),
ppgender = factor(ppgender),
ppeduc5 = factor(ppeduc5),
ppmarit5 = factor(ppmarit5)
)
# ------------------------------------------------------------
# 2. Helper for BMW5 checkbox fields
#
# Current use:
# selected = 1
# not selected = 0
# skipped/missing = NA
# ------------------------------------------------------------
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
)
}
# ------------------------------------------------------------
# 3. Helper for BMW6 willingness fields
#
# IMPORTANT:
# This assumes BMW6 is also stored as checkbox-style selected /
# not-selected variables.
#
# A missing BMW6 value is NOT automatically treated as unwilling.
# If they CURRENTLY USE the method, they are classified as
# "Currently uses" before BMW6 is considered.
# ------------------------------------------------------------
willing_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
)
}
# ------------------------------------------------------------
# 4. Labels
# ------------------------------------------------------------
storage_labels <- c(
"1" = "Keep all firearms unloaded",
"2" = "Use a secure container",
"3" = "Use a cable or trigger lock",
"4" = "Use an access alarm or notification",
"5" = "Use a firearm sensor",
"6" = "Lock ammunition separately",
"7" = "Disassemble firearms",
"8" = "Entrust keys or parts to another",
"9" = "Use another safety measure",
"10" = "None of the above"
)
# ------------------------------------------------------------
# 5. Create three-level outcome for each method
# ------------------------------------------------------------
for (j in 1:10) {
use_var <- paste0("BMW5_", j)
willing_var <- paste0("BMW6_", j)
use_binary <- checkbox_binary(
storage3_data[[use_var]]
)
willing_binary_j <- willing_binary(
storage3_data[[willing_var]]
)
new_var <- paste0("Storage_Status_", j)
storage3_data[[new_var]] <- dplyr::case_when(
# Current user always goes here,
# regardless of BMW6 being skipped/not shown
use_binary == 1 ~ "Currently uses",
# Not a current user, but selected willingness option
use_binary == 0 &
willing_binary_j == 1 ~
"Would be willing to use",
# Not a current user and did not select willingness option
use_binary == 0 &
willing_binary_j == 0 ~
"Would not be willing to use",
# Anything genuinely unresolved stays missing
TRUE ~ NA_character_
)
storage3_data[[new_var]] <- factor(
storage3_data[[new_var]],
levels = c(
"Currently uses",
"Would be willing to use",
"Would not be willing to use"
)
)
}
# ------------------------------------------------------------
# 6. DIAGNOSTIC CHECK:
# Make sure we are NOT accidentally deleting current users
# ------------------------------------------------------------
diagnostic_table <- purrr::map_dfr(
1:10,
function(j) {
status_var <- paste0("Storage_Status_", j)
storage3_data %>%
dplyr::group_by(Military) %>%
dplyr::summarise(
Method = storage_labels[as.character(j)],
Total_N = dplyr::n(),
Currently_Uses =
sum(
.data[[status_var]] == "Currently uses",
na.rm = TRUE
),
Willing =
sum(
.data[[status_var]] ==
"Would be willing to use",
na.rm = TRUE
),
Not_Willing =
sum(
.data[[status_var]] ==
"Would not be willing to use",
na.rm = TRUE
),
Missing =
sum(
is.na(.data[[status_var]])
),
.groups = "drop"
)
}
)
diagnostic_table
# ============================================================
# 7. Convert to long format
# ============================================================
storage3_long <- storage3_data %>%
dplyr::select(
Military,
Weights,
dplyr::starts_with("Storage_Status_")
) %>%
tidyr::pivot_longer(
cols = dplyr::starts_with("Storage_Status_"),
names_to = "Method_Number",
values_to = "Storage_Status"
) %>%
dplyr::mutate(
Method_Number =
gsub(
"Storage_Status_",
"",
Method_Number
),
Method =
unname(
storage_labels[Method_Number]
),
Storage_Status = factor(
Storage_Status,
levels = c(
"Currently uses",
"Would be willing to use",
"Would not be willing to use"
)
)
)
# ============================================================
# 8. Survey-weighted three-level distribution
# Robust version
# ============================================================
weighted_storage3 <- purrr::map_dfr(
1:10,
function(j) {
status_var <- paste0("Storage_Status_", j)
method_name <- storage_labels[
as.character(j)
]
# Keep everyone whose status for THIS method can be determined
temp <- storage3_data %>%
dplyr::filter(
!is.na(.data[[status_var]])
) %>%
droplevels()
# Survey design for this method
design_temp <- survey::svydesign(
ids = ~1,
weights = ~Weights,
data = temp
)
# Create survey-weighted table
tab <- survey::svytable(
as.formula(
paste0(
"~ Military + ",
status_var
)
),
design = design_temp
)
# Convert to within-military percentages
pct <- prop.table(
tab,
margin = 1
) * 100
# Convert to dataframe
out <- as.data.frame(pct)
# Rename by POSITION rather than assuming Var1 / Var2
names(out)[1:3] <- c(
"Military",
"Storage_Status",
"Percent"
)
# Add method information
out %>%
dplyr::mutate(
Method_Number = j,
Method = method_name
)
}
)
weighted_storage3
dplyr::glimpse(weighted_storage3)
## Rows: 60
## Columns: 5
## $ Military <fct> No Military Experience, Military Experience, No Militar…
## $ Storage_Status <fct> Currently uses, Currently uses, Would be willing to use…
## $ Percent <dbl> 50.4815382, 43.7832845, 12.1998976, 9.8951404, 37.31856…
## $ Method_Number <int> 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4…
## $ Method <chr> "Keep all firearms unloaded", "Keep all firearms unload…
weighted_storage3 %>%
dplyr::select(
Method,
Military,
Storage_Status,
Percent
) %>%
print()
## Method Military
## 1 Keep all firearms unloaded No Military Experience
## 2 Keep all firearms unloaded Military Experience
## 3 Keep all firearms unloaded No Military Experience
## 4 Keep all firearms unloaded Military Experience
## 5 Keep all firearms unloaded No Military Experience
## 6 Keep all firearms unloaded Military Experience
## 7 Use a secure container No Military Experience
## 8 Use a secure container Military Experience
## 9 Use a secure container No Military Experience
## 10 Use a secure container Military Experience
## 11 Use a secure container No Military Experience
## 12 Use a secure container Military Experience
## 13 Use a cable or trigger lock No Military Experience
## 14 Use a cable or trigger lock Military Experience
## 15 Use a cable or trigger lock No Military Experience
## 16 Use a cable or trigger lock Military Experience
## 17 Use a cable or trigger lock No Military Experience
## 18 Use a cable or trigger lock Military Experience
## 19 Use an access alarm or notification No Military Experience
## 20 Use an access alarm or notification Military Experience
## 21 Use an access alarm or notification No Military Experience
## 22 Use an access alarm or notification Military Experience
## 23 Use an access alarm or notification No Military Experience
## 24 Use an access alarm or notification Military Experience
## 25 Use a firearm sensor No Military Experience
## 26 Use a firearm sensor Military Experience
## 27 Use a firearm sensor No Military Experience
## 28 Use a firearm sensor Military Experience
## 29 Use a firearm sensor No Military Experience
## 30 Use a firearm sensor Military Experience
## 31 Lock ammunition separately No Military Experience
## 32 Lock ammunition separately Military Experience
## 33 Lock ammunition separately No Military Experience
## 34 Lock ammunition separately Military Experience
## 35 Lock ammunition separately No Military Experience
## 36 Lock ammunition separately Military Experience
## 37 Disassemble firearms No Military Experience
## 38 Disassemble firearms Military Experience
## 39 Disassemble firearms No Military Experience
## 40 Disassemble firearms Military Experience
## 41 Disassemble firearms No Military Experience
## 42 Disassemble firearms Military Experience
## 43 Entrust keys or parts to another No Military Experience
## 44 Entrust keys or parts to another Military Experience
## 45 Entrust keys or parts to another No Military Experience
## 46 Entrust keys or parts to another Military Experience
## 47 Entrust keys or parts to another No Military Experience
## 48 Entrust keys or parts to another Military Experience
## 49 Use another safety measure No Military Experience
## 50 Use another safety measure Military Experience
## 51 Use another safety measure No Military Experience
## 52 Use another safety measure Military Experience
## 53 Use another safety measure No Military Experience
## 54 Use another safety measure Military Experience
## 55 None of the above No Military Experience
## 56 None of the above Military Experience
## 57 None of the above No Military Experience
## 58 None of the above Military Experience
## 59 None of the above No Military Experience
## 60 None of the above Military Experience
## Storage_Status Percent
## 1 Currently uses 50.4815382
## 2 Currently uses 43.7832845
## 3 Would be willing to use 12.1998976
## 4 Would be willing to use 9.8951404
## 5 Would not be willing to use 37.3185642
## 6 Would not be willing to use 46.3215751
## 7 Currently uses 57.3461847
## 8 Currently uses 67.0979187
## 9 Would be willing to use 22.2942573
## 10 Would be willing to use 9.3710213
## 11 Would not be willing to use 20.3595580
## 12 Would not be willing to use 23.5310600
## 13 Currently uses 18.5255303
## 14 Currently uses 26.5453903
## 15 Would be willing to use 19.6427084
## 16 Would be willing to use 18.2682663
## 17 Would not be willing to use 61.8317613
## 18 Would not be willing to use 55.1863434
## 19 Currently uses 1.4174014
## 20 Currently uses 8.0292487
## 21 Would be willing to use 19.5974070
## 22 Would be willing to use 8.6915016
## 23 Would not be willing to use 78.9851915
## 24 Would not be willing to use 83.2792497
## 25 Currently uses 0.4510429
## 26 Currently uses 2.3106119
## 27 Would be willing to use 8.9210292
## 28 Would be willing to use 3.3076279
## 29 Would not be willing to use 90.6279279
## 30 Would not be willing to use 94.3817602
## 31 Currently uses 36.6049370
## 32 Currently uses 35.7117386
## 33 Would be willing to use 15.4533688
## 34 Would be willing to use 8.6244134
## 35 Would not be willing to use 47.9416943
## 36 Would not be willing to use 55.6638480
## 37 Currently uses 5.8774830
## 38 Currently uses 8.8429876
## 39 Would be willing to use 5.4937484
## 40 Would be willing to use 6.0472553
## 41 Would not be willing to use 88.6287686
## 42 Would not be willing to use 85.1097571
## 43 Currently uses 5.5493362
## 44 Currently uses 3.6273959
## 45 Would be willing to use 3.8905389
## 46 Would be willing to use 4.4005773
## 47 Would not be willing to use 90.5601249
## 48 Would not be willing to use 91.9720268
## 49 Currently uses 3.0077071
## 50 Currently uses 1.8443660
## 51 Would be willing to use 4.5921476
## 52 Would be willing to use 1.7508919
## 53 Would not be willing to use 92.4001453
## 54 Would not be willing to use 96.4047421
## 55 Currently uses 16.4831935
## 56 Currently uses 10.6497086
## 57 Would be willing to use 27.1706712
## 58 Would be willing to use 40.1375046
## 59 Would not be willing to use 56.3461353
## 60 Would not be willing to use 49.2127867
diagnostic_check <- diagnostic_table %>%
dplyr::mutate(
Accounted_For =
Currently_Uses +
Willing +
Not_Willing +
Missing,
All_Respondents_Accounted_For =
Accounted_For == Total_N
) %>%
dplyr::arrange(Method, Military)
#View(diagnostic_check)
storage3_summary <- weighted_storage3 %>%
dplyr::select(
Method,
Military,
Storage_Status,
Percent
) %>%
tidyr::pivot_wider(
names_from = Storage_Status,
values_from = Percent
) %>%
dplyr::arrange(Method, Military)
storage3_summary
storage3_differences <- weighted_storage3 %>%
dplyr::select(
Method,
Military,
Storage_Status,
Percent
) %>%
tidyr::pivot_wider(
names_from = Military,
values_from = Percent
) %>%
dplyr::mutate(
Difference =
`Military Experience` -
`No Military Experience`
) %>%
dplyr::arrange(
Storage_Status,
dplyr::desc(abs(Difference))
)
storage3_differences
# ============================================================
# ADJUSTED MODELS:
# POTENTIAL ADOPTION OF EACH STORAGE METHOD
#
# Outcome:
# 1 = Would be willing to use
# 0 = Currently uses OR would not be willing to use
#
# NOTE:
# Current users remain in the denominator.
# ============================================================
library(dplyr)
library(purrr)
library(survey)
library(broom)
# ------------------------------------------------------------
# 1. Create potential-adopter variables
# ------------------------------------------------------------
for (j in 1:9) {
status_var <- paste0("Storage_Status_", j)
new_var <- paste0("Potential_Adopter_", j)
storage3_data[[new_var]] <- dplyr::case_when(
storage3_data[[status_var]] ==
"Would be willing to use" ~ 1L,
storage3_data[[status_var]] %in%
c(
"Currently uses",
"Would not be willing to use"
) ~ 0L,
TRUE ~ NA_integer_
)
}
# ------------------------------------------------------------
# 2. Function to run adjusted survey-weighted logistic model
# ------------------------------------------------------------
fit_potential_adopter_model <- function(j) {
outcome <- paste0("Potential_Adopter_", j)
model_data <- storage3_data %>%
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()
# Keep track of military/nonmilitary sample sizes
n_nonmilitary <- sum(
model_data$Military == "No Military Experience"
)
n_military <- sum(
model_data$Military == "Military Experience"
)
# Survey design
design_temp <- survey::svydesign(
ids = ~1,
weights = ~Weights,
data = model_data
)
# Adjusted model
fit <- survey::svyglm(
as.formula(
paste0(
outcome,
" ~ Military + age.cat + ppgender + ppeduc5 + ppmarit5"
)
),
design = design_temp,
family = quasibinomial()
)
# Extract military coefficient
result <- broom::tidy(
fit,
exponentiate = TRUE,
conf.int = TRUE
) %>%
dplyr::filter(
term == "MilitaryMilitary Experience"
)
tibble::tibble(
Method = storage_labels[as.character(j)],
Total_N = nrow(model_data),
Nonmilitary_N = n_nonmilitary,
Military_N = n_military,
OR = result$estimate,
CI_Lower = result$conf.low,
CI_Upper = result$conf.high,
p_value = result$p.value
)
}
# ------------------------------------------------------------
# 3. Run models for all actual safety practices
#
# 1:9 only
# BMW5/6_10 = "None of the above", so we do not treat that
# as a safety practice someone could adopt.
# ------------------------------------------------------------
potential_adopter_models <- purrr::map_dfr(
1:9,
fit_potential_adopter_model
)
# ------------------------------------------------------------
# 4. FDR correction across the 9 models
# ------------------------------------------------------------
potential_adopter_models <- potential_adopter_models %>%
dplyr::mutate(
p_FDR = p.adjust(
p_value,
method = "BH"
)
)
# ------------------------------------------------------------
# 5. Display results
# ------------------------------------------------------------
potential_adopter_models %>%
dplyr::mutate(
OR = round(OR, 2),
CI_Lower = round(CI_Lower, 2),
CI_Upper = round(CI_Upper, 2),
p_value = round(p_value, 4),
p_FDR = round(p_FDR, 4)
) %>%
dplyr::arrange(p_value)
# ============================================================
# STRATIFIED MODELS:
# Friend/family motivation ~ Military experience
#
# Run separately among:
# 1. Those who endorse keeping ≥1 firearm unlocked
# 2. Those who do not endorse keeping ≥1 firearm unlocked
# ============================================================
library(dplyr)
library(survey)
library(broom)
# ------------------------------------------------------------
# Function to fit model within each unlocked-storage group
# ------------------------------------------------------------
run_stratified_friend_model <- function(unlocked_value) {
model_data <- friend_family_data %>%
dplyr::filter(
Unlocked_Binary == unlocked_value,
!is.na(Friend_Family),
!is.na(Military),
!is.na(age.cat),
!is.na(ppgender),
!is.na(ppeduc5),
!is.na(ppmarit5),
!is.na(Weights)
) %>%
droplevels()
design_temp <- survey::svydesign(
ids = ~1,
weights = ~Weights,
data = model_data
)
model <- survey::svyglm(
Friend_Family ~
Military +
age.cat +
ppgender +
ppeduc5 +
ppmarit5,
design = design_temp,
family = quasibinomial()
)
result <- broom::tidy(
model,
exponentiate = TRUE,
conf.int = TRUE
) %>%
dplyr::filter(
term == "MilitaryMilitary Experience"
)
tibble::tibble(
Unlocked_Group = ifelse(
unlocked_value == 1,
"Endorses keeping ≥1 firearm unlocked",
"Does not endorse keeping ≥1 firearm unlocked"
),
Total_N = nrow(model_data),
Nonmilitary_N = sum(
model_data$Military == "No Military Experience"
),
Military_N = sum(
model_data$Military == "Military Experience"
),
aOR = result$estimate,
CI_Lower = result$conf.low,
CI_Upper = result$conf.high,
p_value = result$p.value
)
}
# ------------------------------------------------------------
# Run both models
# ------------------------------------------------------------
friend_stratified_results <- dplyr::bind_rows(
run_stratified_friend_model(0),
run_stratified_friend_model(1)
)
res=friend_stratified_results %>%
dplyr::mutate(
aOR = round(aOR, 2),
CI_Lower = round(CI_Lower, 2),
CI_Upper = round(CI_Upper, 2),
p_value = round(p_value, 4)
)
res
#16. Recommended publication structure
Based on the original paper, the descriptive military comparison
could be organized as follows:
Main manuscript
Table 1. Participant characteristics and
secure-container storage prevalence by military experience.
Figure 1. Secure-container storage by carry
frequency and military experience.
Table 2. Current use of and willingness to adopt
firearm-safety practices among military and nonmilitary firearm owners,
shown for the full sample and secure-container nonusers.
Figure 2. Motivators for changing firearm-storage
behavior, stratified by military experience and endorsement of always
keeping at least one firearm unlocked.
Supplemental material
Supplemental Table 1. Full demographic-specific
secure-container estimates.
Supplemental Table 2. Exploratory adjusted
individual-method models with FDR correction.
Supplemental Figure 1. Collapsed carry-frequency
comparison.
Supplemental Figure 2. Children-in-the-household
interaction.