Background

Collection Methods

This resulted in a total of 1,017 completed surveys.

#Read in data

## [1] 1017

For this report, we are particularly interested in those participants who (1) report that they own or currently have access to a firearm, and (2) report their current storage methods. After subsetting for participants that meet these criteria, there was a total of 363 surveys used in this analysis.

#Comparison is always between those that practice safe storage (case/safe/lockbox) and those that do not
    # In order to get this population of interest, we first subsetted data to only include those that respond "Yes" to "Do you currently own or have access to a firearm". This was question "BMW1"

data_owners=subset(data, BMW1 == "Yes") #this totals 363 participants
nrow(data_owners)

## [1] 363

Survey Considerations

The data also has a weighting value assigned to it. The variable “Weights” contains the proportional weighted value for each individual participant. This value will need to be applied to each measure variable evaluated below.

There is a “Skipped” option for each question asked, which is recorded if the individual did not respond. This will need to be accounted for in each question as well.

Results

#Question BMW3 is "How often do you carry a loaded firearm". 

BMW3.only= data_owners %>% 
  #Keep only the ID, Question/Responses for this question, Demographic information, and weighting score
      dplyr::select(1, 17, 15,69,47:65) %>%
    filter(BMW3 != "Skipped") %>% #filter out anyone who skipped the question 
  mutate(BMW5_2 = case_when(BMW5_2== "A safe, lockbox, or gun case (with a key, access code, or biometric lock)" ~ "1", TRUE ~ BMW5_2)) 
 

#Make columns with 0/1 binary data numerical variables in dataframe 
BMW3.only$BMW5_2=as.numeric(BMW3.only$BMW5_2)

#Recreate categorical variables within this analysis
BMW3.only$ppkid017 <- ifelse(BMW3.only$ppkid017 > 0, "yes", ifelse(BMW3.only$ppkid017 == 0, "no", BMW3.only$ppkid017))
BMW3.only$age.cat <- cut(BMW3.only$ppage,
                         breaks = c(18, 25, 40, 60, Inf),
                         labels = c("18-25", "26-40", "41-60", ">61"),
                         right = FALSE)


#get the number of participants that chose each response for graphing purposes
BMW3.only %>%
  dplyr::count(BMW3, BMW5_2) %>%
  tidyr::spread(key = BMW5_2, value = n, fill = 0)

##                                        BMW3  0  1 <NA>
## 1                                     Daily 14 23    0
## 2   At least once a week, but not every day  6 22    1
## 3 At least once a month, but not every week  6 13    0
## 4 At least once a year, but not every month  6  9    0
## 5                              Almost never 37 51    2
## 6 I own a firearm but never carry it loaded 75 90    1

# Group data by BMW3, calculate the weighted mean of BMW5_2 within each group using the survey weights, keeps one row per group, and outputs a summary table with the group label (response options) renamed as parameter.

chart=BMW3.only %>%
  group_by(BMW3) %>%
  dplyr::summarize(weighted_mean_BMW5_2 = weighted.mean(BMW5_2, Weights, na.rm = TRUE)) %>%
  ungroup() %>%
  filter(BMW3 != "Skipped")

chart

## # A tibble: 6 × 2
##   BMW3                                      weighted_mean_BMW5_2
##   <fct>                                                    <dbl>
## 1 Daily                                                    0.640
## 2 At least once a week, but not every day                  0.774
## 3 At least once a month, but not every week                0.681
## 4 At least once a year, but not every month                0.617
## 5 Almost never                                             0.586
## 6 I own a firearm but never carry it loaded                0.540

#Rename categoires for graphing
chart <- chart %>%
mutate(BMW3 = case_when(
    BMW3 == "I own a firearm but never carry it loaded" ~ "Never",
    BMW3 == "At least once a year, but not every month" ~ "Annually",
    BMW3 == "At least once a month, but not every week" ~ "Monthly",
    BMW3 == "At least once a week, but not every day" ~ "Weekly",
    TRUE ~ BMW3  # Keeps all other values as they are
  ))

#order levels for graphing
response_levels = c("Never", "Almost never", "Annually", "Monthly", "Weekly", "Daily")
chart$BMW3 <- factor(chart$BMW3, levels = response_levels)

#create plot showing frequency of carry among those with access to a firearm
ggplot(chart, aes(x = BMW3, y = weighted_mean_BMW5_2)) + 
  geom_line(group=1) + 
  geom_point() +  
  geom_label(aes(label = sprintf("%.3f", weighted_mean_BMW5_2)), nudge_y = 0.005, check_overlap = TRUE, 
            color = "black", size = 3, fontface = "bold", 
            label.padding = unit(0.2, "lines"),  
            label.size = 0.2) +  
  labs(x = "BMW3", y = "Weighted Mean BMW5_2", title = "Line Plot of BMW3 vs Weighted Mean BMW5_2") +
  theme_minimal()  +
  xlab("Frequency of Carry") +
  ggtitle("Proportion of Owners Using a Secure Container by Frequency of Carry") +
  ylab("Proportion Using Secure Container")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
          axis.title.x = element_text(size = 10, face = "bold"),
        axis.title.y = element_text(size = 10, face = "bold"),
        plot.title = element_text(size = 12, face = "bold", hjust = 0.5))

# Recode BMW3 levels for readability
BMW3.only <- BMW3.only %>%
  mutate(BMW3 = case_when(
    BMW3 == "I own a firearm but never carry it loaded" ~ "Never (47%)",
    BMW3 == "At least once a year, but not every month" ~ "Annually (4%)",
    BMW3 == "At least once a month, but not every week" ~ "Monthly (5%)",
    BMW3 == "At least once a week, but not every day" ~ "Weekly (7%)",
    BMW3 == "Almost every day" ~ "Daily (11%)",
    BMW3 == "Almost never" ~ "Almost never (26%)",
    BMW3 == "Daily" ~ "Daily (11%)",
    TRUE ~ BMW3
  ))

response_levels <- c(
  "Never (47%)", "Almost never (26%)", "Annually (4%)",
  "Monthly (5%)", "Weekly (7%)", "Daily (11%)"
)

BMW3.only$BMW3 <- factor(BMW3.only$BMW3, levels = response_levels)

# NEW: remove anything not in those levels (prevents NA categories)
BMW3.only <- BMW3.only %>% dplyr::filter(!is.na(BMW3))

chart <- BMW3.only %>%
  dplyr::group_by(BMW3) %>%
  dplyr::summarize(
    weighted_mean = weighted.mean(BMW5_2, Weights, na.rm = TRUE),
    n = dplyr::n(),              # CHANGED
    total_n = nrow(BMW3.only),
    .groups = "drop"
  ) %>%
  dplyr::mutate(prop_n = n / total_n)

# Plot with bubble sizes representing proportion of sample
ggsave(
  filename = "Figure 1.png",
  dpi = 300,
  width = 7,
  height = 5,
  units = "in",
  plot =
    ggplot(chart, aes(x = BMW3, y = weighted_mean, group = 1)) +
  geom_line(color = "#3F326C", size = 1) +
  geom_point(aes(size = prop_n), color = "#5C7750", alpha = 0.7) +
  geom_label(aes(label = scales::percent(weighted_mean, accuracy = 0.1)),
             nudge_y = 0.02,
             size = 3,
             fontface = "bold",
             label.padding = unit(0.2, "lines"),
             label.size = 0.2) +
  scale_size_continuous(range = c(3, 10), name = "Sample Proportion") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    x = "Frequency of Carry",
    y = "Proportion Using a Secure Container",
    title = "Proportion of Owners Using a Secure Container by Frequency of Carry"
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.title = element_text(size = 10, face = "bold"),
    plot.title = element_text(size = 12, face = "bold", hjust = 0.5),
    legend.position = "none"
  ))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

bubble_values <- chart %>%
  dplyr::select(BMW3, weighted_mean, n, total_n, prop_n) %>%
  dplyr::arrange(BMW3)

print(bubble_values)

## # A tibble: 6 × 5
##   BMW3               weighted_mean     n total_n prop_n
##   <fct>                      <dbl> <int>   <int>  <dbl>
## 1 Never (47%)                0.540   166     356 0.466 
## 2 Almost never (26%)         0.586    90     356 0.253 
## 3 Annually (4%)              0.617    15     356 0.0421
## 4 Monthly (5%)               0.681    19     356 0.0534
## 5 Weekly (7%)                0.774    29     356 0.0815
## 6 Daily (11%)                0.640    37     356 0.104

bubble_values %>%
  dplyr::mutate(
    weighted_mean_label = scales::percent(weighted_mean, accuracy = 0.1),
    prop_n_label = scales::percent(prop_n, accuracy = 0.1)
  ) %>%
  dplyr::select(BMW3, weighted_mean_label, n, prop_n_label) %>%
  print()

## # A tibble: 6 × 4
##   BMW3               weighted_mean_label     n prop_n_label
##   <fct>              <chr>               <int> <chr>       
## 1 Never (47%)        54.0%                 166 46.6%       
## 2 Almost never (26%) 58.6%                  90 25.3%       
## 3 Annually (4%)      61.7%                  15 4.2%        
## 4 Monthly (5%)       68.1%                  19 5.3%        
## 5 Weekly (7%)        77.4%                  29 8.1%        
## 6 Daily (11%)        64.0%                  37 10.4%

What other storage methods are used?

#All findings are subsetted by whether or not participants use secure storage, so first, BMW5 is process to get those populations
#Question BMW5 is "Do you use any of the following measures to safely store your firearms" and BMW5_2 specifies responses that do report the use of "A safe, lockbox, or gun case". Anyone who did not report the use of these tools has a value of "0". Those that do use these measures for safety have "A safe, lockbox, or gun case" in the cell.

BMW5.only= data_owners %>% 
  #Keep only the ID, Question/Responses for this question, Demographic information, and weighting score
    select(c(1, 69, 16:25, 47:65))  %>% 
  #filter out anyone who skipped the question 
  filter(BMW5_1 != "Skipped") %>% #this drops the sample size to 355 from 363
  #replace all text responses with a "1" value to allow adjusting by weight
  mutate(BMW5_1 = case_when(BMW5_1== "Keeping all firearms unloaded" ~ "1", TRUE ~ BMW5_1))  %>%
  mutate(BMW5_2 = case_when(BMW5_2== "A safe, lockbox, or gun case (with a key, access code, or biometric lock)" ~ "1", TRUE ~ BMW5_2))  %>%
  mutate(BMW5_3 = case_when(BMW5_3== "A cable or trigger lock" ~ "1", TRUE ~ BMW5_3))  %>%
  mutate(BMW5_4 = case_when(BMW5_4== "Motion alarm or notifications for firearm access" ~ "1", TRUE ~ BMW5_4))  %>%
  mutate(BMW5_5 = case_when(BMW5_5== "Firearm sensor for close range or self-aiming" ~ "1", TRUE ~ BMW5_5))  %>%
  mutate(BMW5_6 = case_when(BMW5_6== "Locking up ammunition separate from firearms" ~ "1", TRUE ~ BMW5_6))  %>%
  mutate(BMW5_7 = case_when(BMW5_7== "Disassembling firearms and storing parts separately" ~ "1", TRUE ~ BMW5_7))  %>%
  mutate(BMW5_8 = case_when(BMW5_8== "Entrusting keys or firearm parts to a trusted person or organization" ~ "1", TRUE ~ BMW5_8))  %>%
  mutate(BMW5_9 = case_when(BMW5_9== "Using other safety measure not listed here (please specify)" ~ "1", TRUE ~ BMW5_9))  %>%
  mutate(BMW5_10 = case_when(BMW5_10== "None of the above" ~ "1", TRUE ~ BMW5_10))

#Make columns with 0/1 binary data numerical variables in dataframe 
for (i in 3:12) {
  BMW5.only[[i]] <- as.numeric(BMW5.only[[i]])
}
#BMW5_2 is the lockbox/safe/case. Create a column that is yes/no to the use of safe/lockbox/case

BMW5.only$`Lock.Safe` <- ifelse(BMW5.only$BMW5_2 > 0, "yes", "no")

#Create categorical yes/no responses for kids present in home at over and under 18 years of age to use in analysis. Respondents reported the number of children in the home. We chose to use yes/no designation here.
BMW5.only$ppkid017.cat <- ifelse(BMW5.only$ppkid017 > 0, "yes", ifelse(BMW5.only$ppkid017 == 0, "no", BMW5.only$ppkid017))
BMW5.only$ppt18ov.cat <- ifelse(BMW5.only$ppt18ov > 0, "yes", ifelse(BMW5.only$ppt18ov == 0, "no", BMW5.only$ppt18ov))

#Create categorical variables for age within the dataset. 
BMW5.only$age.cat <- cut(BMW5.only$ppage,
                         breaks = c(18, 25, 40, 60, Inf),
                         labels = c("1. 18-25", "2. 26-40", "3. 41-60", "4. >61"),
                         right = FALSE)

library(reshape2)

#Question BMW5 is "Do you use any of the following measures to safely store your firearms" 
#BMW5_1 contains data for "Keeping all firearms unloaded"
#MBW5_2 A safe, lockbox, or gun case
#BMW5_3 contains data for "A cable or trigger lock"
#BMW4_4 contains data for "Motion alarm or notifications for firearm access"
#BMW5_5 contains data for "Firearm sensor for close range or self-aiming" 
#BMW5_6 contains data for "Locking up ammunition separate from firearms" 
#BMW5_7 contains data for "Disassembling firearms and storing parts separately" 
#BMW5_8 contains data for "Entrusting keys or firearm parts to a trusted person or organization" 
#BMW5_9 contains data for "Using other safety measures not listed here (please specify)" 
#BMW5_10 contains data for "None of the above" 

#BMW5_9_TEXT contains self reported methods of storage by those who selected "other" 

item_labels <- tibble::tribble(
  ~item,      ~measure,
  "BMW5_1",   "Keeping all firearms unloaded",
  "BMW5_3",   "Using a cable or trigger lock",
  "BMW5_4",   "Using motion alarm or access notifications",
  "BMW5_5",   "Firearm sensor for close-range or self-aiming",
  "BMW5_6",   "Locking up ammunition separate from firearms",
  "BMW5_7",   "Disassembling firearms and storing parts separately",
  "BMW5_8",   "Entrusting keys or parts to a trusted person or organization",
  "BMW5_9",   "Using other safety measures not listed here",
  "BMW5_10",  "None of the above"
)

#all users
BMW5_all=BMW5.only
sum(BMW5_all$Weights) #333.9897

## [1] 333.9897

#get sum of weights in each column. first multiply by 0, or 1. This will replace a weight with 0 for an individual who said no. And sum the individuals who said yes. 
sumBMW5all=BMW5_all %>%
  dplyr::summarise(across(starts_with("BMW5_"),
                   ~ sum(.x * Weights)))

#convert to proportions

percentBMW5all <- 100 * sumBMW5all / sum(BMW5_all$Weights)
percentBMW5all

##     BMW5_1   BMW5_2   BMW5_3   BMW5_4   BMW5_5   BMW5_6   BMW5_7   BMW5_8
## 1 51.19077 59.06157 22.06666 2.760355 1.148505 38.13443 6.562758 5.701756
##     BMW5_9  BMW5_10
## 1 3.138727 14.27945

#Only those who do not use a lockbox/gun safe
BMW5_secure= subset(BMW5.only, Lock.Safe == "no")
sum(BMW5_secure$Weights) #136.7302

## [1] 136.7302

sumBMW5sec=BMW5_secure %>%
  dplyr::summarise(across(starts_with("BMW5_"),
                   ~ sum(.x * Weights)))

percentBMW5sec <- 100 * sumBMW5sec / sum(BMW5_secure$Weights)
percentBMW5sec

##     BMW5_1 BMW5_2  BMW5_3    BMW5_4    BMW5_5   BMW5_6   BMW5_7   BMW5_8
## 1 44.55068      0 17.2489 0.5774933 0.8490075 27.12474 5.824486 1.966883
##     BMW5_9 BMW5_10
## 1 4.581305 34.8803

What other storage methods are being considered?

#Question BMW6 is "Which of the following methods for storing firearms would you consider using" 

#BMW6_1 contains data for "Keeping all firearms unloaded"
#BMW6_3 contains data for "A cable or trigger lock"
#BMW6_4 contains data for "Motion alarm or notifications for firearm access"
#BMW6_5 contains data for "Firearm sensor for close range or self-aiming" 
#BMW6_6 contains data for "Locking up ammunition separate from firearms" 
#BMW6_7 contains data for "Disassembling firearms and storing parts separately" 
#BMW6_8 contains data for "Entrusting keys or firearm parts to a trusted person or organization" 
#BMW6_9 contains data for "Using other safety measures not listed here (please specify)" 
#BMW6_10 contains data for "None of the above" 

#BMW6_9_TEXT contains self reported methods of storage by those who selected "other" 


BMW6.only= data_owners %>% 
  dplyr::select(c(1,17, 16:25, 27:36, 47:65, 69)) %>%
  filter(BMW6_1 != "Skipped" | is.na(BMW6_1)) #filter skipped by keep NA values

BMW6.only= BMW6.only %>% 
    mutate(BMW5_2 = case_when(BMW5_2== "A safe, lockbox, or gun case (with a key, access code, or biometric lock)" ~ "1", TRUE ~ BMW5_2)) %>%
  mutate(BMW6_1 = case_when(BMW6_1== "Keeping all firearms unloaded" ~ "1", TRUE ~ BMW6_1))  %>%
  mutate(BMW6_2 = case_when(BMW6_2== "A safe, lockbox, or gun case (with a key, access code, or biometric lock)" ~ "1", TRUE ~ BMW6_2))  %>%
  mutate(BMW6_3 = case_when(BMW6_3== "A cable or trigger lock" ~ "1", TRUE ~ BMW6_3))  %>%
  mutate(BMW6_4 = case_when(BMW6_4== "Motion alarm or notifications for firearm access" ~ "1", TRUE ~ BMW6_4))  %>%
  mutate(BMW6_5 = case_when(BMW6_5== "Firearm sensor for close range or self-aiming" ~ "1", TRUE ~ BMW6_5))  %>%
  mutate(BMW6_6 = case_when(BMW6_6== "Locking up ammo separate from firearms" ~ "1", TRUE ~ BMW6_6))  %>%
  mutate(BMW6_7 = case_when(BMW6_7== "Disassembling firearms and storing parts separately" ~ "1", TRUE ~ BMW6_7))  %>%
  mutate(BMW6_8 = case_when(BMW6_8== "Entrusting keys or firearm parts to a trusted person or organization" ~ "1", TRUE ~ BMW6_8))  %>%
  mutate(BMW6_9 = case_when(BMW6_9== "Using other safety measure not listed here (please specify)" ~ "1", TRUE ~ BMW6_9))  %>%
  mutate(BMW6_10 = case_when(BMW6_10== "None of the above" ~ "1", TRUE ~ BMW6_10))



#Make columns with 0/1 binary data numerical variables in dataframe 
for (i in 2:21) {
  BMW6.only[[i]] <- as.numeric(BMW6.only[[i]])
}

## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion

BMW6.only$ppkid017 <- ifelse(BMW6.only$ppkid017 > 0, "yes", ifelse(BMW6.only$ppkid017 == 0, "no", BMW6.only$ppkid017))

BMW6.only$age.cat <- cut(BMW6.only$ppage,
                         breaks = c(18, 25, 40, 60, Inf),
                         labels = c("18-25", "26-40", "41-60", ">61"),
                         right = FALSE)


BMW6.only$`Lock.Safe` <- ifelse(BMW6.only$BMW5_2 > 0, "yes", "no")

#because of survey logic, if they said "yes" to an option in BMW5 they were not shown it in BMW6. So we are replacing NA with "0". As they are not potential future adopters. 

BMW6.only <- BMW6.only %>%
  mutate(across(starts_with("BMW6_"),
                ~ replace_na(.x, 0)))


#all users
BMW6_all=BMW6.only
sum(BMW6_all$Weights) #333.9897

## [1] 339.6666

#get sum of weights in each column. first multiply by 0, or 1. This will replace a weight with 0 for an individual who said no. And sum the individuals who said yes. 
sumBMW6all=BMW6_all %>%
  dplyr::summarise(across(starts_with("BMW6_"),
                   ~ sum(.x * Weights)))

#convert to proportions

percentBMW6all <- 100 * sumBMW6all / sum(BMW6_all$Weights)
percentBMW6all

##     BMW6_1   BMW6_2   BMW6_3   BMW6_4   BMW6_5   BMW6_6   BMW6_7   BMW6_8
## 1 12.05202 21.09652 19.83761 18.82692 8.728971 14.21035 7.125905 4.844154
##     BMW6_9  BMW6_10
## 1 3.734571 35.38603

#Only those who do not use a lockbox/gun safe
BMW6_secure= subset(BMW6.only, Lock.Safe == "no")
sum(BMW6_secure$Weights) #136.7302

## [1] 136.7302

sumBMW6sec=BMW6_secure %>%
  dplyr::summarise(across(starts_with("BMW6_"),
                   ~ sum(.x * Weights)))

percentBMW6sec <- 100 * sumBMW6sec / sum(BMW6_secure$Weights)
percentBMW6sec

##     BMW6_1   BMW6_2   BMW6_3   BMW6_4   BMW6_5 BMW6_6   BMW6_7   BMW6_8
## 1 15.03943 50.35738 12.27476 11.45228 5.482342  15.78 3.552276 2.588313
##     BMW6_9  BMW6_10
## 1 4.117316 34.87712

What would motivate a change in gun safety behavior?

New Code:

We will remove anyone who said “strongly agree” or “somewhat agree” to option 7, because we are assuming that it should have been read as mutually exclusive. So we are looking for motivators of change in people who stated that they did not agree they would always keep at least one unlocked.

#Question BMW4 is "how much do you agree or disagree with each of the following". A series of potential motivators for gun safety are then presented. 

#BMW4_1 contains data for "If I had kids at home I would lock up my firearms"
#BMW4_2 contains data for "If someone in my home was struggling with mental or physical health, I would lock up my firearm(s)."
#BMW4_3 contains data for "If a close friend or family asked me to, I would lock up my firearm(s).
#BMW4_4 contains data for "To prevent my firearm(s) from being used, stolen, or damaged, I would lock it up."
#BMW4_5 contains data for "To prevent accidental injury from my firearm(s), I would lock it up."
#BMW4_6 contains data for "To prevent suicide for me or others, I would lock up my firearm(s)."
#BMW4_7 contains data for "I plan to always keep at least one firearm unlocked."


BMW4.only= data_owners %>% 
  #Keep only the ID, Question/Responses for this question, Demographic information, and weighting score
      dplyr::select(1, 17, 38:44,69,47:65) %>%
#filter out anyone who skipped the question 
  mutate(BMW5_2 = case_when(BMW5_2== "A safe, lockbox, or gun case (with a key, access code, or biometric lock)" ~ "1", TRUE ~ BMW5_2)) 
 

#Make columns with 0/1 binary data numerical variables in dataframe 
BMW4.only$BMW5_2=as.numeric(BMW4.only$BMW5_2)

## Warning: NAs introduced by coercion

BMW4.only$ppkid017 <- ifelse(BMW4.only$ppkid017 > 0, "yes", ifelse(BMW4.only$ppkid017 == 0, "no", BMW4.only$ppkid017))

BMW4.only$age.cat <- cut(BMW4.only$ppage,
                         breaks = c(18, 25, 40, 60, Inf),
                         labels = c("18-25", "26-40", "41-60", ">61"),
                         right = FALSE)



#Subset data to be those who disagree and those who agree with option 7 for planning to always keeping one unlocked 

BMW4_7T <- BMW4.only %>%
  filter(BMW4_7 %in% c("Strongly agree", "Somewhat agree"))

likert_long <- BMW4_7T %>%
  select(starts_with("BMW4_"), -BMW4_7, Weights) %>%
  pivot_longer(
    cols = starts_with("BMW4_"),
    names_to = "Question",
    values_to = "Response"
  )

plot_df <- likert_long %>%
  mutate(
    strong = Response == "Strongly agree",
    somewhat = Response == "Somewhat agree"
  ) %>%
  group_by(Question) %>%
  summarise(
    strong_pct = 100 * sum(Weights * strong, na.rm = TRUE) / sum(Weights, na.rm = TRUE),
    top2_pct = 100 * sum(Weights * (strong | somewhat), na.rm = TRUE) / sum(Weights, na.rm = TRUE)
  ) %>%
  ungroup()

plot_df <- plot_df %>%
  mutate(
    p = top2_pct / 100,
    n_eff = sum(BMW4_7T$Weights)^2 / sum(BMW4_7T$Weights^2),  # Kish effective N
    se = sqrt(p * (1 - p) / n_eff),
    low = 100 * (p - 1.96 * se),
    high = 100 * (p + 1.96 * se)
  )

plot_df <- plot_df %>%
  mutate(Question = case_when(
    Question == "BMW4_1" ~ "If I had kids at home...",
    Question == "BMW4_2" ~ "If someone in my home was struggling with mental or physical health...",
    Question == "BMW4_3" ~ "If a close friend or family asked me to...",
    Question == "BMW4_4" ~ "To prevent my firearm(s) from being used, stolen, or damaged...",
    Question == "BMW4_5" ~ "To prevent accidental injury from my firearm(s)...",
    Question == "BMW4_6" ~ "To prevent suicide for me or others...",
    TRUE ~ Question
  ))


motiv <- ggplot(plot_df, aes(y = reorder(Question, top2_pct))) +
  
  # Top 2 (dark red)
  geom_col(aes(x = top2_pct), fill = "#7E1417", width = 0.7) +
  
  # Strongly agree (dark blue overlay)
  geom_col(aes(x = strong_pct), fill = "#001042", width = 0.7) +
  
  geom_errorbarh(aes(xmin = low, xmax = high), height = 0.2) +
  
  geom_text(
    aes(x = -3, label = sprintf("%.1f%%", top2_pct)),
    hjust = 1,
    size = 4
  ) +
  
  scale_x_continuous(
    limits = c(-8, max(plot_df$high) + 8)
  ) +
  
  labs(
    title = "Agreement with Motivators Among Those Who Will Keep One Unlocked",
    subtitle = "Proportion of respondents who 
<span style='color:#001042; font-weight:bold;'>Somewhat Agree</span> 
or 
<span style='color:#7E1417; font-weight:bold;'>Strongly Agree</span>
<br>(Survey-weighted estimates with 95% confidence intervals)",
    y = NULL,
    x = "Percent"
  ) +
  
  theme_minimal(base_size = 14) +
  theme(
    plot.background  = element_rect(fill = "white", color = NA),
    panel.background = element_rect(fill = "white", color = NA),
    plot.title = element_text(size = 20, face = "bold", hjust = 0.5),
    plot.subtitle = ggtext::element_markdown(size = 15, hjust = 0.5),
    axis.text.y = element_text(size = 14),
    axis.text.x = element_text(size = 14),
    axis.title.x = element_text(size = 15),
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank()
  )

## Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
## ℹ Please use the `orientation` argument of `geom_errorbar()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

BMW4_7F <- BMW4.only %>%
  filter(BMW4_7 %in% c("Strongly disagree", "Somewhat disagree"))



likert_long2 <- BMW4_7F %>%
  select(starts_with("BMW4_"), -BMW4_7, Weights) %>%
  pivot_longer(
    cols = starts_with("BMW4_"),
    names_to = "Question",
    values_to = "Response"
  )

plot_df2 <- likert_long2 %>%
  mutate(
    strong = Response == "Strongly agree",
    somewhat = Response == "Somewhat agree"
  ) %>%
  group_by(Question) %>%
  summarise(
    strong_pct = 100 * sum(Weights * strong, na.rm = TRUE) / sum(Weights, na.rm = TRUE),
    top2_pct = 100 * sum(Weights * (strong | somewhat), na.rm = TRUE) / sum(Weights, na.rm = TRUE)
  ) %>%
  ungroup()

plot_df2 <- plot_df2 %>%
  mutate(
    p = top2_pct / 100,
    n_eff = sum(BMW4_7F$Weights)^2 / sum(BMW4_7F$Weights^2),  # Kish effective N
    se = sqrt(p * (1 - p) / n_eff),
    low = 100 * (p - 1.96 * se),
    high = 100 * (p + 1.96 * se)
  )

plot_df2 <- plot_df2 %>%
  mutate(Question = case_when(
    Question == "BMW4_1" ~ "If I had kids at home...",
    Question == "BMW4_2" ~ "If someone in my home was struggling with mental or physical health...",
    Question == "BMW4_3" ~ "If a close friend or family asked me to...",
    Question == "BMW4_4" ~ "To prevent my firearm(s) from being used, stolen, or damaged...",
    Question == "BMW4_5" ~ "To prevent accidental injury from my firearm(s)...",
    Question == "BMW4_6" ~ "To prevent suicide for me or others...",
    TRUE ~ Question
  ))


motiv2 <- ggplot(plot_df2, aes(y = reorder(Question, top2_pct))) +
  
  # Top 2 (dark red)
  geom_col(aes(x = top2_pct), fill = "#7E1417", width = 0.7) +
  
  # Strongly agree (dark blue overlay)
  geom_col(aes(x = strong_pct), fill = "#001042", width = 0.7) +
  
  geom_errorbarh(aes(xmin = low, xmax = high), height = 0.2) +
  
  geom_text(
    aes(x = -3, label = sprintf("%.1f%%", top2_pct)),
    hjust = 1,
    size = 4
  ) +
  
  scale_x_continuous(
    limits = c(-8, max(plot_df$high) + 8)
  ) +
  
  labs(
    title = "Agreement with Motivators Among Those Who Will Not Keep One Unlocked",
    subtitle = "Proportion of respondents who 
<span style='color:#001042; font-weight:bold;'>Somewhat Agree</span> 
or 
<span style='color:#7E1417; font-weight:bold;'>Strongly Agree</span>
<br>(Survey-weighted estimates with 95% confidence intervals)",
    y = NULL,
    x = "Percent"
  ) +
  
  theme_minimal(base_size = 14) +
  theme(
    plot.background  = element_rect(fill = "white", color = NA),
    panel.background = element_rect(fill = "white", color = NA),
    plot.title = element_text(size = 20, face = "bold", hjust = 0.5),
    plot.subtitle = ggtext::element_markdown(size = 15, hjust = 0.5),
    axis.text.y = element_text(size = 14),
    axis.text.x = element_text(size = 14),
    axis.title.x = element_text(size = 15),
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank()
  )

motiv2 <- motiv2 +
  labs(y = NULL) +
  theme(
    axis.text.y = element_blank(),
     plot.title = element_blank(),
    plot.subtitle = element_blank(),
  )

library(patchwork)

two_panel <- motiv + motiv2 +
  plot_layout(ncol = 2)

two_panel

## `height` was translated to `width`.
## `height` was translated to `width`.

ggsave(
  "motiv_two_panel.png",
  two_panel,
  width = 16,
  height = 8,
  dpi = 300
)

## `height` was translated to `width`.
## `height` was translated to `width`.

ggsave(
  "motiv.png",
  motiv,
  width = 12,
  height = 7,
  dpi = 600
)

## `height` was translated to `width`.

ggsave(
  "motiv_two.png",
  motiv2,
  width = 12,
  height = 7,
  dpi = 600
)

## `height` was translated to `width`.

RQ1a. Secure Container Utilization

Weighted percentages of respondents reporting the use of a safe, lockbox, or gun case were calculated using survey weights to account for sampling design. Weighted means and their 95% confidence intervals were estimated using the survey package in R, which incorporates respondent-level weights. Confidence intervals were derived using standard errors obtained from the weighted survey design, providing estimates that reflect the population distribution of the sampled respondents.

Survey Summaries

#This function computes survey-weighted outcome percentages with standard errors and confidence intervals by each demographic factor, then fits a survey-weighted logistic regression of the outcome (odds of secure container use) on that predictor and returns both results.

survey_summary <- function(var_name, design, outcome = "BMW5_2") {
  
  # Weighted percentages with SE and CI
  ci_table <- svyby(
    formula = as.formula(paste0("~", outcome)),
    by = as.formula(paste0("~", var_name)),
    design = design,
    FUN = svymean,
    vartype = c("se", "ci"),
    na.rm = TRUE
  ) %>%
    as.data.frame() %>%
    rename_with(~ c(var_name, "Weighted_Percentage", "SE", "CI_Lower", "CI_Upper"), everything()) %>%
    mutate(
      Weighted_Percentage = Weighted_Percentage * 100,
      CI_Lower = CI_Lower * 100,
      CI_Upper = CI_Upper * 100
    )
  
  # Survey-weighted logistic regression
  formula_str <- paste0(outcome, " ~ ", var_name)
  fit <- svyglm(as.formula(formula_str), design = design, family = quasibinomial())
  
  list(
    CI_Table = ci_table,
    Regression = list(
      Summary = summary(fit),
      CI = confint(fit)
    )
  )
}

# Create survey design object once
design <- svydesign(ids = ~1, data = BMW5.only, weights = ~Weights)

# List of variables for which we will run a regression
vars <- c("ppgender", "age.cat", "xparty4", "xurbanicity", "ppeducat", 
          "ppeduc5", "ppethm", "pphouse4", "ppinc7", "ppmarit5", 
          "ppmsacat", "ppreg4", "ppemploy", "ppkid017.cat", "ppt18ov")

# Apply function to all variables
results_list <- lapply(vars, function(v) survey_summary(v, design))

# Access example:
# Weighted CI table for age
results_list$age.cat$CI_Table   # won't work directly since lapply returns unnamed list
# Better to name the list elements:
names(results_list) <- vars

for (var in names(results_list)) {
  cat("\n==============================\n")
  cat("Variable:", var, "\n")
  cat("==============================\n")
  
  cat("\nWeighted Percentages + CI:\n")
  print(results_list[[var]]$CI_Table)
  
  cat("\nSurvey-weighted Logistic Regression:\n")
  print(results_list[[var]]$Regression$Summary)
  
  cat("\nConfidence Intervals for coefficients:\n")
  print(results_list[[var]]$Regression$CI)
  
  cat("\n\n")  # extra space between variables
}

Demographic Regresssion Models

ppgender

Weighted % (SE, CI)

ppgender	Weighted_Percentage	SE	CI_Lower	CI_Upper
Male	59.07602	0.03518634	52.17963	65.97242
Female	59.04026	0.04319632	50.57394	67.50658

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.44	0.37	0.15	0.08	0.65	0.012
ppgenderFemale	1.00	0.00	0.23	-0.45	0.45	0.995

age.cat

Weighted % (SE, CI)

age.cat	Weighted_Percentage	SE	CI_Lower	CI_Upper
18-25	73.65499	0.08509242	56.97718	90.33280
26-40	66.13982	0.05446808	55.46428	76.81537
41-60	60.11025	0.04339029	51.60591	68.61459
>61	50.53462	0.04963460	40.80641	60.26282

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	2.80	1.03	0.44	0.17	1.89	0.020
age.cat2. 26-40	0.70	-0.36	0.50	-1.34	0.63	0.475
age.cat3. 41-60	0.54	-0.62	0.47	-1.55	0.31	0.193
age.cat4. >61	0.37	-1.01	0.48	-1.95	-0.06	0.037

xparty4

Weighted % (SE, CI)

xparty4	Weighted_Percentage	SE	CI_Lower	CI_Upper
Republican	61.58013	0.04287909	53.17599	69.98428
Democrat	56.71812	0.06896706	43.20082	70.23541
Independent	57.15847	0.04612307	48.11852	66.19843
Something else	59.65833	0.09011666	41.99579	77.32087

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.60	0.47	0.18	0.12	0.83	0.010
xparty4Democrat	0.82	-0.20	0.33	-0.86	0.46	0.547
xparty4Independent	0.83	-0.18	0.26	-0.70	0.33	0.483
xparty4Something else	0.92	-0.08	0.42	-0.90	0.74	0.847

xurbanicity

Weighted % (SE, CI)

xurbanicity	Weighted_Percentage	SE	CI_Lower	CI_Upper
Urban	57.96007	0.05443195	47.29161	68.62854
Rural	58.72802	0.05132055	48.66938	68.78667
Suburban	59.88209	0.03995091	52.05185	67.71232

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.38	0.32	0.22	-0.12	0.76	0.151
xurbanicityRural	1.03	0.03	0.31	-0.57	0.64	0.918
xurbanicitySuburban	1.08	0.08	0.28	-0.47	0.63	0.776

ppeducat

Weighted % (SE, CI)

ppeducat	Weighted_Percentage	SE	CI_Lower	CI_Upper
Less than HS	33.17008	0.11649244	10.33798	56.00218
HS	56.91554	0.05112485	46.89526	66.93583
Some college	56.18232	0.04972332	46.43673	65.92791
Bachelor or higher	68.49396	0.04210514	60.24151	76.74642

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	0.50	-0.70	0.53	-1.73	0.33	0.183
ppeducatHS	2.66	0.98	0.57	-0.13	2.09	0.084
ppeducatSome college	2.58	0.95	0.56	-0.16	2.06	0.093
ppeducatBachelor or higher	4.38	1.48	0.56	0.37	2.58	0.009

ppeduc5

Weighted % (SE, CI)

ppeduc5	Weighted_Percentage	SE	CI_Lower	CI_Upper
No high school diploma or GED	33.17008	0.11649244	10.33798	56.00218
High school graduate (high school diploma or the equivalent GED)	56.91554	0.05112485	46.89526	66.93583
Some college or Associate degree	56.18232	0.04972332	46.43673	65.92791
Bachelor’s degree	64.11467	0.05686230	52.96987	75.25948
Master’s degree or above	74.44852	0.06225743	62.24629	86.65075

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	0.50	-0.70	0.53	-1.73	0.33	0.183
ppeduc5High school graduate (high school diploma or the equivalent GED)	2.66	0.98	0.57	-0.13	2.09	0.084
ppeduc5Some college or Associate degree	2.58	0.95	0.56	-0.16	2.06	0.093
ppeduc5Bachelor’s degree	3.60	1.28	0.58	0.14	2.42	0.028
ppeduc5Master’s degree or above	5.87	1.77	0.62	0.55	2.99	0.005

ppethm

Weighted % (SE, CI)

ppethm	Weighted_Percentage	SE	CI_Lower	CI_Upper
White, Non-Hispanic	58.37785	0.03004375	52.48939	64.26632
Black or African American, Non-Hispanic	53.98844	0.10690988	33.03449	74.94239
Other, Non-Hispanic	78.57190	0.13820299	51.48461	105.65919
Hispanic	63.04836	0.09353498	44.71584	81.38088
2+ races, Non-Hispanic	46.41445	0.13193034	20.55658	72.27232

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.40	0.34	0.12	0.10	0.58	0.007
ppethmBlack or African American, Non-Hispanic	0.84	-0.18	0.45	-1.06	0.70	0.691
ppethmOther, Non-Hispanic	2.61	0.96	0.83	-0.67	2.59	0.248
ppethmHispanic	1.22	0.20	0.42	-0.63	1.02	0.641
ppethm2+ races, Non-Hispanic	0.62	-0.48	0.54	-1.55	0.59	0.377

pphouse4

Weighted % (SE, CI)

pphouse4	Weighted_Percentage	SE	CI_Lower	CI_Upper
One-family house detached from any other house	61.06382	0.02940491	55.30056	66.82707
One-family condo or townhouse attached to other units	56.08947	0.12023152	32.52452	79.65441
Building with 2 or more apartments	38.98875	0.10243937	18.91100	59.06649
Other (mobile home, boat, RV, van, etc.)	49.01316	0.14027631	21.51951	76.50681

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.57	0.45	0.12	0.21	0.69	0.000
pphouse4One-family condo or townhouse attached to other units	0.81	-0.21	0.50	-1.20	0.79	0.684
pphouse4Building with 2 or more apartments	0.41	-0.90	0.45	-1.78	-0.02	0.046
pphouse4Other (mobile home, boat, RV, van, etc.)	0.61	-0.49	0.57	-1.62	0.64	0.395

ppinc7

Weighted % (SE, CI)

ppinc7	Weighted_Percentage	SE	CI_Lower	CI_Upper
Under $10,000	56.82471	0.19122233	19.34582	94.30360
$10,000 to $24,999	40.15789	0.10239880	20.08810	60.22769
$25,000 to $49,999	44.98968	0.07934691	29.43797	60.54139
$50,000 to $74,999	52.18102	0.06580500	39.28348	65.07856
$75,000 to $99,999	54.79359	0.07507078	40.07999	69.50719
$100,000 to $149,999	67.81940	0.05472496	57.09350	78.54529
$150,000 or more	70.16699	0.04924848	60.51446	79.81951

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.32	0.27	0.78	-1.26	1.81	0.725
ppinc7$10,000 to $24,999	0.51	-0.67	0.89	-2.42	1.07	0.449
ppinc7$25,000 to $49,999	0.62	-0.48	0.84	-2.13	1.18	0.573
ppinc7$50,000 to $74,999	0.83	-0.19	0.82	-1.81	1.43	0.820
ppinc7$75,000 to $99,999	0.92	-0.08	0.84	-1.73	1.56	0.922
ppinc7$100,000 to $149,999	1.60	0.47	0.82	-1.14	2.08	0.566
ppinc7$150,000 or more	1.79	0.58	0.81	-1.02	2.18	0.476

ppmarit5

Weighted % (SE, CI)

ppmarit5	Weighted_Percentage	SE	CI_Lower	CI_Upper
Now married	63.06410	0.03413009	56.37472	69.75347
Widowed	35.22435	0.12035218	11.63576	58.81294
Divorced	56.40292	0.08177034	40.37623	72.42962
Separated	57.33557	0.22982424	12.29084	102.38029
Never married	54.91611	0.06018871	43.11934	66.71288

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.71	0.53	0.15	0.25	0.82	0.000
ppmarit5Widowed	0.32	-1.14	0.55	-2.22	-0.07	0.037
ppmarit5Divorced	0.76	-0.28	0.36	-0.99	0.44	0.446
ppmarit5Separated	0.79	-0.24	0.95	-2.11	1.63	0.801
ppmarit5Never married	0.71	-0.34	0.28	-0.90	0.22	0.235

ppmsacat

Weighted % (SE, CI)

ppmsacat	Weighted_Percentage	SE	CI_Lower	CI_Upper
Non-metro area	57.7626	0.06326162	45.36355	70.16165
Metro area	59.3854	0.03021909	53.46257	65.30824

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.37	0.31	0.26	-0.2	0.82	0.228
ppmsacatMetro area	1.07	0.07	0.29	-0.5	0.63	0.816

ppreg4

Weighted % (SE, CI)

ppreg4	Weighted_Percentage	SE	CI_Lower	CI_Upper
Northeast	65.25129	0.07496860	50.55772	79.94487
Midwest	65.81647	0.05368496	55.29441	76.33853
South	54.03173	0.04101341	45.99325	62.07021
West	59.11143	0.06474228	46.42218	71.80069

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.88	0.63	0.33	-0.02	1.28	0.058
ppreg4Midwest	1.03	0.03	0.41	-0.78	0.83	0.951
ppreg4South	0.63	-0.47	0.37	-1.20	0.26	0.206
ppreg4West	0.77	-0.26	0.43	-1.10	0.58	0.539

ppemploy

Weighted % (SE, CI)

ppemploy	Weighted_Percentage	SE	CI_Lower	CI_Upper
Working full-time	63.26282	0.03662169	56.08510	70.44054
Working part-time	66.40362	0.07817543	51.08152	81.72572
Not working	51.09609	0.04640956	41.99998	60.19220

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.72	0.54	0.16	0.23	0.85	0.001
ppemployWorking part-time	1.15	0.14	0.38	-0.62	0.89	0.720
ppemployNot working	0.61	-0.50	0.24	-0.98	-0.02	0.041

ppkid017.cat

Weighted % (SE, CI)

ppkid017.cat	Weighted_Percentage	SE	CI_Lower	CI_Upper
no	54.62022	0.03302592	48.14726	61.09318
yes	69.74103	0.04693341	60.54225	78.93981

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.20	0.19	0.13	-0.08	0.45	0.165
ppkid017.catyes	1.91	0.65	0.26	0.14	1.16	0.013

ppt18ov

Weighted % (SE, CI)

ppt18ov	Weighted_Percentage	SE	CI_Lower	CI_Upper
1	49.96925	0.06990993	36.26716	63.67135
2	61.10978	0.03448161	54.35151	67.86805
3	63.16869	0.07174878	49.10619	77.23120
4	52.85846	0.09927923	33.40009	72.31683
5	46.03779	0.31157239	-15.02928	107.10485
6	100.00000	0.00000000	100.00000	100.00000

Logistic regression (OR)

term	OR	Estimate	SE	CI_Lower	CI_Upper	p_value
(Intercept)	1.23	0.21	0.31	-0.41	0.83	0.504
ppt18ov	1.07	0.07	0.13	-0.19	0.34	0.600

Associated Materials: Supplemental Figure 1

design <- svydesign(ids = ~1, data = BMW5.only, weights = ~Weights)

# Survey-weighted logistic regression
fit_age <- svyglm(BMW5_2 ~ age.cat, design = design, family = quasibinomial())

# Estimated marginal means (predicted probabilities) by age group
emm_age <- emmeans(fit_age, ~ age.cat, type = "response")  # response scale = probabilities

emm_df <- as.data.frame(emm_age) %>%
  transmute(
    age.cat,
    p = prob,
    lower = asymp.LCL,
    upper = asymp.UCL,
    pct = 100 * p,
    pct_low = 100 * lower,
    pct_up  = 100 * upper
  )

pairs_age <- pairs(emm_age, adjust = "holm")  # or "tukey" if you prefer
pairs_df <- as.data.frame(pairs_age)

pairs_df

##  contrast                odds.ratio        SE  df null z.ratio p.value
##  (1. 18-25) / (2. 26-40)   1.431298 0.7177270 Inf    1   0.715  0.7841
##  (1. 18-25) / (3. 41-60)   1.855310 0.8801432 Inf    1   1.303  0.5918
##  (1. 18-25) / 4. >61       2.736631 1.3173605 Inf    1   2.091  0.2190
##  (2. 26-40) / (3. 41-60)   1.296243 0.3929562 Inf    1   0.856  0.7841
##  (2. 26-40) / 4. >61       1.911993 0.6003166 Inf    1   2.064  0.2190
##  (3. 41-60) / 4. >61       1.475026 0.3962663 Inf    1   1.447  0.5918
## 
## P value adjustment: holm method for 6 tests 
## Tests are performed on the log odds ratio scale

p1 <- ggplot(emm_df, aes(x = age.cat, y = pct, group = 1)) +
  geom_line() +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = pct_low, ymax = pct_up), width = 0.12) +
   geom_text(
    aes(label = sprintf("%.1f%%", pct)),
    vjust = -1.1,
    hjust= -0.3,
    size = 3.5
  ) +
  labs(
    x = "Age category",
    y = "Predicted probability of using a secure container (%)",
    title = "Secure Container storage by age category",
    subtitle = "Predicted probabilities from survey-weighted logistic regression (95% CI)"
  ) +
  theme_classic()

p1

ggsave(
  filename = "Supp Fig 1.png",
  plot=p1,
  dpi = 600,
  width = 7,
  height = 5,
  units = "in")

Secure Container Storage

Abby Beatty

2026-02-28

Background

Collection Methods

Survey Considerations

Results

What other storage methods are used?

What other storage methods are being considered?

What would motivate a change in gun safety behavior?

New Code:

RQ1a. Secure Container Utilization

Survey Summaries

Demographic Regresssion Models

ppgender

age.cat

xparty4

xurbanicity

ppeducat

ppeduc5

ppethm

pphouse4

ppinc7

ppmarit5

ppmsacat

ppreg4

ppemploy

ppkid017.cat

ppt18ov

Associated Materials: Supplemental Figure 1

RQ1b. Frequency of Carry