Analysis of the Motivations, Preferences, and Behaviors of Anglers During the COVID-19 Pandemic

This is an analysis of a Texas anglers during the COVID-19 pandemic. The data set includes a group of anglers that recruited during the 4 years prior to the pandemic and a group of anglers that recruited during the pandemic. Additionally, a portion of these anglers identified COVID-19 as a primary driver in their motivation to start fishing. The purpose of this analysis is to characterize angler motivations, preferences, and behaviors during the pandemic, with special emphasis on anglers that recruited due to the pandemic.

library(dplyr) 
library(pander)
library(knitr)
library(tidyverse)
library(readxl)
library(naniar)
library(vegan)
library(ggplot2)
library(effectsize)
library(survey)
library(car)
library(caret) 
library(regclass) 
library(generalhoslem)
library(effects)
library(kableExtra)
library(ade4)
library(ggrepel)
library(corrplot)
library(sjPlot)
library(ltm)
library(psych)
library(wordcloud)
library(RColorBrewer)
library(tm)

Import data and recode factors

### Import the data and revise variables ###
raw1 <- as.data.frame(unclass(read_xlsx("REOrg TPWD NEW ANGLERS Includes INCOMPLETES.xlsx", sheet = "Data")),stringsAsFactors = TRUE) # COVID anglers data (includes incompletes)
rand_draw1 <- as.data.frame(unclass(read_xlsx("RAND062021_angler_rand_draw.xlsx")),stringsAsFactors = TRUE) # random draw sample for COVID and preCOVID anglers data
raw2 <- as.data.frame(unclass(read_xlsx("REOrg TPWD PreCOVID ANGLERS includes INCOMPLETES.xlsx", sheet = "Data")),stringsAsFactors = TRUE) # preCOVID anglers data (includes incompletes)
rawA <- as.data.frame(unclass(read_xlsx("MATCHING_ANGLERS_EXPORT.xlsx", sheet = "UNIQUE_IDS")),stringsAsFactors = TRUE) # License database information to join by emails

# This key indicating whether anglers were coastal or inland anglers was created from 'COVID_anglers_join_POSREP_V2.R'
coastal_inland_key <- as.data.frame(unclass(read_csv("Coastal_inland_anglers_key.csv")),stringsAsFactors = TRUE) 

# This list indicating whether anglers purchased a license in 2021 and 2022 was created from 'COVID_anglers_join_POSREP_V2.R'
LCYR_2021_2022 <- as.data.frame(unclass(read_xlsx("License_buyers_2021_2022.xlsx")),stringsAsFactors = TRUE) 

# Concatenate COVID and PreCOVID datasets
raw1 <- raw1 %>% 
  mutate(New.Old = rep("New",length(RespID))) # New anglers that started after COVID
    
raw2 <- raw2 %>% 
  mutate(New.Old = rep("Old",length(RespID))) # Old anglers that started preCOVID

raw3 <- rbind(raw1,raw2)

# Make emails uppercase to join with POSREP database (separate file brought into SAS)
raw3 <- raw3 %>%
  mutate(Q0_Email = as.character(toupper(Q0_Email)))

write.csv(raw3,"Incompletes_new_and_retained.csv") # save an object to csv

# Add in coastal vs inland identifiers for anglers
raw4 <- raw3 %>%
  left_join(coastal_inland_key[2:3], by = "RespID") %>%
  distinct()

# Add in whether customers purchased a license in 2021 or 2022
LCYR_2021_2022 <- LCYR_2021_2022 %>% replace_with_na_all(condition = ~.x %in% c("NA")) # replaces missing values
LCYR_2021_2022 <- LCYR_2021_2022 %>%
  mutate(across(LCYR_2021:LCYR_2022,~ str_replace_all(., "NA", "None"))) # select and replace strings with parentheses

raw5 <- raw4 %>%
  left_join(LCYR_2021_2022, by = "RespID") 

# Clean independent variables
rev1 <- raw5 %>%
  replace_with_na_at(.vars = "Q28_Gender",condition = ~.x %in% c("Other/nonbinary")) %>%
  replace_with_na_all(condition = ~.x %in% c("Prefer not to answer")) %>%
  mutate(Children = as.factor(case_when(Q29B %in% c("Children under age 5","Children age 5 to 12","Children age 13 to 17") ~ "Under 17", 
                              Q29E == "Children 18 years or older" ~ "18 or over",
                              Q29A == "No. I do not have any children." ~ "No children")),
         Race = as.factor(case_when(Q26C == "Black or African American" ~ "Black or African American",
                               Q26D == "Native American or Alaskan Native" ~ "Native American or Alaskan Native",
                               Q26E == "Asian or Pacific Islander" ~ "Asian or Pacific Islander",
                               Q26B == "White" ~ "White",
                               str_length(Q26Other)>0 ~ "Other")),
         Hispanic = as.factor(case_when(Q25_Hisp %in% c("Yes, Mexican, Mexican American, Chicano","Yes, other Spanish/Hispanic group or mixed heritage (Please specify.)") ~ "Hispanic",
                              Q25_Hisp == "No, not Spanish/Hispanic" ~ "Not Hispanic")),
         R_E = as.factor(paste(Race,Hispanic)),
         Race_Ethn = as.factor(case_when(R_E %in% c("Other Hispanic","Asian or Pacific Islander Hispanic","Black or African American Hispanic","Native American or Alaskan Native Hispanic") ~ "Other Hispanic",
                                R_E %in% c("Asian or Pacific Islander NA","Asian or Pacific Islander Not Hispanic") ~ "Asian or Pacific Islander",
                                R_E %in% c("Black or African American NA","Black or African American Not Hispanic") ~ "Black or African American",
                                R_E %in% c("Native American or Alaskan Native NA","Native American or Alaskan Native Not Hispanic") ~ "Native American or Alaskan Native",
                                R_E %in% c("NA Hispanic") ~ "Hispanic",
                                R_E %in% c("White NA","White Not Hispanic") ~ "White",
                                R_E %in% c("White Hispanic") ~ "White Hispanic")),
         SW_FW1 = paste(Q7A,Q7B,Q7C,Q7D,Q7E,sep = "_"),
         SW_FW2 = as.factor(case_when(grepl("Saltwater",SW_FW1) & grepl("Freshwater", SW_FW1) ~ "SW.FW",
                                                     grepl("Saltwater", SW_FW1) ~ "SW",
                                                     grepl("Freshwater", SW_FW1) ~ "FW")),
         Income = as.factor(Q30_Income),
         Gender = as.factor(Q28_Gender),
         Age = Q27_Age,
         New.Old = as.factor(New.Old),
         COVID = as.factor(Q4_COVID_Fished),
         Coastal = as.factor(COASTAL),
         UniqueID = paste(rep("A",length(RespID)),RespID,sep = "_")) %>%
  droplevels()

Calculate non-response bias by gender and age

# Survey response proportions
responses <- rev1 %>%
  filter(!(is.na(Gender) | is.na(Age))) %>% # remove NA values from Gender & Age
  count(Gender,Age) %>%
  mutate(freq = prop.table(n)) # get proportions for Gender and Age

# Random draw proportions
rand_draw2 <- rand_draw1 %>%
    mutate(Age = Age_Grp) %>%
  filter(!(is.na(Gender) | is.na(Age))) %>% # remove NA values from Gender & Age
  count(Gender,Age) %>%
  mutate(freq = prop.table(n),
         Gender = str_replace_all(Gender,"F", "Female"), # get proportions for Gender and Age 
         Gender = str_replace_all(Gender,"M", "Male"))

# Weights for Age and Gender survey data
weights1 <- responses %>%
  inner_join(rand_draw2, by = c("Gender","Age"),suffix = c("_survey", "_random")) %>%
  mutate(weights = freq_random/freq_survey)

weights_F <- weights1 %>%
  filter(Gender == "Female")

weights_M <- weights1 %>%
  filter(Gender == "Male")

weights2 <- weights_F %>%
  inner_join(weights_M,by = "Age",suffix = c("Female","Male")) %>%
  mutate(Survey_Female = freq_surveyFemale*100,
         Survey_Male = freq_surveyMale*100,
         Database_Female = freq_randomFemale*100,
         Database_Male = freq_randomMale*100,
         Weights_Female = weightsFemale,
         Weights_Male = weightsMale) %>%
  dplyr::select(Age,Survey_Female,Database_Female,Weights_Female,Survey_Male,Database_Male,Weights_Male)

weights2 %>% 
  kbl(caption = "Gender and Age Frequencies and Weights",digits = 3) %>%
  kable_classic(full_width = F, html_font = "Cambria")

Gender and Age Frequencies and Weights

Age	Survey_Female	Database_Female	Weights_Female	Survey_Male	Database_Male	Weights_Male
18-24	2.072	5.291	2.553	2.710	13.663	5.042
25-34	7.333	8.071	1.101	6.589	14.207	2.156
35-44	9.086	8.517	0.937	12.593	14.849	1.179
45-54	8.555	6.456	0.755	16.206	11.566	0.714
55-64	8.130	3.971	0.488	13.496	7.326	0.543
65 or older	4.091	1.983	0.485	9.139	4.100	0.449

# Weights for Coastal and Inland counties
weights_coast1 <- rev1 %>%
  filter(!(is.na(Coastal))) %>% # remove NA values from Coastal variable
  count(Coastal) %>%
  mutate(freq = prop.table(n)) # get proportions for Gender and Age
# 56% Inland counties and 44% Coastal counties
# There is a response bias because the original split was 50% Inland, 50% Coast

# Coastal weights table
weights_coast2 <- weights_coast1 %>%
  mutate(weights = 0.50/freq) %>%
  dplyr::select(-freq,-n)

## Adjust responses using Ageweights
rev1_wts <- rev1 %>%
  left_join(weights1, by = c("Gender","Age")) %>%
  left_join(weights_coast2, by = "Coastal") %>%
  replace_na(list(weights.x = 1, weights.y = 1)) %>%
  mutate(Gender = as.factor(Gender),
         weights = weights.x * weights.y) %>%
  dplyr::select(-weights.x,-weights.y,-n_survey,-freq_survey,-n_random,-freq_random)

The survey was sent to a random sample from the fishing license database that was 34% female and 66% male, while 39% of survey respondents were female and 61% were male. Ages 18-44 were underrepresented in the survey respondents, while ages 45-65+ were over-represented and weighted accordingly.The most heavily weighted groups were male respondents aged 18-24 yrs followed by female respondents aged 18-24 yrs. The groups that were most over-represented in the survey respondents were males and females aged 65 or older.

The survey was stratified to include a 50:50 split of inland and coastal counties. The survey respondents were split 56% from inland counties and 44% from coastal counties. Thus, coastal county responses were more heavily weighted to take into account the non-response bias.

Weighted demographics of survey respondents

rev2_wts <- rev1_wts
rev2_wts$Income <- factor(rev2_wts$Income,levels = c("Less than $20,000", "Between $20,000 and $39,999", "Between $40,000 and $59,999","Between $60,000 and $79,999","Between $80,000 and $99,999","Between $100,000 and $149,999","Over $150,000"))

Income <- tab_xtab(rev2_wts$Income,rev2_wts$New.Old,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)

Children <- tab_xtab(rev2_wts$Children,rev2_wts$New.Old,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)

Gender <- tab_xtab(rev2_wts$Gender,rev2_wts$New.Old,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)

Age <- tab_xtab(rev2_wts$Age,rev2_wts$New.Old,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)

Race_Ethn <- tab_xtab(rev2_wts$Race_Ethn,rev2_wts$New.Old,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)

Children

Children	New.Old		Total
	New	Old
18 or over	347 22.6 %	208 13.5 %	555 36.1 %
No children	547 35.6 %	179 11.6 %	726 47.2 %
Under 17	172 11.2 %	85 5.5 %	257 16.7 %
Total	1066 69.3 %	472 30.7 %	1538 100 %

Income

Income	New.Old		Total
	New	Old
Less than $20,000	69 4.6 %	22 1.5 %	91 6.1 %
Between $20,000 and $39,999	121 8.2 %	52 3.5 %	173 11.7 %
Between $40,000 and $59,999	141 9.5 %	49 3.3 %	190 12.8 %
Between $60,000 and $79,999	186 12.5 %	67 4.5 %	253 17 %
Between $80,000 and $99,999	146 9.8 %	59 4 %	205 13.8 %
Between $100,000 and $149,999	220 14.8 %	98 6.6 %	318 21.4 %
Over $150,000	168 11.3 %	86 5.8 %	254 17.1 %
Total	1051 70.8 %	433 29.2 %	1484 100 %

Race_Ethn

Race_Ethn	New.Old		Total
	New	Old
Asian or Pacific Islander	66 3.7 %	10 0.6 %	76 4.3 %
Black or African American	62 3.5 %	16 0.9 %	78 4.4 %
Hispanic	31 1.8 %	17 1 %	48 2.8 %
Native American or Alaskan Native	35 2 %	11 0.6 %	46 2.6 %
Other Hispanic	70 4 %	45 2.5 %	115 6.5 %
White	753 42.6 %	331 18.7 %	1084 61.3 %
White Hispanic	226 12.8 %	93 5.3 %	319 18.1 %
Total	1243 70.4 %	523 29.6 %	1766 100 %

Age 

Age	New.Old		Total
	New	Old
18-24	274 14.4 %	81 4.3 %	355 18.7 %
25-34	312 16.4 %	110 5.8 %	422 22.2 %
35-44	310 16.3 %	132 7 %	442 23.3 %
45-54	217 11.4 %	125 6.6 %	342 18 %
55-64	141 7.4 %	83 4.4 %	224 11.8 %
65 or older	68 3.6 %	45 2.4 %	113 6 %
Total	1322 69.7 %	576 30.3 %	1898 100 %

Gender

Gender	New.Old		Total
	New	Old
Female	509 26.8 %	137 7.2 %	646 34 %
Male	815 43 %	435 22.9 %	1250 65.9 %
Total	1324 69.8 %	572 30.2 %	1896 100 %

We surveyed a total of 2,238 anglers with 68% recruited during the COVID-19 pandemic (hereafter new anglers), while 32% were anglers with fishing licenses during the 4 years prior to the start of the pandemic (hereafter retained anglers). Approximately half of the anglers surveyed (47%) had no children, 17% had children under 17 yrs old, and 36% had adult children. The survey respondents were primarily composed of anglers that identified as White (61%), followed by White Hispanic (18%), Black or African American (4%), Other Hispanic (7%), Asian or Pacific Islander (4%), Hispanic (3%), and Native American or Alaskan Native (3%). Finally, the majority of survey respondents had incomes over $80K (52%), while 42% of respondents had incomes between $20K-$79K, and only 6% had incomes less than $20K. Most retained anglers were male (66%), 35 years or older (66.8%), White (69%), with incomes over $80,000 (56.1%), and had adult children (44%).

Q0. Are there differences between the new and retained anglers?

# Set reference levels to typical angler (i.e. White, Male, etc.)
rev2_wts$New.Old <- relevel(rev2_wts$New.Old, ref = "Old")
rev2_wts$COVID <- relevel(rev2_wts$COVID, ref = "No, I was already planning to go fishing.")
rev2_wts$Income <- relevel(rev2_wts$Income, ref = "Between $100,000 and $149,999")
rev2_wts$Race_Ethn <- relevel(rev2_wts$Race_Ethn, ref = "White")
rev2_wts$Age <- relevel(rev2_wts$Age, ref = "45-54")
rev2_wts$SW_FW2 <- relevel(rev2_wts$SW_FW2, ref = "FW")
rev2_wts$Children2 <- relevel(rev2_wts$Children, ref = "18 or over")
rev2_wts$Gender <- relevel(rev2_wts$Gender, ref = "Male")

# Survey logistic regression 
svy_desn_Q0 <- svydesign(ids = ~0,  weights = ~weights,data = rev2_wts) # Survey design
Q0_mod1 <- svyglm(New.Old ~ Gender + Income + Children + Race_Ethn + SW_FW2 + Age + COVID  + Coastal, design=svy_desn_Q0, family=quasibinomial)
summary(Q0_mod1) # examine model results
varImp(Q0_mod1) # Variable importance
VIF(Q0_mod1) # collinearity
Q0_mod2 <- svyglm(New.Old ~ Gender + Race_Ethn + SW_FW2 + Age, design=svy_desn_Q0, family=quasibinomial)
summary(Q0_mod2)
VIF(Q0_mod2) # collinearity
car::Anova(Q0_mod2, type = 3, test.statistic = "F",error.df = degf(Q0_mod2$survey.design)) # Evaluate variable significance

# Estimate and odds ratio table
Q0_coef_table <- coef(summary(Q0_mod2))
Q0_AOR <- standardize_parameters(Q0_mod2, exp = TRUE) # Odds_ratio and 95% CI
Q0_summary_table <- cbind(Q0_coef_table,Odds_Ratio = Q0_AOR$Std_Odds_Ratio,CI_low = Q0_AOR$CI_low,CI_high = Q0_AOR$CI_high)

# Check model accuracy
logitgof(obs = Q0_mod2$y, fitted(Q0_mod2)) #run hoslem test to test goodness of fit
glm.probs <- predict(Q0_mod2,type = "response") #Predict() assigns fitted model values for each participants based on the model
glm.pred <- as.factor(ifelse(glm.probs > 0.5, 1, 0)) #If prediction is above 50% assign 1 and assign 0 if prediction is below 50%
prop.table(table(glm.pred == Q0_mod2$y)) 

# Survey response proportions
Q0_prop <- rev2_wts %>%
  filter(!(is.na(New.Old))) %>% # remove NA values from Gender & Age
  count(New.Old) %>%
  mutate(freq = prop.table(n)) # get proportions

# Table of responses
Q0_prop %>% 
  kbl(caption = "New vs Retained Anglers",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") # 81% were already planning to go fishing

New vs Retained Anglers

New.Old	n	freq
Old	721	0.32
New	1517	0.68

# Model odds ratio results
Q0_summary_table %>% 
  kbl(caption = "Q0 Model Odds Ratios",digits = 3) %>%
  kable_classic(full_width = F, html_font = "Cambria")

Q0 Model Odds Ratios

	Estimate	Std. Error	t value	Pr(>|t|)	Odds_Ratio	CI_low	CI_high
(Intercept)	0.501	0.151	3.328	0.001	1.651	1.229	2.218
GenderFemale	0.644	0.140	4.592	0.000	1.903	1.446	2.506
Race_EthnAsian or Pacific Islander	1.033	0.365	2.833	0.005	2.809	1.374	5.742
Race_EthnBlack or African American	0.623	0.320	1.948	0.052	1.864	0.996	3.488
Race_EthnHispanic	-0.113	0.596	-0.189	0.850	0.893	0.278	2.876
Race_EthnNative American or Alaskan Native	0.245	0.610	0.401	0.688	1.278	0.386	4.229
Race_EthnOther Hispanic	-0.343	0.358	-0.957	0.339	0.710	0.352	1.433
Race_EthnWhite Hispanic	0.182	0.208	0.875	0.382	1.200	0.797	1.805
SW_FW2SW	-0.680	0.184	-3.695	0.000	0.507	0.353	0.727
SW_FW2SW.FW	-0.346	0.166	-2.082	0.038	0.707	0.511	0.980
Age18-24	1.000	0.325	3.081	0.002	2.719	1.438	5.142
Age25-34	0.518	0.200	2.592	0.010	1.678	1.134	2.483
Age35-44	0.343	0.166	2.066	0.039	1.409	1.017	1.950
Age55-64	0.147	0.165	0.888	0.375	1.158	0.838	1.600
Age65 or older	0.078	0.200	0.390	0.696	1.081	0.730	1.600

# Barplots of responses
model.frame(Q0_mod2) %>% ggplot(aes(x = New.Old, fill = Age)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "New vs Retained Anglers")

model.frame(Q0_mod2) %>% ggplot(aes(x = New.Old, fill = Gender)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "New vs Retained Anglers")

model.frame(Q0_mod2) %>% ggplot(aes(x = New.Old, fill = Race_Ethn)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "New vs Retained Anglers")

model.frame(Q0_mod2) %>% ggplot(aes(x = New.Old, fill = SW_FW2)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "New vs Retained Anglers")

# Effect plot
plot(Effect(focal.predictors = c("Age"),Q0_mod2))

plot(Effect(focal.predictors = c("Gender"),Q0_mod2))

plot(Effect(focal.predictors = c("SW_FW2"),Q0_mod2))

plot(Effect(focal.predictors = c("Race_Ethn"),Q0_mod2))

A total of 68% of survey respondents were new (recruited during the COVID-19 pandemic), while 32% of respondents were retained from the 4 years prior to the pandemic.

For all of the logistic regressions that follow in this analysis, the reference levels will be set for the typical angler represented in this data set as follows: White, male, aged 45-54, with adult children, freshwater angler, income between $100-$150K, retained angler, and not prompted by COVID-19 to fish.

Compared to the typical angler (i.e. White, Male, aged 45-54, etc.), new anglers that recruited during the pandemic were: 1.903 times more likely to be female, 1.974 times more likely to fish in freshwater compared to saltwater, 1.413 times more likely to fish in freshwater compared to fishing in both saltwater and freshwater, 2.809 and 1.864 times more likely to identify as Asian or Pacific Islander and Black or African American, respectively, and 1.409 to 2.719 times more likely to be between the ages of 18-44.

Q1. How likely are you to go fishing in the next 6 months (Very likely to very unlikely)?

rev2_wts$Q1_Fish_6months <- recode_factor(rev2_wts$Q1_Fish_6months,"Somewhat likely" = "Likely","Very likely" = "Likely","Very unlikely" = "Unlikely")

# Survey response proportions
Q1_prop <- rev2_wts %>%
  filter(!(is.na(Q1_Fish_6months))) %>% # remove NA values from Gender & Age
  count(Q1_Fish_6months) %>%
  mutate(freq = prop.table(n)) # 99% report they are likely to go fishing in the next 6 months

Responses are summarized into likely or unlikely. Most anglers (99%) are likely to go fishing in the next 6 months, further analysis is not needed.

Q2. Why are you unlikely to go fishing in the next 6 months?

Q2_responses <- rev1_wts %>%
  mutate(Q2_all = as.factor(case_when(Q2B == "It’s too hard to find a place to fish." ~ "Hard to find fishing place", 
                              Q2C == "No one in my family or friend circle wants to go with me." ~ "No one to fish with",
                              Q2D == "I don’t know what to do with fish when I catch them." ~ "Don't know what to do with fish",
                              Q2E == "I decided I just don’t like fishing." ~ "Don't like fishing",
                              str_length(Q2Other)>0 ~ "Other"))) %>%
  filter(!is.na(Q2_all))

ggplot(Q2_responses,aes(Q2_all)) + geom_histogram(stat = "count") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8),axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + labs(y = "Counts", x = "Q2. Reasons Unlikely to Go Fishing") # histogram counts of levels

Only 19 responses, not enough data for robust analysis

Q3. Did you go fishing between mid-March 2020 and today?

svy_desn_Q3 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts) # Survey design
Q3_mod1 <- svyglm(Q3_Fished ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age + Coastal, design=svy_desn_Q3, family=quasibinomial)
summary(Q3_mod1) # examine model results
varImp(Q3_mod1) # Variable importance
VIF(Q3_mod1) # collinearity
Q3_mod2 <- svyglm(Q3_Fished ~ Gender + Age, design=svy_desn_Q3, family=quasibinomial)
summary(Q3_mod2)
Q3_mod3 <- svyglm(Q3_Fished ~ Age, design=svy_desn_Q3, family=quasibinomial)
summary(Q3_mod3)
car::Anova(Q3_mod3, type = 3, test.statistic = "F",error.df = degf(Q3_mod3$survey.design)) # Evaluate variable significance

# Estimate and odds ratio table
Q3_coef_table <- coef(summary(Q3_mod3))
Q3_AOR <- standardize_parameters(Q3_mod3, exp = TRUE) # Odds_ratio and 95% CI
Q3_summary_table <- cbind(Q3_coef_table,Odds_Ratio = Q3_AOR$Std_Odds_Ratio,CI_low = Q3_AOR$CI_low,CI_high = Q3_AOR$CI_high)

# Check model accuracy
logitgof(obs = Q3_mod3$y, fitted(Q3_mod3)) #run hoslem test to test goodness of fit
glm.probs <- predict(Q3_mod3,type = "response") #Predict() assigns fitted model values for each participants based on the model
glm.pred <- as.factor(ifelse(glm.probs > 0.5, 1, 0)) #If prediction is above 50% assign 1 and assign 0 if prediction is below 50%
prop.table(table(glm.pred == Q3_mod3$y)) # Correct assignment 88% of the time 
Q3_age_xtab <- xtabs(~Q3_Fished + Age, data = rev2_wts)
Q3_AOR_18 <- 1/0.27

# Survey response proportions
Q3_prop <- rev2_wts %>%
  filter(!(is.na(Q3_Fished))) %>% # remove NA values from Gender & Age
  count(Q3_Fished) %>%
  mutate(freq = prop.table(n)) # get proportions

# Table of responses
Q3_prop %>% 
  kbl(caption = "Q3 Survey Responses",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") # 88% went fishing in the last 6 months

Q3 Survey Responses

Q3_Fished	n	freq
No	256	0.12
Yes	1965	0.88

# Model odds ratios
Q3_summary_table %>% 
  kbl(caption = "Q3 Model Odds Ratios",digits = 3) %>%
  kable_classic(full_width = F, html_font = "Cambria")

Q3 Model Odds Ratios

	Estimate	Std. Error	t value	Pr(>|t|)	Odds_Ratio	CI_low	CI_high
(Intercept)	1.997	0.146	13.678	0.000	7.369	5.534	9.813
Age18-24	0.336	0.436	0.770	0.441	1.399	0.595	3.289
Age25-34	0.970	0.337	2.879	0.004	2.637	1.362	5.104
Age35-44	0.486	0.237	2.050	0.041	1.626	1.021	2.588
Age55-64	-0.025	0.209	-0.119	0.906	0.976	0.647	1.470
Age65 or older	-0.927	0.207	-4.476	0.000	0.396	0.264	0.594

# Histogram of responses
model.frame(Q3_mod3) %>% ggplot(aes(x = Q3_Fished, fill = Age)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q3. Fished in Last Year")

# Model effect plot
plot(Effect(focal.predictors = c("Age"),Q3_mod3))

Most respondents (88%) went fishing between mid-March 2020 and 2021. Compared to the typical angler aged 45-54 yrs, younger individuals (25-34 yrs) were 2.637 times more likely to have fished, while older individuals (65+ yrs) were 2.527 times less likely to have fished in the last year.

Q4. Did the effects of the COVID-19 pandemic prompt you to go fishing?

svy_desn_Q4 <- svydesign(ids = ~0,  weights = ~weights,data = rev2_wts) # Survey design
Q4_mod1 <- svyglm(Q4_COVID_Fished ~ Gender + Income + Children + New.Old + Race_Ethn + SW_FW2 + Age + Coastal, design=svy_desn_Q4, family=quasibinomial)
summary(Q4_mod1) # examine model results
varImp(Q4_mod1) # Variable importance
VIF(Q4_mod1) # collinearity
Q4_mod2 <- svyglm(Q4_COVID_Fished ~ Age + Children + Gender + Income + New.Old, design=svy_desn_Q4, family=quasibinomial)
summary(Q4_mod2)
Q4_mod3 <- svyglm(Q4_COVID_Fished ~ Children + Gender + New.Old, design=svy_desn_Q4, family=quasibinomial)
summary(Q4_mod3)
varImp(Q4_mod3) # Variable importance
VIF(Q4_mod3) # collinearity
car::Anova(Q4_mod3, type = 3, test.statistic = "F",error.df = degf(Q4_mod3$survey.design)) # Evaluate variable significance

# Estimate and odds ratio table
Q4_coef_table <- coef(summary(Q4_mod3))
Q4_AOR <- standardize_parameters(Q4_mod3, exp = TRUE) # Odds_ratio and 95% CI
Q4_summary_table <- cbind(Q4_coef_table,Odds_Ratio = Q4_AOR$Std_Odds_Ratio,CI_low = Q4_AOR$CI_low,CI_high = Q4_AOR$CI_high)

# Check model accuracy
logitgof(obs = Q4_mod3$y, fitted(Q4_mod3)) #run hoslem test to test goodness of fit
glm.probs <- predict(Q4_mod3,type = "response") #Predict() assigns fitted model values for each participants based on the model
glm.pred <- as.factor(ifelse(glm.probs > 0.5, 1, 0)) #If prediction is above 50% assign 1 and assign 0 if prediction is below 50%
prop.table(table(glm.pred == Q4_mod3$y)) 
Q4_xtab <- xtabs(~Q4_COVID_Fished + Age, data = rev2_wts)
Q4_xtab <- Q4_xtab%>%ftable(row.vars=c("Age"))

# Survey response proportions
Q4_prop <- rev2_wts %>%
  filter(!(is.na(Q4_COVID_Fished))) %>% # remove NA values from Gender & Age
  count(Q4_COVID_Fished) %>%
  mutate(freq = prop.table(n)) # get proportions

# Table of responses
Q4_prop %>% 
  kbl(caption = "Q4 Survey Responses",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") # 81% were already planning to go fishing

Q4 Survey Responses

Q4_COVID_Fished	n	freq
No, I was already planning to go fishing.	1583	0.81
Yes	381	0.19

# Model odds ratio results
Q4_summary_table %>% 
  kbl(caption = "Q4 Model Odds Ratios",digits = 3) %>%
  kable_classic(full_width = F, html_font = "Cambria")

Q4 Model Odds Ratios

	Estimate	Std. Error	t value	Pr(>|t|)	Odds_Ratio	CI_low	CI_high
(Intercept)	-2.417	0.209	-11.565	0.000	0.089	0.059	0.134
ChildrenNo children	0.762	0.181	4.202	0.000	2.143	1.501	3.059
ChildrenUnder 17	0.509	0.235	2.162	0.031	1.663	1.048	2.639
GenderFemale	0.637	0.180	3.545	0.000	1.892	1.329	2.692
New.OldNew	0.401	0.212	1.891	0.059	1.494	0.985	2.265

# Barplot of responses
model.frame(Q4_mod3) %>% ggplot(aes(x = Q4_COVID_Fished, fill = Gender)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q4. Fished due to COVID-19")

model.frame(Q4_mod3) %>% ggplot(aes(x = Q4_COVID_Fished, fill = Children)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q4. Fished due to COVID-19")

model.frame(Q4_mod3) %>% ggplot(aes(x = Q4_COVID_Fished, fill = New.Old)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q4. Fished due to COVID-19")

# Effect plot
plot(Effect(focal.predictors = c("Children"),Q4_mod3))

plot(Effect(focal.predictors = c("Gender"),Q4_mod3))

plot(Effect(focal.predictors = c("New.Old"),Q4_mod3))

Most anglers (81%) were already planning on going fishing and were not prompted by COVID-19 to go fishing.

The odds of going fishing due to COVID were: 2.143 times more likely for anglers with no children compared to those with adult children, 1.663 times more likely for anglers with children under 17 years compared to those with adult children, 1.892 times more likely for female anglers compared to males, 1.494 times more likely for new anglers compared to retained anglers.

Q5. How important were each of the following reasons why you went fishing? Please select no more than 2 of these for the “Very Important” column.

# Clean Q5 survey answers
rev2_wts <- rev2_wts %>%  
  mutate_at(vars(Q5A:Q5I), list(~recode_factor(.,"Not Important" = "1","Somewhat Important" = "2","Important" = "3","Very Important (please limit to 2 reasons)" = "4"))) %>%
  mutate(across(Q5A:Q5I, ~as.numeric(as.character(.)))) %>%
  mutate_at(vars(Q5A:Q5I), list(~.*weights)) %>%
  mutate(Food = Q5A,
         Sport = Q5B,
         Outdoors = Q5C,
         Nature = Q5D,
         Friends = Q5E,
         Children1 = Q5F,
         Stress = Q5G,
         Relax = Q5H,
         Get_away = Q5I)

# Remove rows with missing independent variable data
Q5_df1 <- rev2_wts %>%
  filter_at(vars(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2,Coastal),all_vars(!is.na(.))) 
Q5_df2 <- Q5_df1 %>%
  dplyr::select(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2,Coastal,Food,Sport,Outdoors,Nature,Friends,Children1,Stress,Relax,Get_away)
# Select independent variables
Q5_indep1 <- Q5_df1 %>%
  dplyr::select(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2,Coastal)
# Select dependent variables
Q5_dep1 <- Q5_df1 %>%
  dplyr::select(Food,Sport,Outdoors,Nature,Friends,Children1,Stress,Relax,Get_away)
# Transform dependent data
Q5_dep2 <- decostand(Q5_dep1, method = 'normalize') 
# MCA, keep 2 axes, on independent data
Q5_indep2 <- dudi.mix(Q5_indep1, scannf=F, nf=2)$tab # MCA, keep 2 axes

# Run RDA
Q5_rda_mod1 <- rda(Q5_dep2 ~ ., data = Q5_indep2) 
envfit(Q5_rda_mod1,Q5_indep2) # select variables at alpha = 0.001 for next model
Q5_indep3 <- Q5_indep2 %>%
  dplyr::select(Gende.Female,Gende.Male,Age.18.24,Age.25.34,Age.65.or.older,Child.18.or.over, Child.No.children,Child.Under.17,COVID.Yes,New.O.Old)
Q5_rda_mod2 <- rda(Q5_dep2 ~ ., data = Q5_indep3)
envfit(Q5_rda_mod2,Q5_indep3) # select variables at alpha = 0.001 for next model
RsquareAdj(Q5_rda_mod2) 
sqrt(vif.cca(Q5_rda_mod2))

# Extract scores for visualization
vare_spp_sco1 <- scores(Q5_rda_mod2, display = "species")
vare_sam_sco1 <- scores(Q5_rda_mod2, display = "sites")
vare_env_sco1 <- scores(Q5_rda_mod2, display = "bp")
vare_spp_tbl1 <- as_tibble(vare_spp_sco1)
vare_sam_tbl1 <- as_tibble(vare_sam_sco1)
vare_env_tbl1 <- as_tibble(vare_env_sco1)
vare_spp_tbl1 <- mutate(vare_spp_tbl1, vgntxt=rownames(vare_spp_sco1), ccatype = "species")
vare_sam_tbl1 <- mutate(vare_sam_tbl1, vgntxt=rownames(vare_sam_sco1), ccatype = "sites")
vare_env_tbl1 <- mutate(vare_env_tbl1, vgntxt=rownames(vare_env_sco1), ccatype = "bp")
vare_env_tbl1 <- mutate(vare_env_tbl1, vgntxt=c("Female","Age 18-24","Age 25-34","Age 65 over","Children 18 over","No Children","COVID-Yes","Retained"),ccatype = "bp")
rescaled1 <- vare_env_tbl1 %>% 
            dplyr::select(RDA1, RDA2) %>%
            as.matrix() *0.75
vare_tbl <- dplyr::select(vare_env_tbl1, vgntxt, ccatype) %>%
            bind_cols(as_tibble(rescaled1)) %>%
            bind_rows(vare_spp_tbl1)

# Get mean, SD, and calculate cronbach's alpha
Q5_summary1 <- rev2_wts %>%
  dplyr::select(Food,Sport,Outdoors,Nature,Friends,Children1,Stress,Relax,Get_away)

Q5_summary2 <- describe(Q5_summary1,na.rm = TRUE,IQR = TRUE) 
Q5_summary2$Variables <- row.names(Q5_summary2) 
Q5_summary3 <- Q5_summary2[,c("Variables","median","IQR")]
Q5_summary3$Variables<-str_replace_all(Q5_summary3$Variables, c("_" = " ")) # replaces spaces in names with periods

# Calculate cronbach's alpha
Q5_cronbach.alpha <- cronbach.alpha(abs(Q5_summary1),na.rm = TRUE,standardized = TRUE)
Q5_cronbach.alpha$alpha

# Pull out significant variable terms for each axis
Q5_envfit1 <- envfit(Q5_rda_mod2,Q5_indep3) # select variables at alpha = 0.001 for next model
Q5_envfit2 <- data.frame(Q5_envfit1$vectors[1],Q5_envfit1$vectors[2],Q5_envfit1$vectors[4])
names(Q5_envfit2) <- c("RDA 1","RDA 2","r2","p")
Q5_envfit2$Variables <- row.names(Q5_envfit2) 
Q5_envfit2 <- Q5_envfit2 %>% relocate(Variables, .before = "RDA 1")

# Descriptive stats table
tab_df(Q5_summary3)

Variables	median	IQR
Food	1.24	1.48
Sport	1.86	1.87
Outdoors	2.67	2.02
Nature	2.58	2.02
Friends	2.54	1.82
Children1	2.04	2.15
Stress	2.51	1.90
Relax	2.58	1.85
Get away	1.91	2.13

# Reorganize data to plot
Q5_a <- rev1_wts[c(which(colnames(rev1_wts)=="Q5A"):which(colnames(rev1_wts)=="Q5I"))]
colnames(Q5_a) <- c("Food","Sport","Outdoors","Nature","Friends","Children1","Stress","Relax","Get_away")
Q5_b <- Q5_a %>% 
  pivot_longer(cols = Food:Get_away,names_to = "Q5_answers",values_to = "Agreement")
Q5_c <- Q5_b %>%
  group_by(Q5_answers,Agreement) %>%
  count() %>%
  na.omit()
Q5_c$Agreement <- ordered(Q5_c$Agreement,levels = c("Not Important", "Somewhat Important", "Important","Very Important (please limit to 2 reasons)"))

# All anglers
ggplot(Q5_c, aes(fill=Agreement, y=n, x=Q5_answers)) + geom_bar(position="fill", stat="identity") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8),axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + labs(y = "Percent", x = "Q5. Reasons for Fishing") 

# PERMANOVA
anova.cca(Q5_rda_mod2, step = 1000, by = "axis",permutations = 999) # 4 sig axes

## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = Q5_dep2 ~ Gende.Female + Gende.Male + Age.18.24 + Age.25.34 + Age.65.or.older + Child.18.or.over + Child.No.children + Child.Under.17 + COVID.Yes + New.O.Old, data = Q5_indep3)
##           Df Variance       F Pr(>F)    
## RDA1       1 0.005048 59.7942  0.001 ***
## RDA2       1 0.001506 17.8377  0.001 ***
## RDA3       1 0.000475  5.6203  0.008 ** 
## RDA4       1 0.000223  2.6446  0.351    
## RDA5       1 0.000082  0.9767  0.987    
## RDA6       1 0.000036  0.4306  1.000    
## RDA7       1 0.000009  0.1015  1.000    
## RDA8       1 0.000001  0.0076  1.000    
## Residual 915 0.077253                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_df(Q5_envfit2)

Variables	RDA.1	RDA.2	r2	p
Gende.Female	-0.44	0.90	0.09	0.00
Gende.Male	0.44	-0.90	0.09	0.00
Age.18.24	-0.96	0.27	0.02	0.00
Age.25.34	-0.98	0.21	0.05	0.00
Age.65.or.older	0.81	-0.58	0.01	0.00
Child.18.or.over	1.00	-0.10	0.14	0.00
Child.No.children	-1.00	0.09	0.26	0.00
Child.Under.17	1.00	-0.05	0.02	0.00
COVID.Yes	-0.89	0.46	0.03	0.00
New.O.Old	0.67	-0.74	0.03	0.00

# Basic plot for visualization
ggplot() +
  geom_point(aes(x=RDA1, y=RDA2), data=filter(vare_tbl, ccatype=="species"))  +
  geom_text_repel(aes(x=RDA1, y=RDA2, label=vgntxt, size=3,colour=ccatype),
                  data=vare_tbl, seed=123) + 
  geom_segment(aes(x=0, y=0, xend=RDA1, yend=RDA2), arrow=arrow(length = unit(0.2,"cm")),
               data=filter(vare_tbl, ccatype=="bp"), color="blue") +
  coord_fixed() +
  theme_classic() +
  scale_colour_manual(values = c("blue", "black")) +
  theme(legend.position="none")

Most anglers indicated that being outdoors in nature, relaxing, relieving stress, and spending time with friends were all important reasons for fishing, while fishing for food, sport, or spending time with children were less important reasons to go fishing.

Anglers over 65 and retained anglers were relatively more likely to indicate an increased preference for fishing for food or sport. While most anglers indicated that being outdoors, enjoying nature, and relieving stress were important factors, anglers between the ages of 18-34 were the most likely to list these reasons as important factors. Anglers that recruited during COVID and female anglers were more likely to fish to spend time with friends and family. Finally, anglers with adult children were more likely than anglers without children to report that they fish to spend time with children.

Cronbach’s alpha for Q5 is 0.969647

Q6 Since mid-March 2020, how many days have you gone fishing?

# Order the factor
rev2_wts$Q6_DaysFished <- as.factor(rev2_wts$Q6_DaysFished)
rev2_wts$Q6_DaysFished  <- ordered(rev2_wts$Q6_DaysFished, levels = c("Once (1 day)", "2–3 days", "4–10 days","More than 10 days"))
svy_desn_Q6 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts,na.rm = TRUE) # New survey design

# Ordinal regression
Q6_mod1 <- svyolr(Q6_DaysFished ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age + Coastal, design=svy_desn_Q6)
summary(Q6_mod1) # examine model results
Q6_mod2 <- svyolr(Q6_DaysFished ~ SW_FW2 + Gender + Age + Income + Race_Ethn, design=svy_desn_Q6)
summary(Q6_mod2) # examine model results
Q6_mod3 <- svyolr(Q6_DaysFished ~ SW_FW2 + Gender + Age + Race_Ethn, design=svy_desn_Q6)
summary(Q6_mod3) # examine model results
VIF(Q6_mod3) # collinearity

# Estimate and odds ratio table
Q6_coef_table <- coef(summary(Q6_mod3))
Q6_pval <- pnorm(abs(Q6_coef_table[, "t value"]),lower.tail = FALSE)* 2
Q6_AOR <- exp(coef(Q6_mod3))
Q6_summary_table <- cbind(Q6_coef_table, "p value" = round(Q6_pval,3),Odds_Ratio = Q6_AOR)

# Survey response proportions
Q6_prop <- rev2_wts %>%
  filter(!(is.na(Q6_DaysFished))) %>% # remove NA values from Gender & Age
  count(Q6_DaysFished) %>%
  mutate(freq = prop.table(n)) # get proportions

# Table of responses
Q6_prop %>% 
  kbl(caption = "Q6 Days Fished",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") # 81% were already planning to go fishing

Q6 Days Fished

Q6_DaysFished	n	freq
Once (1 day)	97	0.05
2–3 days	441	0.23
4–10 days	630	0.33
More than 10 days	726	0.38

# Model odds ratio
Q6_summary_table %>% 
  kbl(caption = "Q6 Model Odds Ratios",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") # 81% were already planning to go fishing

Q6 Model Odds Ratios

	Value	Std. Error	t value	p value	Odds_Ratio
SW_FW2SW	-0.77	0.16	-4.88	0.00	0.46
SW_FW2SW.FW	0.49	0.14	3.59	0.00	1.64
GenderFemale	-0.32	0.12	-2.72	0.01	0.73
Age18-24	0.30	0.25	1.21	0.23	1.36
Age25-34	0.21	0.17	1.28	0.20	1.24
Age35-44	-0.02	0.14	-0.17	0.87	0.98
Age55-64	-0.11	0.14	-0.78	0.43	0.89
Age65 or older	-0.43	0.18	-2.37	0.02	0.65
Race_EthnAsian or Pacific Islander	-0.13	0.26	-0.51	0.61	0.88
Race_EthnBlack or African American	-0.65	0.27	-2.46	0.01	0.52
Race_EthnHispanic	-0.86	0.59	-1.47	0.14	0.42
Race_EthnNative American or Alaskan Native	0.72	0.52	1.37	0.17	2.06
Race_EthnOther Hispanic	0.24	0.34	0.71	0.48	1.28
Race_EthnWhite Hispanic	-0.14	0.17	-0.83	0.41	0.87
Once (1 day)|2–3 days	-3.33	0.19	-17.46	0.00	0.04
2–3 days|4–10 days	-1.22	0.15	-8.29	0.00	0.29
4–10 days|More than 10 days	0.24	0.14	1.73	0.08	1.27

# Histogram of responses
model.frame(Q6_mod3) %>% ggplot(aes(x = Q6_DaysFished, fill = Gender)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q6. Days Fished")

model.frame(Q6_mod3) %>% ggplot(aes(x = Q6_DaysFished, fill = Age)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q6. Days Fished")

model.frame(Q6_mod3) %>% ggplot(aes(x = Q6_DaysFished, fill = Race_Ethn)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q6. Days Fished")

model.frame(Q6_mod3) %>% ggplot(aes(x = Q6_DaysFished, fill = SW_FW2)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q6. Days Fished")

Most anglers were likely to have gone fishing more than 1 day (94%).

The odds of fishing the most often (>10 days) were: 2.152 times less likely for anglers that fish in saltwater compared to freshwater anglers, 1.639 times more likely for anglers that fish in both freshwater and saltwater compared to just freshwater anglers, 1.376 times less likely for females compared to males, 1.535 times less likely for the oldest anglers (65 yrs+) compared to anglers aged 45-54 yrs, and 1.918 times less likely for Black or African American anglers compared to White anglers.

Q7 From mid-March 2020 until today, select all the places you have fished.

# Proportions by waterbody type
SW_FW2 <- rev2_wts %>%
  filter(!(is.na(SW_FW2))) %>% 
  count(SW_FW2) %>%
  mutate(freq = prop.table(n)) # get proportions

# Chi-square tests
svy_desn_Q7 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts) # Survey design

# SW/FW
Q7_New.Old_SW_FW2_chisq <- svychisq(~New.Old + SW_FW2, svy_desn_Q7, statistic = "adjWald") # YES significant
Q7_Gender_SW_FW2_chisq <- svychisq(~Gender + SW_FW2, svy_desn_Q7, statistic = "adjWald") # YES significant
Q7_Age_SW_FW2_chisq <- svychisq(~Age + SW_FW2, svy_desn_Q7, statistic = "adjWald") # YES significant
Q7_Children_SW_FW2_chisq <- svychisq(~Children + SW_FW2, svy_desn_Q7, statistic = "adjWald") # YES significant
Q7_Race_Ethn_SW_FW2_chisq <- svychisq(~Race_Ethn + SW_FW2, svy_desn_Q7, statistic = "adjWald") # YES significant
Q7_Income_SW_FW2_chisq <- svychisq(~Income + SW_FW2, svy_desn_Q7, statistic = "adjWald") # NOT significant
Q7_Coastal_SW_FW2_chisq <- svychisq(~Coastal + SW_FW2, svy_desn_Q7, statistic = "adjWald") # YES significant

# Contingency tables for significant factors
Q7_New.Old_xtab <- tab_xtab(rev2_wts$New.Old,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q7_Gender_xtab <- tab_xtab(rev2_wts$Gender,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q7_Age_xtab <- tab_xtab(rev2_wts$Age,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q7_Children_xtab <- tab_xtab(rev2_wts$Children,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q7_Race_Ethn_xtab <- tab_xtab(rev2_wts$Race_Ethn,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q7_Income_xtab <- tab_xtab(rev2_wts$Income,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q7_Coastal_xtab <- tab_xtab(rev2_wts$Coastal,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE)

# correlation plots of residuals
Q7corr_New.Old_SW_FW2 <- corrplot(Q7_New.Old_SW_FW2_chisq$residuals, is.cor = FALSE)

Q7corr_Gender_SW_FW2 <- corrplot(Q7_Gender_SW_FW2_chisq$residuals, is.cor = FALSE)

Q7corr_Age_SW_FW2 <- corrplot(Q7_Age_SW_FW2_chisq$residuals, is.cor = FALSE)

Q7corr_Children_SW_FW2 <- corrplot(Q7_Children_SW_FW2_chisq$residuals, is.cor = FALSE)

Q7corr_Race_Ethn_SW_FW2 <- corrplot(Q7_Race_Ethn_SW_FW2_chisq$residuals, is.cor = FALSE)

Q7corr_Coastal_SW_FW2 <- corrplot(Q7_Coastal_SW_FW2_chisq$residuals, is.cor = FALSE)

# Combine chi-squared test results into dataframe
Q7_Gender_SW_FW <- data.frame(cbind("SW/FW vs Gender",Q7_Gender_SW_FW2_chisq[[c(2,1)]],Q7_Gender_SW_FW2_chisq$statistic,Q7_Gender_SW_FW2_chisq$p.value))
Q7_Age_SW_FW <- data.frame(cbind("SW/FW vs Age",Q7_Age_SW_FW2_chisq[[c(2,1)]],Q7_Age_SW_FW2_chisq$statistic,Q7_Age_SW_FW2_chisq$p.value))
Q7_Children_SW_FW <- data.frame(cbind("SW/FW vs Children",Q7_Children_SW_FW2_chisq[[c(2,1)]],Q7_Children_SW_FW2_chisq$statistic,Q7_Children_SW_FW2_chisq$p.value))
Q7_Race_Ethn_SW_FW <- data.frame(cbind("SW/FW vs Race/Ethnicity",Q7_Race_Ethn_SW_FW2_chisq[[c(2,1)]],Q7_Race_Ethn_SW_FW2_chisq$statistic,Q7_Race_Ethn_SW_FW2_chisq$p.value))
Q7_New.Old_SW_FW <- data.frame(cbind("SW/FW vs New/Retained",Q7_New.Old_SW_FW2_chisq[[c(2,1)]],Q7_New.Old_SW_FW2_chisq$statistic,Q7_New.Old_SW_FW2_chisq$p.value))
Q7_Coastal_SW_FW <- data.frame(cbind("SW/FW vs Coastal",Q7_Coastal_SW_FW2_chisq[[c(2,1)]],Q7_Coastal_SW_FW2_chisq$statistic,Q7_Coastal_SW_FW2_chisq$p.value))

Q7_chi_sq_SW_FW2 <- data.frame(rbind(Q7_New.Old_SW_FW,Q7_Age_SW_FW,Q7_Gender_SW_FW,Q7_Children_SW_FW,Q7_Race_Ethn_SW_FW,Q7_Coastal_SW_FW))
names(Q7_chi_sq_SW_FW2) <- c("Variable", "df","F","p")
Q7_chi_sq_SW_FW2$F <- as.numeric(Q7_chi_sq_SW_FW2$F)
Q7_chi_sq_SW_FW2$p <- as.numeric(Q7_chi_sq_SW_FW2$p)

# Proportions of anglers by waterbody type
SW_FW2 %>% 
  kbl(caption = "Saltwater or Freshwater Angler",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Saltwater or Freshwater Angler

SW_FW2	n	freq
FW	802	0.43
SW	448	0.24
SW.FW	631	0.34

# Significant chi-square test results for Fw/SW
Q7_New.Old_xtab

New.Old	SW_FW2			Total
	FW	SW	SW.FW
Old	180 9.4 %	167 8.7 %	236 12.3 %	583 30.4 %
New	575 30.1 %	264 13.8 %	491 25.7 %	1330 69.6 %
Total	755 39.5 %	431 22.5 %	727 38 %	1913 100 %

Q7_New.Old_SW_FW2_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~New.Old + SW_FW2, svy_desn_Q7, statistic = "adjWald")
## F = 10.583, ndf = 2, ddf = 2236, p-value = 2.662e-05

Q7_Gender_xtab

Gender	SW_FW2			Total
	FW	SW	SW.FW
Male	436 25.6 %	225 13.2 %	461 27 %	1122 65.8 %
Female	248 14.5 %	146 8.6 %	190 11.1 %	584 34.2 %
Total	684 40.1 %	371 21.7 %	651 38.2 %	1706 100 %

Q7_Gender_SW_FW2_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Gender + SW_FW2, svy_desn_Q7, statistic = "adjWald")
## F = 4.4308, ndf = 2, ddf = 2236, p-value = 0.01201

Q7_Age_xtab

Age	SW_FW2			Total
	FW	SW	SW.FW
45-54	115 6.7 %	78 4.6 %	106 6.2 %	299 17.5 %
18-24	127 7.4 %	47 2.8 %	144 8.4 %	318 18.6 %
25-34	157 9.2 %	71 4.2 %	174 10.2 %	402 23.6 %
35-44	162 9.5 %	90 5.3 %	156 9.1 %	408 23.9 %
55-64	81 4.7 %	60 3.5 %	56 3.3 %	197 11.5 %
65 or older	42 2.5 %	22 1.3 %	19 1.1 %	83 4.9 %
Total	684 40.1 %	368 21.6 %	655 38.4 %	1707 100 %

Q7_Age_SW_FW2_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Age + SW_FW2, svy_desn_Q7, statistic = "adjWald")
## F = 4.0948, ndf = 10, ddf = 2228, p-value = 1.302e-05

Q7_Children_xtab

Children	SW_FW2			Total
	FW	SW	SW.FW
18 or over	197 14.4 %	144 10.5 %	152 11.1 %	493 36 %
No children	245 17.9 %	110 8 %	287 21 %	642 46.9 %
Under 17	99 7.2 %	49 3.6 %	85 6.2 %	233 17 %
Total	541 39.5 %	303 22.1 %	524 38.3 %	1368 100 %

Q7_Children_SW_FW2_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Children + SW_FW2, svy_desn_Q7, statistic = "adjWald")
## F = 5.6924, ndf = 4, ddf = 2234, p-value = 0.0001478

Q7_Race_Ethn_xtab

Race_Ethn	SW_FW2			Total
	FW	SW	SW.FW
White	453 28.4 %	193 12.1 %	329 20.6 %	975 61.1 %
Asian or Pacific Islander	25 1.6 %	11 0.7 %	27 1.7 %	63 4 %
Black or African American	32 2 %	7 0.4 %	29 1.8 %	68 4.2 %
Hispanic	6 0.4 %	13 0.8 %	22 1.4 %	41 2.6 %
Native American or Alaskan Native	15 0.9 %	11 0.7 %	17 1.1 %	43 2.7 %
Other Hispanic	34 2.1 %	23 1.4 %	51 3.2 %	108 6.7 %
White Hispanic	80 5 %	93 5.8 %	124 7.8 %	297 18.6 %
Total	645 40.4 %	351 22 %	599 37.6 %	1595 100 %

Q7_Race_Ethn_SW_FW2_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Race_Ethn + SW_FW2, svy_desn_Q7, statistic = "adjWald")
## F = 3.7609, ndf = 12, ddf = 2226, p-value = 1.124e-05

Q7_Coastal_xtab

Coastal	SW_FW2			Total
	FW	SW	SW.FW
N	556 30.3 %	98 5.3 %	275 15 %	929 50.6 %
Y	164 9 %	319 17.4 %	420 22.9 %	903 49.3 %
Total	720 39.3 %	417 22.8 %	695 37.9 %	1832 100 %

Q7_Coastal_SW_FW2_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Coastal + SW_FW2, svy_desn_Q7, statistic = "adjWald")
## F = 135.15, ndf = 2, ddf = 2236, p-value < 2.2e-16

tab_df(Q7_chi_sq_SW_FW2)

Variable	df	F	p
SW/FW vs New/Retained	2	10.58	0.00
SW/FW vs Age	10	4.09	0.00
SW/FW vs Gender	2	4.43	0.01
SW/FW vs Children	4	5.69	0.00
SW/FW vs Race/Ethnicity	12	3.76	0.00
SW/FW vs Coastal	2	135.15	0.00

Anglers most commonly fish in freshwater (43%), followed by both freshwater and saltwater (34%), and finally in saltwater only (24%). This sample was stratified evenly among residents of coastal and inland counties, although the customer database when the survey was distributed reflected a split of 77% inland and 23% coastal counties. Anglers located in coastal counties were more likely to fish in saltwater, while those in inland counties were more likely to fish in freshwater. New anglers were more likely to fish in freshwater, while retained anglers were more likely to fish in saltwater. Females were more likely to fish in either saltwater or freshwater (but not both), while males were more likely to fish both in freshwater and saltwater. Anglers between 18-34 years were most likely to fish both in freshwater and saltwater, anglers between 45-65 were most likely to fish in saltwater, and anglers older than 65 years were most likely to fish in freshwater. Anglers with adult children were most likely to fish in saltwater, while anglers without children were most likely to fish in both saltwater and freshwater. White anglers were most likely to fish in freshwater, White Hispanic anglers were most likely to fish in saltwater, and Hispanic or Other Hispanic anglers were most likely to fish in both saltwater and freshwater.

Q8 From mid-March 2020 until today, did you go fishing from a boat?

# Relevel income variable
rev2_wts$Income <- relevel(rev2_wts$Income, ref = "Between $100,000 and $149,999")

# Revise data to evaluate whether fished from a boat and owned boat
rev2_wts <- rev2_wts %>%
  mutate(Boat1 = paste(Q8A,Q8B,Q8C,sep = "_"),
         Boat_Y_N = as.factor(case_when(grepl("Yes",Boat1) ~ "Yes",
                             grepl("No", Boat1) ~ "No")),
         Boat_Own = as.factor(case_when(grepl("owned",Boat1) ~ "Borrow",
                                        grepl("own", Boat1) ~ "Own")))

# Anglers that fished from a boat
responses_Q8_Y_N <- rev2_wts %>%
  filter(!(is.na(Boat_Y_N))) %>% # remove NA values
  count(Boat_Y_N) %>%
  mutate(freq = prop.table(n)) # 56% of anglers fished from a boat

# Anglers that owned the boat they fished from
responses_Q8_Own <- rev2_wts %>%
  filter(!(is.na(Boat_Own))) %>% # remove NA values 
  count(Boat_Own) %>%
  mutate(freq = prop.table(n)) # Of those anglers that fished from a boat, 36% owned the boat

# Run logistic regression
svy_desn_Q8 <- svydesign(ids = ~0,  weights = ~weights,data = rev2_wts) # Survey design
Q8_mod1 <- svyglm(Boat_Y_N ~ Gender + Age + Income + Children + Race_Ethn + SW_FW2 + COVID + New.Old + Coastal, design=svy_desn_Q8, family=quasibinomial)
summary(Q8_mod1) # examine model results
VIF(Q8_mod1) # collinearity
Q8_mod2 <- svyglm(Boat_Y_N  ~ Income + Race_Ethn + SW_FW2 + COVID + New.Old, design=svy_desn_Q8, family=quasibinomial)
summary(Q8_mod2)
VIF(Q8_mod2) # collinearity
car::Anova(Q8_mod2, type = 3, test.statistic = "F",error.df = degf(Q8_mod2$survey.design)) # Evaluate variable significance

# Estimate and odds ratio table
Q8_coef_table <- coef(summary(Q8_mod2))
Q8_AOR <- standardize_parameters(Q8_mod2, exp = TRUE) # Odds_ratio and 95% CI
Q8_summary_table <- cbind(Q8_coef_table,Odds_Ratio = Q8_AOR$Std_Odds_Ratio,CI_low = Q8_AOR$CI_low,CI_high = Q8_AOR$CI_high)

# Check model accuracy
logitgof(obs = Q8_mod2$y, fitted(Q8_mod2)) #run hoslem test to test goodness of fit
glm.probs <- predict(Q8_mod2,type = "response") #Predict() assigns fitted model values for each participants based on the model
glm.pred <- as.factor(ifelse(glm.probs > 0.5, 1, 0)) #If prediction is above 50% assign 1 and assign 0 if prediction is below 50%
prop.table(table(glm.pred == Q8_mod2$y))  # 64% True

# Table of responses
responses_Q8_Y_N %>% 
  kbl(caption = "Q8A. Did you fish from a boat?",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q8A. Did you fish from a boat?

Boat_Y_N	n	freq
No	823	0.44
Yes	1053	0.56

responses_Q8_Own %>% 
  kbl(caption = "Q8B. Did you own the boat you fished on?",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q8B. Did you own the boat you fished on?

Boat_Own	n	freq
Borrow	674	0.64
Own	379	0.36

# Model odds ratio results
Q8_summary_table %>% 
  kbl(caption = "Q8 Model Odds Ratios",digits = 3) %>%
  kable_classic(full_width = F, html_font = "Cambria")

Q8 Model Odds Ratios

	Estimate	Std. Error	t value	Pr(>|t|)	Odds_Ratio	CI_low	CI_high
(Intercept)	0.724	0.214	3.377	0.001	2.062	1.354	3.140
IncomeLess than $20,000	-1.173	0.379	-3.093	0.002	0.309	0.147	0.651
IncomeBetween $20,000 and $39,999	-0.997	0.303	-3.291	0.001	0.369	0.204	0.669
IncomeBetween $40,000 and $59,999	-0.611	0.262	-2.335	0.020	0.543	0.325	0.907
IncomeBetween $60,000 and $79,999	-0.461	0.240	-1.926	0.054	0.630	0.394	1.009
IncomeBetween $80,000 and $99,999	-0.179	0.248	-0.723	0.470	0.836	0.514	1.360
IncomeOver $150,000	0.040	0.235	0.172	0.864	1.041	0.657	1.651
Race_EthnAsian or Pacific Islander	-0.526	0.394	-1.337	0.182	0.591	0.273	1.279
Race_EthnBlack or African American	-1.625	0.353	-4.608	0.000	0.197	0.099	0.393
Race_EthnHispanic	-0.713	0.452	-1.577	0.115	0.490	0.202	1.190
Race_EthnNative American or Alaskan Native	0.681	0.474	1.437	0.151	1.976	0.780	5.009
Race_EthnOther Hispanic	-1.259	0.360	-3.497	0.000	0.284	0.140	0.575
Race_EthnWhite Hispanic	-0.490	0.209	-2.343	0.019	0.613	0.407	0.923
SW_FW2SW	0.649	0.189	3.427	0.001	1.914	1.320	2.775
SW_FW2SW.FW	0.787	0.177	4.444	0.000	2.197	1.552	3.110
COVIDYes	-0.378	0.178	-2.128	0.034	0.685	0.484	0.971
New.OldNew	-0.384	0.163	-2.359	0.018	0.681	0.495	0.937

# Barplot of responses
model.frame(Q8_mod2) %>% ggplot(aes(x = Boat_Y_N, fill = Income)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q8. Did you fish from a boat?")

model.frame(Q8_mod2) %>% ggplot(aes(x = Boat_Y_N, fill = Race_Ethn)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q8. Did you fish from a boat?")

model.frame(Q8_mod2) %>% ggplot(aes(x = Boat_Y_N, fill = SW_FW2)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q8. Did you fish from a boat?")

model.frame(Q8_mod2) %>% ggplot(aes(x = Boat_Y_N, fill = COVID)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q8. Did you fish from a boat?")

model.frame(Q8_mod2) %>% ggplot(aes(x = Boat_Y_N, fill = New.Old)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q8. Did you fish from a boat?")

# Effect plot
plot(Effect(focal.predictors = c("COVID","New.Old"),Q8_mod2))

Anglers fished from a boat 56% of the time and 36% of those that fished from a boat owned the boat.

The odds that that an angler fished from a boat was: 3.233 to 1.586 times less likely for anglers with an income <$20K to $79K compared to anglers with an income between $100K-$150K, 5.077 times less likely for Black or African American anglers, 3.523 times less likely for Other Hispanic anglers, and 1.632 times less likely for White Hispanic anglers compared to White anglers, 1.914 times more likely for anglers that fish in saltwater compared to those that fish in freshwater, 2.197 times more likely for anglers that fish in both saltwater and freshwater compared to those that fish in freshwater, 1.459 times less likely for anglers that were prompted by COVID-19 to start fishing, and 1.468 times more likely for new anglers compared to retained anglers.

Q9. The last time you went fishing, was fishing the main purpose of the trip, or was it just one of the activities but not the main purpose of the trip?

rev2_wts$Q9_Purpose <- recode_factor(rev2_wts$Q9_Purpose,"Fishing was one activity, but not the main purpose of the trip." = "No","Fishing was the main purpose of the trip." = "Yes")
svy_desn_Q9 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts) # Survey design
Q9_mod1 <- svyglm(Q9_Purpose ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age + Coastal, design=svy_desn_Q9, family=quasibinomial)
summary(Q9_mod1) # examine model results
varImp(Q9_mod1) # Variable importance
VIF(Q9_mod1) # collinearity
Q9_mod2 <- svyglm(Q9_Purpose ~ Gender + Age + Income + Race_Ethn + New.Old + Coastal, design=svy_desn_Q9, family=quasibinomial)
summary(Q9_mod2)
Q9_mod3 <- svyglm(Q9_Purpose ~ Gender + Age + Income + Race_Ethn + Coastal, design=svy_desn_Q9, family=quasibinomial)
summary(Q9_mod3)
VIF(Q9_mod3) # collinearity
car::Anova(Q9_mod3, type = 3, test.statistic = "F",error.df = degf(Q9_mod2$survey.design)) # Evaluate variable significance

# Check model accuracy
logitgof(obs = Q9_mod3$y, fitted(Q9_mod3)) #run hoslem test to test goodness of fit
glm.probs <- predict(Q9_mod3,type = "response") #Predict() assigns fitted model values for each participants based on the model
glm.pred <- as.factor(ifelse(glm.probs > 0.5, 1, 0)) #If prediction is above 50% assign 1 and assign 0 if prediction is below 50%
prop.table(table(glm.pred == Q9_mod3$y)) # Correct assignment 75% of the time 

# Estimate and odds ratio table
Q9_coef_table <- coef(summary(Q9_mod3))
Q9_AOR <- standardize_parameters(Q9_mod3, exp = TRUE) # Odds_ratio and 95% CI
Q9_summary_table <- cbind(Q9_coef_table,Odds_Ratio = Q9_AOR$Std_Odds_Ratio,CI_low = Q9_AOR$CI_low,CI_high = Q9_AOR$CI_high)

# Survey response proportions
Q9_prop <- rev2_wts %>%
  filter(!(is.na(Q9_Purpose))) %>% # remove NA values from Gender & Age
  count(Q9_Purpose) %>%
  mutate(freq = prop.table(n)) # get proportions

# Table of responses
Q9_prop %>% 
  kbl(caption = "Q9 Was Fishing the Main Purpose?",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") # 81% were already planning to go fishing

Q9 Was Fishing the Main Purpose?

Q9_Purpose	n	freq
No	468	0.25
Yes	1398	0.75

# Summary table
Q9_summary_table %>% 
  kbl(caption = "Q9 Model Odds Ratios",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q9 Model Odds Ratios

	Estimate	Std. Error	t value	Pr(>|t|)	Odds_Ratio	CI_low	CI_high
(Intercept)	0.93	0.22	4.15	0.00	2.53	1.63	3.92
GenderFemale	-0.34	0.16	-2.09	0.04	0.71	0.52	0.98
Age18-24	0.81	0.47	1.71	0.09	2.24	0.89	5.66
Age25-34	-0.21	0.24	-0.84	0.40	0.81	0.50	1.32
Age35-44	-0.16	0.21	-0.79	0.43	0.85	0.57	1.28
Age55-64	-0.10	0.22	-0.44	0.66	0.91	0.59	1.40
Age65 or older	-0.64	0.26	-2.44	0.01	0.53	0.31	0.88
IncomeLess than $20,000	-0.59	0.42	-1.39	0.16	0.55	0.24	1.28
IncomeBetween $20,000 and $39,999	1.23	0.36	3.37	0.00	3.42	1.67	7.00
IncomeBetween $40,000 and $59,999	0.54	0.29	1.85	0.06	1.71	0.97	3.02
IncomeBetween $60,000 and $79,999	0.39	0.27	1.45	0.15	1.48	0.87	2.52
IncomeBetween $80,000 and $99,999	0.10	0.26	0.38	0.71	1.10	0.66	1.83
IncomeOver $150,000	-0.27	0.24	-1.12	0.26	0.77	0.48	1.22
Race_EthnAsian or Pacific Islander	0.09	0.44	0.21	0.84	1.09	0.47	2.57
Race_EthnBlack or African American	0.79	0.42	1.91	0.06	2.21	0.98	5.00
Race_EthnHispanic	0.95	0.70	1.35	0.18	2.58	0.65	10.23
Race_EthnNative American or Alaskan Native	1.44	0.77	1.88	0.06	4.24	0.94	19.23
Race_EthnOther Hispanic	1.26	0.46	2.71	0.01	3.51	1.41	8.71
Race_EthnWhite Hispanic	0.34	0.25	1.33	0.18	1.40	0.85	2.29
CoastalY	0.41	0.17	2.48	0.01	1.51	1.09	2.09

# Histogram of responses
model.frame(Q9_mod3) %>% ggplot(aes(x = Q9_Purpose, fill = Age)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q9. Was Fishing the Main Purpose?")

model.frame(Q9_mod3) %>% ggplot(aes(x = Q9_Purpose, fill = Gender)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q9. Was Fishing the Main Purpose?")

model.frame(Q9_mod3) %>% ggplot(aes(x = Q9_Purpose, fill = Income)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q9. Was Fishing the Main Purpose?")

model.frame(Q9_mod3) %>% ggplot(aes(x = Q9_Purpose, fill = Race_Ethn)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q9. Was Fishing the Main Purpose?")

model.frame(Q9_mod3) %>% ggplot(aes(x = Q9_Purpose, fill = Coastal)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q9. Was Fishing the Main Purpose?")

# Effect plot
plot(Effect(focal.predictors = c("Gender","Age"),Q9_mod2))

plot(Effect(focal.predictors = c("Income"),Q9_mod2))

plot(Effect(focal.predictors = c("Coastal"),Q9_mod2))

Most anglers (75%) reported that fishing was the main purpose of their trip.

The odds that the main purpose of the trip was fishing was: 1.402 times less likely for female anglers compared to males, 2.241 times more likely for anglers ages 18-24 compared to anglers aged 45-54 yrs, 1.898 times less likely for anglers over 65 yrs compared to anglers aged 45-54 yrs, 3.422 times more likely for anglers with an income between $20K-$39K, 3.422 times more likely for anglers with an income between $20K-$39K, 3.509 times more likely for anglers that identified as ‘Other Hispanic’, 1.511 times more likely for anglers that were from a coastal county.

Q10. Who do you typically fish with?

# Evaluate frequencies by categorical responses
Q10_df1 <- rev2_wts %>% 
  dplyr::select(c("Q10A":"Q10J")) %>%
  pivot_longer(cols = "Q10A":"Q10J",names_to = "Q10_answers",values_to = "Agreement") %>%
  group_by(Q10_answers,Agreement) %>%
  count() %>%
  na.omit()
freq = prop.table(Q10_df1$n)
Q10_df1 <- data.frame(cbind(Q10_df1,freq=freq))

# Format data for analysis
Q10_df2 <- rev2_wts %>% 
    dplyr::select(Q10A,Q10F,Q10G,Q10H,Q10I,Q10J,Gender,Age,Race_Ethn,Income,SW_FW2,New.Old,Children,COVID,Coastal,weights) %>%
  pivot_longer(cols = "Q10A":"Q10J",names_to = "Q10_answers",values_to = "Agreement") %>%
  group_by(Q10_answers,Agreement,Gender,Age,Race_Ethn,Income,SW_FW2,New.Old,Children,COVID,Coastal,weights) %>%
  filter(!is.na(Agreement)) %>%
  mutate(Agreement = as.factor(Agreement)) %>%
  data.frame()

Q10_df2$Age <- ordered(Q10_df2$Age,levels = c("18-24", "25-34", "35-44","45-54","55-64","65 or older"))

Q10_df2$Income <- ordered(Q10_df2$Income,levels = c("Less than $20,000", "Between $20,000 and $39,999", "Between $40,000 and $59,999","Between $60,000 and $79,999","Between $80,000 and $99,999","Between $100,000 and $149,999","Over $150,000"))

# Chi-square tests
svy_desn_Q10 <- svydesign(ids = ~0, weights = ~weights, data = Q10_df2) # Survey design

Q10_Gender_chisq <- svychisq(~Gender + Agreement, svy_desn_Q10, statistic = "adjWald") # significant
Q10_Age_chisq <- svychisq(~Age + Agreement, svy_desn_Q10, statistic = "adjWald") # significant
Q10_Race_Ethn_chisq <- svychisq(~Race_Ethn + Agreement, svy_desn_Q10, statistic = "adjWald") # significant
Q10_Income_chisq <- svychisq(~Income + Agreement, svy_desn_Q10, statistic = "adjWald") # significant
Q10_SW_FW2_chisq <- svychisq(~SW_FW2 + Agreement, svy_desn_Q10, statistic = "adjWald") # significant
Q10_New.Old_chisq <- svychisq(~New.Old + Agreement, svy_desn_Q10, statistic = "adjWald") # not significant
Q10_Children_chisq <- svychisq(~Children + Agreement, svy_desn_Q10, statistic = "adjWald") # significant
Q10_COVID_chisq <- svychisq(~COVID + Agreement, svy_desn_Q10, statistic = "adjWald") # not very significant
Q10_Coastal_chisq <- svychisq(~Coastal + Agreement, svy_desn_Q10, statistic = "adjWald") # NOT significant

# Contingency tables for significant factors
Q10_Gender_xtab <- tab_xtab(Q10_df2$Gender,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q10_Age_xtab <- tab_xtab(Q10_df2$Age,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q10_Race_Ethn_xtab <- tab_xtab(Q10_df2$Race_Ethn,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q10_Income_xtab <- tab_xtab(Q10_df2$Income,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q10_SW_FW2_xtab <- tab_xtab(Q10_df2$SW_FW2,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q10_Children_xtab <- tab_xtab(Q10_df2$Children,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)

# Correlation plots of residuals to examine results
Q10_corr_Gender <- corrplot(Q10_Gender_chisq$residuals, is.cor = FALSE,method = 'number')

Q10_corr_Age <- corrplot(Q10_Age_chisq$residuals, is.cor = FALSE,method = 'number')

Q10_corr_Race_Ethn <- corrplot(Q10_Race_Ethn_chisq$residuals, is.cor = FALSE,method = 'number')

Q10_corr_Income <- corrplot(Q10_Income_chisq$residuals, is.cor = FALSE,method = 'number')

Q10_corr_SW_FW2 <- corrplot(Q10_SW_FW2_chisq$residuals, is.cor = FALSE,method = 'number')

Q10_corr_Children <- corrplot(Q10_Children_chisq$residuals, is.cor = FALSE,method = 'number')

# Combine chi-squared test results into dataframe
Q10_Gender <- data.frame(cbind("Gender",Q10_Gender_chisq[[c(2,1)]],Q10_Gender_chisq$statistic,Q10_Gender_chisq$p.value))
Q10_Age <- data.frame(cbind("Age",Q10_Age_chisq[[c(2,1)]],Q10_Age_chisq$statistic,Q10_Age_chisq$p.value))
Q10_Race_Ethn <- data.frame(cbind("Race/Ethnicity",Q10_Race_Ethn_chisq[[c(2,1)]],Q10_Race_Ethn_chisq$statistic,Q10_Race_Ethn_chisq$p.value))
Q10_Income <- data.frame(cbind("Income",Q10_Income_chisq[[c(2,1)]],Q10_Income_chisq$statistic,Q10_Income_chisq$p.value))
Q10_SW_FW2 <- data.frame(cbind("SW/FW",Q10_SW_FW2_chisq[[c(2,1)]],Q10_SW_FW2_chisq$statistic,Q10_SW_FW2_chisq$p.value))
Q10_Children <- data.frame(cbind("Children",Q10_Children_chisq[[c(2,1)]],Q10_Children_chisq$statistic,Q10_Children_chisq$p.value))

Q10_chi_sq <- data.frame(rbind(Q10_Age,Q10_Race_Ethn,Q10_Income,Q10_SW_FW2,Q10_Gender,Q10_Children))
names(Q10_chi_sq) <- c("Variable", "df","F","p")
Q10_chi_sq$F <- as.numeric(Q10_chi_sq$F)
Q10_chi_sq$p <- as.numeric(Q10_chi_sq$p)

# Response percentages graph
ggplot(Q10_df1,aes(x=Agreement,y=freq)) + geom_bar(stat="identity") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q10. Who do you go fishing with?") 

# Significant Chi-square test results
Q10_Gender_xtab

Gender	Agreement						Total
	No one, I go by myself	With my kids who are 12 years or younger	My family goes with extended family	With friends	With my teenagers, ages 13–19	With my grown children 20 years and older
Male	181 6.9 %	264 10.1 %	354 13.6 %	591 22.7 %	178 6.8 %	138 5.3 %	1706 65.4 %
Female	37 1.4 %	156 6 %	239 9.2 %	279 10.7 %	96 3.7 %	96 3.7 %	903 34.7 %
Total	218 8.4 %	420 16.1 %	593 22.7 %	870 33.3 %	274 10.5 %	234 9 %	2609 100 %

Q10_Gender_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Gender + Agreement, svy_desn_Q10, statistic = "adjWald")
## F = 7.2119, ndf = 5, ddf = 2854, p-value = 1.017e-06

Q10_Age_xtab 

Age	Agreement						Total
	No one, I go by myself	With my kids who are 12 years or younger	My family goes with extended family	With friends	With my teenagers, ages 13–19	With my grown children 20 years and older
18-24	47 1.8 %	9 0.3 %	121 4.6 %	229 8.8 %	11 0.4 %	23 0.9 %	440 16.8 %
25-34	56 2.1 %	124 4.8 %	144 5.5 %	242 9.3 %	23 0.9 %	6 0.2 %	595 22.8 %
35-44	40 1.5 %	183 7 %	143 5.5 %	169 6.5 %	117 4.5 %	28 1.1 %	680 26.1 %
45-54	37 1.4 %	70 2.7 %	87 3.3 %	131 5 %	92 3.5 %	67 2.6 %	484 18.5 %
55-64	30 1.1 %	24 0.9 %	67 2.6 %	70 2.7 %	23 0.9 %	69 2.6 %	283 10.8 %
65 or older	9 0.3 %	13 0.5 %	32 1.2 %	30 1.1 %	6 0.2 %	38 1.5 %	128 4.8 %
Total	219 8.4 %	423 16.2 %	594 22.8 %	871 33.4 %	272 10.4 %	231 8.9 %	2610 100 %

Q10_Age_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Age + Agreement, svy_desn_Q10, statistic = "adjWald")
## F = 21.671, ndf = 25, ddf = 2834, p-value < 2.2e-16

Q10_Race_Ethn_xtab

Race_Ethn	Agreement						Total
	No one, I go by myself	With my kids who are 12 years or younger	My family goes with extended family	With friends	With my teenagers, ages 13–19	With my grown children 20 years and older
White	133 5.5 %	233 9.6 %	340 14 %	466 19.2 %	156 6.4 %	139 5.7 %	1467 60.4 %
Asian or Pacific Islander	10 0.4 %	11 0.5 %	21 0.9 %	36 1.5 %	6 0.2 %	3 0.1 %	87 3.6 %
Black or African American	8 0.3 %	18 0.7 %	20 0.8 %	34 1.4 %	4 0.2 %	4 0.2 %	88 3.6 %
Hispanic	1 0 %	14 0.6 %	24 1 %	22 0.9 %	9 0.4 %	6 0.2 %	76 3.1 %
Native American or Alaskan Native	9 0.4 %	12 0.5 %	11 0.5 %	23 0.9 %	9 0.4 %	7 0.3 %	71 3 %
Other Hispanic	14 0.6 %	33 1.4 %	29 1.2 %	62 2.6 %	25 1 %	21 0.9 %	184 7.7 %
White Hispanic	31 1.3 %	60 2.5 %	106 4.4 %	173 7.1 %	41 1.7 %	38 1.6 %	449 18.6 %
Total	206 8.5 %	381 15.7 %	551 22.7 %	816 33.7 %	250 10.3 %	218 9 %	2422 100 %

Q10_Race_Ethn_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Race_Ethn + Agreement, svy_desn_Q10, statistic = "adjWald")
## F = 1.9732, ndf = 30, ddf = 2829, p-value = 0.001257

Q10_Income_xtab 

Income	Agreement						Total
	No one, I go by myself	With my kids who are 12 years or younger	My family goes with extended family	With friends	With my teenagers, ages 13–19	With my grown children 20 years and older
Less than $20,000	22 1.1 %	11 0.5 %	39 1.9 %	36 1.8 %	7 0.3 %	4 0.2 %	119 5.8 %
Between $20,000 and $39,999	25 1.2 %	41 2 %	57 2.8 %	91 4.5 %	15 0.7 %	14 0.7 %	243 11.9 %
Between $40,000 and $59,999	33 1.6 %	41 2 %	62 3.1 %	87 4.3 %	22 1.1 %	17 0.8 %	262 12.9 %
Between $60,000 and $79,999	31 1.5 %	61 3 %	68 3.4 %	117 5.8 %	28 1.4 %	27 1.3 %	332 16.4 %
Between $80,000 and $99,999	26 1.3 %	50 2.5 %	59 2.9 %	85 4.2 %	35 1.7 %	30 1.5 %	285 14.1 %
Between $100,000 and $149,999	28 1.4 %	69 3.4 %	93 4.6 %	138 6.8 %	59 2.9 %	44 2.2 %	431 21.3 %
Over $150,000	26 1.3 %	64 3.2 %	70 3.5 %	103 5.1 %	49 2.4 %	37 1.8 %	349 17.3 %
Total	191 9.5 %	337 16.7 %	448 22.2 %	657 32.5 %	215 10.6 %	173 8.6 %	2021 100 %

Q10_Income_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Income + Agreement, svy_desn_Q10, statistic = "adjWald")
## F = 1.5106, ndf = 30, ddf = 2829, p-value = 0.03712

Q10_SW_FW2_xtab 

SW_FW2	Agreement						Total
	No one, I go by myself	With my kids who are 12 years or younger	My family goes with extended family	With friends	With my teenagers, ages 13–19	With my grown children 20 years and older
FW	124 4.3 %	204 7.1 %	239 8.4 %	323 11.3 %	100 3.5 %	87 3 %	1077 37.6 %
SW	29 1 %	90 3.2 %	146 5.1 %	223 7.8 %	72 2.5 %	73 2.6 %	633 22.2 %
SW.FW	78 2.7 %	175 6.1 %	258 9 %	401 14 %	129 4.5 %	104 3.6 %	1145 39.9 %
Total	231 8.1 %	469 16.4 %	643 22.5 %	947 33.2 %	301 10.5 %	264 9.2 %	2855 100 %

Q10_SW_FW2_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~SW_FW2 + Agreement, svy_desn_Q10, statistic = "adjWald")
## F = 3.3766, ndf = 10, ddf = 2849, p-value = 0.0002146

Q10_Children_xtab 

Children	Agreement						Total
	No one, I go by myself	With my kids who are 12 years or younger	My family goes with extended family	With friends	With my teenagers, ages 13–19	With my grown children 20 years and older
18 or over	47 2.3 %	81 4 %	182 9.1 %	186 9.3 %	129 6.4 %	179 8.9 %	804 40 %
No children	121 6 %	13 0.6 %	210 10.5 %	427 21.3 %	19 0.9 %	32 1.6 %	822 40.9 %
Under 17	25 1.2 %	115 5.7 %	83 4.1 %	118 5.9 %	27 1.3 %	15 0.7 %	383 18.9 %
Total	193 9.6 %	209 10.4 %	475 23.6 %	731 36.4 %	175 8.7 %	226 11.2 %	2009 100 %

Q10_Children_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Children + Agreement, svy_desn_Q10, statistic = "adjWald")
## F = 37.241, ndf = 10, ddf = 2849, p-value < 2.2e-16

Anglers typically fish with friends, extended family, or with kids younger than 12 years. Male anglers are more likely to fish by themselves compared to female anglers. Anglers 45+ are more likely to fish with grown children than younger anglers, anglers between the ages of 35-54 are most likely to fish with young children and teenagers, while anglers 18-24 are most likely to fish with friends. Asian or Pacific Islander anglers and White Hispanic anglers are most likely to fish with friends, Black or African American anglers are most likely to fish with young children, Hispanic anglers are most likely to fish with extended family, Native American or Alaskan Native anglers are most likely to fish by themselves, Other Hispanic anglers are most likely to fish with teenagers, and White anglers are most likely to fish with grown children. Anglers with incomes over $80K are more likely to fish with teenagers and grown children, anglers with incomes between <$20-$80K are more likely to fish with friends, extended family, or by themselves. Freshwater anglers are most likely to fish by themselves, saltwater anglers are most likely to fish with grown children, and anglers that fish in both freshwater and saltwater are most likely to fish with friends. Anglers prompted by COVID to fish were most likely to fish by themselves, while those that were already planning to fish were most likely to fish with grown children. Anglers with adult children were most likely to fish with adult children, anglers with young children were most likely to fish with young children, and anglers with no children were most likely to fish with friends.

Q11. Please indicate how much you agree with each of the following statements about fishing.

# Clean Q11 survey answers
rev2_wts <- rev2_wts %>%  
  mutate_at(vars(Q11A:Q11D), list(~recode_factor(.,"Strongly Disagree" = "-2","Disagree" = "-1","Agree" = "1","Strongly Agree"= "2"))) %>%
  mutate(across(Q11A:Q11D, ~as.numeric(as.character(.)))) %>%
  mutate_at(vars(Q11A:Q11D), list(~.*weights)) %>%
  mutate(Enjoyed_fishing = Q11A,
         Too_expensive = Q11B,
         Fishing_rules_confusing = Q11C,
         Enjoyed_surroundings = Q11D)

# Remove rows with missing independent variable data
Q11_df1 <- rev2_wts %>%
  filter_at(vars(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2,Coastal),all_vars(!is.na(.))) 
# Select independent variables
Q11_indep1 <- Q11_df1 %>%
  dplyr::select(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2,Coastal)
# Select dependent variables
Q11_dep1 <- Q11_df1 %>%
  dplyr::select(Enjoyed_fishing,Too_expensive,Fishing_rules_confusing,Enjoyed_surroundings)
# Transform dependent data
Q11_dep2 <- decostand(Q11_dep1, method = 'standardize') 
Q11_indep2 <- dudi.mix(Q11_indep1, scannf=F, nf=2)$tab # MCA, keep 2 axes

# Run RDA
Q11_rda_mod1 <- rda(Q11_dep2 ~ ., data = Q11_indep2) 
envfit(Q11_rda_mod1,Q11_indep2) # select variables at alpha = 0.001 with R2>0.05
Q11_indep3 <- Q11_indep2 %>%
  dplyr::select(Gende.Female,Gende.Male,Incom.Between..20.000.and..39.999,Incom.Over..150.000, Age.18.24, Age.25.34, Age.45.54, Age.55.64, Age.65.or.older, Child.18.or.over, Child.No.children, Child.Under.17,Race_.White,Race_.White.Hispanic,Race_.Other.Hispanic, SW_FW.SW,SW_FW.SW.FW,New.O.New,New.O.Old,Coast.N,Coast.Y)
Q11_rda_mod2 <- rda(Q11_dep2 ~ ., data = Q11_indep3)
envfit(Q11_rda_mod2,Q11_indep3) # select variables at alpha = 0.001 for next model
Q11_indep4 <- Q11_indep2 %>%
  dplyr::select(Age.18.24, Age.25.34, Age.55.64, Age.65.or.older, Child.18.or.over, Child.No.children)
Q11_rda_mod3 <- rda(Q11_dep2 ~ ., data = Q11_indep4)
envfit(Q11_rda_mod3,Q11_indep4) # select variables at alpha = 0.001 for next model
Q11_rda_mod3$call
RsquareAdj(Q11_rda_mod3) 
sqrt(vif.cca(Q11_rda_mod3))
Q11_rda_mod3

# Extract scores for visualization
# Basic plot for visualization
vare_spp_sco <- scores(Q11_rda_mod3, display = "species")
vare_sam_sco <- scores(Q11_rda_mod3, display = "sites")
vare_env_sco <- scores(Q11_rda_mod3, display = "bp")
vare_spp_tbl <- as_tibble(vare_spp_sco)
vare_sam_tbl <- as_tibble(vare_sam_sco)
vare_env_tbl <- as_tibble(vare_env_sco)
vare_spp_tbl <- mutate(vare_spp_tbl, vgntxt=c("Enjoyed fishing","Too expensive","Fishing rules confusing", "Enjoyed surroundings"), ccatype = "species")
vare_sam_tbl <- mutate(vare_sam_tbl, vgntxt=rownames(vare_sam_sco), ccatype = "sites")
vare_env_tbl <- mutate(vare_env_tbl, vgntxt=rownames(vare_env_sco), ccatype = "bp")

vare_env_tbl <- mutate(vare_env_tbl, vgntxt=c("Age 18-24","Age 25-34","Age 55-64","Age 65 over","Children 18 over","No Children"),ccatype = "bp")
rescaled <- vare_env_tbl %>% 
            dplyr::select(RDA1, RDA2) %>%
            as.matrix() *2.5
vare_tbl <- dplyr::select(vare_env_tbl, vgntxt, ccatype) %>%
            bind_cols(as_tibble(rescaled)) %>%
            bind_rows(vare_spp_tbl)

# Get mean, SD, and calculate cronbach's alpha
# Recode factors 1-5 for comparison to other questions in final table
Q11_summary1 <- rev1_wts %>%  
  mutate_at(vars(Q11A:Q11D), list(~recode_factor(.,"Strongly Disagree" = "1","Disagree" = "2","Agree" = "3","Strongly Agree"= "4"))) %>%
  mutate(across(Q11A:Q11D, ~as.numeric(as.character(.)))) %>%
  mutate_at(vars(Q11A:Q11D), list(~.*weights)) %>%
  mutate(Enjoyed_fishing = Q11A,
         Too_expensive = Q11B,
         Fishing_rules_confusing = Q11C,
         Enjoyed_surroundings = Q11D) %>%
  dplyr::select(Enjoyed_fishing,Too_expensive,Fishing_rules_confusing,Enjoyed_surroundings)

Q11_summary2 <- describe(Q11_summary1,na.rm = TRUE,IQR = TRUE) 
Q11_summary2$Variables <- row.names(Q11_summary2) 
Q11_summary3 <- Q11_summary2[,c("Variables","median","IQR")]
Q11_summary3$Variables<-str_replace_all(Q11_summary3$Variables, c("_" = " ")) # replaces spaces in names with periods

# Calculate Cronbach's alpha
Q11_cronbach.alpha <- cronbach.alpha(abs(Q11_summary1),na.rm = TRUE,standardized = TRUE)
Q11_cronbach.alpha$alpha

# Pull out significant variable terms for each axis
Q11_envfit1 <- envfit(Q11_rda_mod3,Q11_indep4) # select variables at alpha = 0.001 for next model
Q11_envfit2 <- data.frame(Q11_envfit1$vectors[1],Q11_envfit1$vectors[2],Q11_envfit1$vectors[4])
names(Q11_envfit2) <- c("RDA 1","RDA 2","r2","p")
Q11_envfit2$Variables <- row.names(Q11_envfit2) 
Q11_envfit2 <- Q11_envfit2 %>% relocate(Variables, .before = "RDA 1")

# Descriptive stats table
tab_df(Q11_summary3)

Variables	median	IQR
Enjoyed fishing	3.25	1.98
Too expensive	1.67	1.41
Fishing rules confusing	1.63	1.23
Enjoyed surroundings	3.15	1.99

# Reorganize data to plot
Q11_a <- rev1_wts[c(which(colnames(rev1_wts)=="Q11A"):which(colnames(rev1_wts)=="Q11D"))]
colnames(Q11_a) <- c("Enjoyed_fishing","Too_expensive","Fishing_rules_confusing","Enjoyed_surroundings")
Q11_b <- Q11_a %>% 
  pivot_longer(cols = Enjoyed_fishing:Enjoyed_surroundings,names_to = "Q11_answers",values_to = "Agreement")
Q11_c <- Q11_b %>%
  group_by(Q11_answers,Agreement) %>%
  count() %>%
  na.omit()
Q11_c$Agreement <- ordered(Q11_c$Agreement,levels = c("Strongly Agree", "Agree", "Disagree","Strongly Disagree") )
ggplot(Q11_c, aes(fill=Agreement, y=n, x=Q11_answers)) + geom_bar(position="fill", stat="identity") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8),axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + labs(y = "Percent", x = "Q11. Agreement on Fishing Statements") 

# PERMANOVA
anova.cca(Q11_rda_mod3, step = 1000, by = "axis",permutations = 999)

## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = Q11_dep2 ~ Age.18.24 + Age.25.34 + Age.55.64 + Age.65.or.older + Child.18.or.over + Child.No.children, data = Q11_indep4)
##           Df Variance        F Pr(>F)    
## RDA1       1  1.84134 785.7323  0.001 ***
## RDA2       1  0.00266   1.1330  0.981    
## RDA3       1  0.00183   0.7801  0.987    
## RDA4       1  0.00052   0.2239  1.000    
## Residual 919  2.15365                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_df(Q11_envfit2)

Variables	RDA.1	RDA.2	r2	p
Age.18.24	1.00	0.01	0.54	0.00
Age.25.34	1.00	0.01	0.08	0.00
Age.55.64	-1.00	-0.01	0.09	0.00
Age.65.or.older	-1.00	-0.01	0.05	0.00
Child.18.or.over	-1.00	-0.01	0.21	0.00
Child.No.children	1.00	0.02	0.12	0.00

# Basic visualization
ggplot() +
  geom_point(aes(x=RDA1, y=RDA2), data=filter(vare_tbl, ccatype=="species"))  +
  geom_text_repel(aes(x=RDA1, y=RDA2, label=vgntxt, size=3,colour=ccatype),
                  data=vare_tbl, seed=123) + 
  geom_segment(aes(x=0, y=0, xend=RDA1, yend=RDA2), arrow=arrow(length = unit(0.2,"cm")),
               data=filter(vare_tbl, ccatype=="bp"), color="blue") +
  coord_fixed() +
  theme_classic() +
  scale_colour_manual(values = c("blue", "black")) +
  ylim(-1,1) +
  theme(legend.position="none")

Most individuals agreed with the statements that they enjoyed fishing and their surroundings and disagreed with the statements that fishing rules were confusing or fishing was too expensive.

Anglers without children between the ages of 18-34 were more likely to report that they enjoy their surroundings and enjoy fishing, while those with adult children over the age of 55 were more likely to report that the fishing rules were confusing and that fishing was too expensive.

Cronbach’s alpha for Q11 is 0.9660789

Q12. How often did you catch fish during your trips?

# Order the factor
rev2_wts$Q12_CatchOften <- as.factor(rev2_wts$Q12_CatchOften)
rev2_wts$Q12_CatchOften  <- ordered(rev2_wts$Q12_CatchOften, levels = c("Never","Rarely", "Sometimes","Most of the time","Always"))
svy_desn_Q12 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts,na.rm = TRUE) # New survey design

# Ordinal regression
Q12_mod1 <- svyolr(Q12_CatchOften ~  SW_FW2 + Age + Gender + Income + Children + New.Old + COVID + Race_Ethn + Coastal, design=svy_desn_Q12)
summary(Q12_mod1) # examine model results
VIF(Q12_mod1) # collinearity
Q12_mod2 <- svyolr(Q12_CatchOften ~ SW_FW2 + Income + COVID + Race_Ethn, design=svy_desn_Q12)
summary(Q12_mod2) # examine model results
VIF(Q12_mod2) # collinearity
Q12_mod3 <- svyolr(Q12_CatchOften ~ SW_FW2 + COVID + Race_Ethn, design=svy_desn_Q12)
summary(Q12_mod3) # examine model results
VIF(Q12_mod3) # collinearity

# Estimate and odds ratio table
Q12_coef_table <- coef(summary(Q12_mod3))
Q12_pval <- pnorm(abs(Q12_coef_table[, "t value"]),lower.tail = FALSE)* 2
Q12_AOR <- exp(coef(Q12_mod3))
Q12_summary_table <- cbind(Q12_coef_table, "p value" = round(Q12_pval,3),Odds_Ratio = Q12_AOR)

# Survey response proportions
Q12_prop <- rev2_wts %>%
  filter(!(is.na(Q12_CatchOften))) %>% # remove NA values from Gender & Age
  count(Q12_CatchOften) %>%
  mutate(freq = prop.table(n)) # get proportions

# Table of responses
Q12_prop %>% 
  kbl(caption = "Q12 How often do you catch fish?",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q12 How often do you catch fish?

Q12_CatchOften	n	freq
Never	38	0.02
Rarely	181	0.10
Sometimes	541	0.29
Most of the time	829	0.45
Always	254	0.14

# Odds Ratio summary table
Q12_summary_table %>% 
  kbl(caption = "Q12 Model Odds Ratios",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q12 Model Odds Ratios

	Value	Std. Error	t value	p value	Odds_Ratio
SW_FW2SW	0.18	0.16	1.10	0.27	1.19
SW_FW2SW.FW	0.31	0.14	2.15	0.03	1.36
COVIDYes	-0.59	0.16	-3.61	0.00	0.56
Race_EthnAsian or Pacific Islander	-0.38	0.26	-1.46	0.14	0.68
Race_EthnBlack or African American	-0.72	0.24	-2.98	0.00	0.49
Race_EthnHispanic	-1.31	0.67	-1.97	0.05	0.27
Race_EthnNative American or Alaskan Native	0.46	0.16	2.92	0.00	1.59
Race_EthnOther Hispanic	-0.06	0.32	-0.19	0.85	0.94
Race_EthnWhite Hispanic	-0.77	0.17	-4.45	0.00	0.46
Never|Rarely	-3.86	0.28	-13.91	0.00	0.02
Rarely|Sometimes	-2.30	0.14	-16.60	0.00	0.10
Sometimes|Most of the time	-0.52	0.11	-4.74	0.00	0.59
Most of the time|Always	1.80	0.13	13.49	0.00	6.03

# Histogram of responses
model.frame(Q12_mod3) %>% ggplot(aes(x = Q12_CatchOften, fill = COVID)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q12. How often do you catch fish?")

model.frame(Q12_mod3) %>% ggplot(aes(x = Q12_CatchOften, fill = SW_FW2)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q12. How often do you catch fish?")

model.frame(Q12_mod3) %>% ggplot(aes(x = Q12_CatchOften, fill = Race_Ethn)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q12. How often do you catch fish?")

Most anglers (88%) reported that they usually catch fish (either sometimes, most of the time, or always), while only 12% of anglers rarely (10%) or never (2%) catch fish.

The odds of catching fish most of the time were: 1.362 times more likely for for anglers that fish in both freshwater and saltwater compared to just freshwater anglers, 1.795 times less likely for anglers that were prompted by COVID-19 to fish, 2.054 times less likely for Black or African American anglers, 3.712 times less likely for Hispanic anglers, 1.59 times more likely for Native American or Alaskan Native anglers, and 2.151 times less likely for White Hispanic anglers.

Q13. How often did you keep the fish you caught?

# Order the factor
rev2_wts$Q13_KeepOften <- as.factor(rev2_wts$Q13_KeepOften)
rev2_wts$Q13_KeepOften  <- ordered(rev2_wts$Q13_KeepOften, levels = c("Never","Rarely", "Sometimes", "Most of the time","Always"))
svy_desn_Q13 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts,na.rm = TRUE) # New survey design

# Ordinal regression
Q13_mod1 <- svyolr(Q13_KeepOften ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age + Coastal, design=svy_desn_Q13)
summary(Q13_mod1) # examine model results
VIF(Q13_mod1) # collinearity
Q13_mod2 <- svyolr(Q13_KeepOften ~ New.Old + COVID + SW_FW2 + Age + Race_Ethn, design=svy_desn_Q13)
summary(Q13_mod2) # examine model results
VIF(Q13_mod2) # collinearity

# Create odds ratio summary table
Q13_coef_table <- coef(summary(Q13_mod2))
Q13_pval <- pnorm(abs(Q13_coef_table[, "t value"]),lower.tail = FALSE)* 2
Q13_AOR <- exp(coef(Q13_mod2))
Q13_summary_table <- cbind(Q13_coef_table, "p value" = round(Q13_pval,3),Odds_Ratio = Q13_AOR)

# Survey response proportions
Q13_prop <- rev2_wts %>%
  filter(!(is.na(Q13_KeepOften))) %>% # remove NA values from Gender & Age
  count(Q13_KeepOften) %>%
  mutate(freq = prop.table(n)) # get proportions

# Table of responses
Q13_prop %>% 
  kbl(caption = "Q13 How often do you keep the fish you catch?",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q13 How often do you keep the fish you catch?

Q13_KeepOften	n	freq
Never	429	0.24
Rarely	409	0.23
Sometimes	572	0.32
Most of the time	296	0.16
Always	97	0.05

# Odds Ratio summary table
Q13_summary_table %>% 
  kbl(caption = "Q13 Model Odds Ratios",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q13 Model Odds Ratios

	Value	Std. Error	t value	p value	Odds_Ratio
New.OldNew	-0.37	0.12	-3.01	0.00	0.69
COVIDYes	-0.46	0.15	-3.14	0.00	0.63
SW_FW2SW	1.24	0.16	7.96	0.00	3.45
SW_FW2SW.FW	0.79	0.14	5.56	0.00	2.20
Age18-24	-0.17	0.26	-0.66	0.51	0.84
Age25-34	-0.07	0.16	-0.42	0.68	0.94
Age35-44	0.18	0.14	1.25	0.21	1.19
Age55-64	0.35	0.15	2.42	0.02	1.42
Age65 or older	0.80	0.18	4.49	0.00	2.24
Race_EthnAsian or Pacific Islander	0.44	0.25	1.74	0.08	1.56
Race_EthnBlack or African American	0.93	0.43	2.16	0.03	2.52
Race_EthnHispanic	-0.01	0.27	-0.05	0.96	0.99
Race_EthnNative American or Alaskan Native	0.37	0.35	1.07	0.28	1.45
Race_EthnOther Hispanic	-0.17	0.27	-0.62	0.54	0.85
Race_EthnWhite Hispanic	-0.08	0.16	-0.49	0.62	0.93
Never|Rarely	-0.85	0.15	-5.53	0.00	0.43
Rarely|Sometimes	0.36	0.15	2.37	0.02	1.43
Sometimes|Most of the time	1.94	0.17	11.47	0.00	6.93
Most of the time|Always	3.65	0.20	18.26	0.00	38.59

# Histogram of responses
model.frame(Q13_mod2) %>% ggplot(aes(x = Q13_KeepOften, fill = COVID)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q13. How often do you keep fish?")

model.frame(Q13_mod2) %>% ggplot(aes(x = Q13_KeepOften, fill = New.Old)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q13. How often do you keep fish?")

model.frame(Q13_mod2) %>% ggplot(aes(x = Q13_KeepOften, fill = SW_FW2)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q13. How often do you keep fish?")

model.frame(Q13_mod2) %>% ggplot(aes(x = Q13_KeepOften, fill = Age)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q13. How often do you keep fish?")

model.frame(Q13_mod2) %>% ggplot(aes(x = Q13_KeepOften, fill = Race_Ethn)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q13. How often do you keep fish?")

Approximately half of anglers (47%) reported that they do not usually keep the fish they catch (either rarely or never), while 32% of anglers sometimes keep their catch and only 21% of anglers usually keep their catch (most of the time or always).

The odds of keeping the fish that anglers caught were: 1.447 times less likely for for new anglers recruited during the pandemic, 1.577 times less likely for anglers that were prompted by COVID-19 to start fishing, 3.446 times more likely for saltwater anglers, 2.198 times more likely for anglers that fish in both freshwater and saltwater, 1.423to 2.236 times more likely for anglers 55+, and 2.523 times more likely for Black or African American anglers.

Q14. When you went fishing, was there anything you did not enjoy? Choose all that apply.

# Evaluate frequencies by categorical responses
rev2_wts$Income <- factor(rev2_wts$Income,levels = c("Less than $20,000", "Between $20,000 and $39,999", "Between $40,000 and $59,999","Between $60,000 and $79,999","Between $80,000 and $99,999","Between $100,000 and $149,999","Over $150,000"))
Q14_df1 <- rev2_wts %>% 
  dplyr::select(c("Q14A":"Q14F")) %>%
  pivot_longer(cols = "Q14A":"Q14F",names_to = "Q14_answers",values_to = "Agreement") %>%
  group_by(Q14_answers,Agreement) %>%
  count() %>%
  na.omit()
freq = prop.table(Q14_df1$n)
Q14_df1 <- data.frame(cbind(Q14_df1,freq=freq))
Q14_df1$Agreement <- recode_factor(Q14_df1$Agreement,"No, I enjoyed everything." = "Enjoyed everything","It was tough to find a good fishing spot close to home." = "Finding close fishing spot","I didnt have all of the needed equipment to fish." = "Needed equipment","I wasnt sure about the right techniques for casting or reeling in my bait/lure." = "Techniques for casting/reeling","My children had a lot of trouble fishing (casting, paying attention, avoiding anything dangerous, etc.)." = "Children had trouble","I didnt catch anything." = "Didn't catch anything")

# Format data for analysis
Q14_df2 <- rev2_wts %>% 
    dplyr::select(Q14A,Q14B,Q14C,Q14D,Q14E,Q14F,Gender,Age,Race_Ethn,Income,SW_FW2,New.Old,Children,COVID,Coastal,weights) %>%
  pivot_longer(cols = "Q14A":"Q14F",names_to = "Q14_answers",values_to = "Agreement") %>%
  group_by(Q14_answers,Agreement,Gender,Age,Race_Ethn,Income,SW_FW2,New.Old,Children,COVID,Coastal,weights) %>%
  filter(!is.na(Agreement)) %>%
  mutate(Agreement = as.factor(Agreement)) %>%
  data.frame()
Q14_df2$Agreement <- recode_factor(Q14_df2$Agreement,"No, I enjoyed everything." = "Enjoyed everything","It was tough to find a good fishing spot close to home." = "Finding close fishing spot","I didnt have all of the needed equipment to fish." = "Needed equipment","I wasnt sure about the right techniques for casting or reeling in my bait/lure." = "Techniques for casting/reeling","My children had a lot of trouble fishing (casting, paying attention, avoiding anything dangerous, etc.)." = "Children had trouble","I didnt catch anything." = "Didn't catch anything")

# Reorder Age and Income Factors
Q14_df2$Age <- ordered(Q14_df2$Age,levels = c("18-24", "25-34", "35-44","45-54","55-64","65 or older"))

Q14_df2$Income <- ordered(Q14_df2$Income,levels = c("Less than $20,000", "Between $20,000 and $39,999", "Between $40,000 and $59,999","Between $60,000 and $79,999","Between $80,000 and $99,999","Between $100,000 and $149,999","Over $150,000"))

# Chi-square tests
svy_desn_Q14 <- svydesign(ids = ~0, weights = ~weights, data = Q14_df2) # Survey design

Q14_Gender_chisq <- svychisq(~Gender + Agreement, svy_desn_Q14, statistic = "adjWald") # significant
Q14_Age_chisq <- svychisq(~Age + Agreement, svy_desn_Q14, statistic = "adjWald") # significant
Q14_Race_Ethn_chisq <- svychisq(~Race_Ethn + Agreement, svy_desn_Q14, statistic = "adjWald") # significant
Q14_Income_chisq <- svychisq(~Income + Agreement, svy_desn_Q14, statistic = "adjWald") # NOT significant
Q14_SW_FW2_chisq <- svychisq(~SW_FW2 + Agreement, svy_desn_Q14, statistic = "adjWald") # significant
Q14_New.Old_chisq <- svychisq(~New.Old + Agreement, svy_desn_Q14, statistic = "adjWald") # significant
Q14_Children_chisq <- svychisq(~Children + Agreement, svy_desn_Q14, statistic = "adjWald") # significant
Q14_COVID_chisq <- svychisq(~COVID + Agreement, svy_desn_Q14, statistic = "adjWald") # significant
Q14_Coastal_chisq <- svychisq(~Coastal + Agreement, svy_desn_Q14, statistic = "adjWald") # NOT significant

# Correlation plots of residuals to examine results
Q14_corr_Gender <- corrplot(Q14_Gender_chisq$residuals, is.cor = FALSE,method = 'number')

Q14_corr_Age <- corrplot(Q14_Age_chisq$residuals, is.cor = FALSE,method = 'number')

Q14_corr_Race_Ethn <- corrplot(Q14_Race_Ethn_chisq$residuals, is.cor = FALSE,method = 'number')

Q14_corr_SW_FW2 <- corrplot(Q14_SW_FW2_chisq$residuals, is.cor = FALSE,method = 'number')

Q14_corr_Children <- corrplot(Q14_Children_chisq$residuals, is.cor = FALSE,method = 'number')

Q14_corr_COVID <- corrplot(Q14_COVID_chisq$residuals, is.cor = FALSE,method = 'number')

Q14_corr_New.Old <- corrplot(Q14_New.Old_chisq$residuals, is.cor = FALSE,method = 'number')

# Contingency tables for significant factors
Q14_Gender_xtab <- tab_xtab(Q14_df2$Gender,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q14_Age_xtab <- tab_xtab(Q14_df2$Age,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q14_Race_Ethn_xtab <- tab_xtab(Q14_df2$Race_Ethn,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q14_SW_FW2_xtab <- tab_xtab(Q14_df2$SW_FW2,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q14_Children_xtab <- tab_xtab(Q14_df2$Children,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q14_COVID_xtab <- tab_xtab(Q14_df2$COVID,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q14_New.Old_xtab <- tab_xtab(Q14_df2$New.Old,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)

# Combine chi-squared test results into dataframe
Q14_Gender <- data.frame(cbind("Gender",Q14_Gender_chisq[[c(2,1)]],Q14_Gender_chisq$statistic,Q14_Gender_chisq$p.value))
Q14_Age <- data.frame(cbind("Age",Q14_Age_chisq[[c(2,1)]],Q14_Age_chisq$statistic,Q14_Age_chisq$p.value))
Q14_Race_Ethn <- data.frame(cbind("Race/Ethnicity",Q14_Race_Ethn_chisq[[c(2,1)]],Q14_Race_Ethn_chisq$statistic,Q14_Race_Ethn_chisq$p.value))
Q14_New.Old <- data.frame(cbind("New/Retained",Q14_New.Old_chisq[[c(2,1)]],Q14_New.Old_chisq$statistic,Q14_New.Old_chisq$p.value))
Q14_SW_FW2 <- data.frame(cbind("SW/FW",Q14_SW_FW2_chisq[[c(2,1)]],Q14_SW_FW2_chisq$statistic,Q14_SW_FW2_chisq$p.value))
Q14_COVID <- data.frame(cbind("COVID-Prompted",Q14_COVID_chisq[[c(2,1)]],Q14_COVID_chisq$statistic,Q14_COVID_chisq$p.value))
Q14_Children <- data.frame(cbind("Children",Q14_Children_chisq[[c(2,1)]],Q14_Children_chisq$statistic,Q14_Children_chisq$p.value))

Q14_chi_sq <- data.frame(rbind(Q14_COVID,Q14_Age,Q14_Race_Ethn,Q14_New.Old,Q14_SW_FW2,Q14_Gender,Q14_Children))
names(Q14_chi_sq) <- c("Variable", "df","F","p")
Q14_chi_sq$F <- as.numeric(Q14_chi_sq$F)
Q14_chi_sq$p <- as.numeric(Q14_chi_sq$p)

# Response percentages graph
ggplot(Q14_df1,aes(x=Agreement,y=freq)) + geom_bar(stat="identity") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q14. Was there anything about fishing you didn't enjoy?") 

tab_df(Q14_chi_sq,digits = 3,title = "Q14. Was there anything about fishing you didn't enjoy?")

Q14. Was there anything about fishing you didn’t enjoy?

Variable	df	F	p
COVID-Prompted	5	13.978	0.000
Age	25	4.910	0.000
Race/Ethnicity	30	2.033	0.001
New/Retained	5	3.040	0.010
SW/FW	10	3.663	0.000
Gender	5	5.861	0.000
Children	10	10.559	0.000

# Chi-square test results
Q14_Gender_xtab

Gender	Agreement						Total
	Enjoyed everything	Finding close fishing spot	Needed equipment	Techniques for casting/reeling	Children had trouble	Didn’t catch anything
Male	432 20.7 %	363 17.4 %	120 5.8 %	174 8.4 %	61 2.9 %	240 11.5 %	1390 66.7 %
Female	282 13.5 %	123 5.9 %	26 1.2 %	101 4.9 %	32 1.5 %	128 6.1 %	692 33.1 %
Total	714 34.3 %	486 23.3 %	146 7 %	275 13.2 %	93 4.5 %	368 17.7 %	2082 100 %

Q14_Gender_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Gender + Agreement, svy_desn_Q14, statistic = "adjWald")
## F = 5.8613, ndf = 5, ddf = 2086, p-value = 2.193e-05

Q14_Age_xtab

Age	Agreement						Total
	Enjoyed everything	Finding close fishing spot	Needed equipment	Techniques for casting/reeling	Children had trouble	Didn’t catch anything
18-24	85 4.1 %	141 6.8 %	43 2.1 %	62 3 %	6 0.3 %	102 4.9 %	439 21.2 %
25-34	147 7.1 %	129 6.2 %	37 1.8 %	98 4.7 %	32 1.5 %	80 3.8 %	523 25.1 %
35-44	174 8.3 %	101 4.8 %	37 1.8 %	57 2.7 %	37 1.8 %	90 4.3 %	496 23.7 %
45-54	151 7.2 %	68 3.3 %	18 0.9 %	33 1.6 %	11 0.5 %	52 2.5 %	333 16 %
55-64	111 5.3 %	33 1.6 %	9 0.4 %	20 1 %	4 0.2 %	27 1.3 %	204 9.8 %
65 or older	48 2.3 %	14 0.7 %	3 0.1 %	6 0.3 %	3 0.1 %	15 0.7 %	89 4.2 %
Total	716 34.4 %	486 23.3 %	147 7.1 %	276 13.2 %	93 4.5 %	366 17.6 %	2084 100 %

Q14_Age_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Age + Agreement, svy_desn_Q14, statistic = "adjWald")
## F = 4.9096, ndf = 25, ddf = 2066, p-value = 2.314e-14

Q14_Race_Ethn_xtab

Race_Ethn	Agreement						Total
	Enjoyed everything	Finding close fishing spot	Needed equipment	Techniques for casting/reeling	Children had trouble	Didn’t catch anything
White	446 22.8 %	233 11.9 %	73 3.7 %	125 6.4 %	58 3 %	195 10 %	1130 57.8 %
Asian or Pacific Islander	20 1 %	24 1.2 %	10 0.5 %	22 1.1 %	2 0.1 %	15 0.8 %	93 4.7 %
Black or African American	20 1 %	28 1.4 %	14 0.7 %	17 0.9 %	8 0.4 %	24 1.2 %	111 5.6 %
Hispanic	14 0.7 %	16 0.8 %	3 0.2 %	12 0.6 %	6 0.3 %	12 0.6 %	63 3.2 %
Native American or Alaskan Native	22 1.1 %	11 0.6 %	7 0.4 %	11 0.6 %	3 0.2 %	8 0.4 %	62 3.3 %
Other Hispanic	34 1.7 %	43 2.2 %	6 0.3 %	16 0.8 %	4 0.2 %	10 0.5 %	113 5.7 %
White Hispanic	111 5.7 %	92 4.7 %	32 1.6 %	69 3.5 %	7 0.4 %	74 3.8 %	385 19.7 %
Total	667 34.1 %	447 22.8 %	145 7.4 %	272 13.9 %	88 4.5 %	338 17.3 %	1957 100 %

Q14_Race_Ethn_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Race_Ethn + Agreement, svy_desn_Q14, statistic = "adjWald")
## F = 2.0326, ndf = 30, ddf = 2061, p-value = 0.0008004

Q14_SW_FW2_xtab

SW_FW2	Agreement						Total
	Enjoyed everything	Finding close fishing spot	Needed equipment	Techniques for casting/reeling	Children had trouble	Didn’t catch anything
FW	291 13 %	215 9.6 %	56 2.5 %	130 5.8 %	42 1.9 %	175 7.8 %	909 40.6 %
SW	228 10.2 %	73 3.3 %	24 1.1 %	54 2.4 %	17 0.8 %	64 2.9 %	460 20.7 %
SW.FW	263 11.8 %	218 9.7 %	72 3.2 %	114 5.1 %	43 1.9 %	159 7.1 %	869 38.8 %
Total	782 34.9 %	506 22.6 %	152 6.8 %	298 13.3 %	102 4.6 %	398 17.8 %	2238 100 %

Q14_SW_FW2_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~SW_FW2 + Agreement, svy_desn_Q14, statistic = "adjWald")
## F = 3.6635, ndf = 10, ddf = 2081, p-value = 7.225e-05

Q14_COVID_xtab

COVID	Agreement						Total
	Enjoyed everything	Finding close fishing spot	Needed equipment	Techniques for casting/reeling	Children had trouble	Didn’t catch anything
No, I was already planning to go fishing.	681 30.3 %	385 17.1 %	118 5.2 %	180 8 %	65 2.9 %	281 12.5 %	1710 76 %
Yes	103 4.6 %	127 5.6 %	34 1.5 %	120 5.3 %	37 1.6 %	118 5.2 %	539 23.8 %
Total	784 34.9 %	512 22.8 %	152 6.8 %	300 13.3 %	102 4.5 %	399 17.7 %	2249 100 %

Q14_COVID_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~COVID + Agreement, svy_desn_Q14, statistic = "adjWald")
## F = 13.978, ndf = 5, ddf = 2086, p-value = 1.827e-13

Q14_Children_xtab

Children	Agreement						Total
	Enjoyed everything	Finding close fishing spot	Needed equipment	Techniques for casting/reeling	Children had trouble	Didn’t catch anything
18 or over	268 15.9 %	96 5.7 %	31 1.8 %	45 2.7 %	17 1 %	92 5.5 %	549 32.6 %
No children	219 13 %	232 13.8 %	77 4.6 %	152 9 %	5 0.3 %	151 9 %	836 49.7 %
Under 17	80 4.8 %	79 4.7 %	19 1.1 %	31 1.8 %	32 1.9 %	57 3.4 %	298 17.7 %
Total	567 33.7 %	407 24.2 %	127 7.5 %	228 13.5 %	54 3.2 %	300 17.8 %	1683 100 %

Q14_Children_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Children + Agreement, svy_desn_Q14, statistic = "adjWald")
## F = 10.559, ndf = 10, ddf = 2081, p-value < 2.2e-16

Most anglers reported that they enjoyed everything about fishing, while some anglers reported they had a hard time finding a good fishing spot close to home or that they didn’t catch anything. Female anglers reported that they enjoyed everything, while male anglers were more likely to report difficulty finding a close fishing spot and needing equipment. Anglers over 45 years were most likely to report that they enjoyed everything, while younger anglers had trouble finding a close fishing spot, had uncertainty about casting/reeling, and reported that their children had trouble fishing. Asian or Pacific Islander anglers and White Hispanic anglers had issues with casting/reeling, Black or African American anglers and Native American or Alaskan Native anglers needed equipment, Hispanic anglers children had trouble fishing, Other Hispanic anglers had trouble finding close fishing spots, and White anglers enjoyed everything. Freshwater anglers reported that they didn’t catch anything, saltwater anglers enjoyed everything, and anglers that fished in both needed equipment. Anglers prompted by COVID-19 to fish needed techniques for casting/reeling and anglers that were already planning to fish enjoyed everything. Finally, anglers with young children reported that their children had trouble fishing, anglers with adult children enjoyed everything, and anglers with no children needed techniques for casting/reeling.

Q15. What is the most important factor in enjoying your fishing experience? Please just list one.

# Create wordclouds for COVID-anglers and non-COVID anglers

# Non-COVID anglers
#Create a vector containing only the text
text <- rev2_wts %>%
    filter(New.Old == "Old") 
# Create a corpus  
text <- text$Q15_MostImp
# Create a corpus  
docs <- Corpus(VectorSource(text))

# Clean the text
docs <- docs %>%
  tm_map(removeNumbers) %>%
  tm_map(removePunctuation) %>%
  tm_map(stripWhitespace)
docs <- tm_map(docs, content_transformer(tolower))
docs <- tm_map(docs, removeWords, stopwords("english"))
docs <- tm_map(docs, removeWords, c("getting", "get","spending","something","able","just","one","day","enjoy","enjoying","fishing","experience","catch","say","immediate","must","factor")) 

# Create Document term matrix
dtm <- TermDocumentMatrix(docs) 
matrix <- as.matrix(dtm) 
words <- sort(rowSums(matrix),decreasing=TRUE) 
df <- data.frame(word = names(words),freq=words)
head(df, 3) # fish, time, catching, family, outdoors

# Find associations
findAssocs(dtm, terms = "time", corlimit = 0.1)
# family        spent      husband       others      quality 
# 0.35         0.21         0.18         0.15         0.15 

findAssocs(dtm, terms = "fish", corlimit = 0.1)
# top 5 words associated with fish
# catching  often      wins      size important
# 0.71      0.21      0.21      0.17      0.14
      

findAssocs(dtm, terms = "catching", corlimit = 0.1)    
# top 5 words associated with catching
# fish  whether     size      air     keep  keepers  species 
# 0.71     0.12     0.12     0.12     0.11     0.11     0.11


# COVID anglers
#Create a vector containing only the text
text_COVID <- rev2_wts %>%
  filter(New.Old == "New") 
text_COVID <- text_COVID$Q15_MostImp
# Create a corpus  
docs_COVID <- Corpus(VectorSource(text_COVID))

# Clean the text
docs_COVID <- docs_COVID %>%
  tm_map(removeNumbers) %>%
  tm_map(removePunctuation) %>%
  tm_map(stripWhitespace)
docs_COVID <- tm_map(docs_COVID, content_transformer(tolower))
docs_COVID <- tm_map(docs_COVID, removeWords, stopwords("english"))
docs_COVID <- tm_map(docs_COVID, removeWords, c("getting", "get","spending","something","able","just","one","day","enjoy","enjoying","fishing","experience","catch","say","immediate","must","factor","youre")) 

# Create Document term matrix
dtm_COVID <- TermDocumentMatrix(docs_COVID) 
matrix_COVID <- as.matrix(dtm_COVID) 
words_COVID <- sort(rowSums(matrix_COVID),decreasing=TRUE) 
df_COVID <- data.frame(word = names(words_COVID),freq=words_COVID)
head(df_COVID, 5) # fish, time, family, catching, relaxing

# Find associations
findAssocs(dtm_COVID, terms = "time", corlimit = 0.1)
# top 5 words associated with time
#family    quality      spend        son      loved       ones 
#0.38       0.23       0.19       0.17       0.15       0.15 

findAssocs(dtm_COVID, terms = "fish", corlimit = 0.1)
# words associated with fish
# catching     big         size        areas    efficient 
# 0.65         0.16         0.16         0.14         0.14 

findAssocs(dtm_COVID, terms = "family", corlimit = 0.1)    
# words associated with family
# time    friends   memories especially     brings
# 0.38       0.29       0.17       0.16       0.15       

# Compare top 20 words from both groups
df_comparison <- data.frame(c(COVID = df_COVID[1:20,],Not_COVID = df[1:20,]))
# Top 20 words are all the same, but ordered slightly differently

### Create one wordcloud for all anglers combined
#Create a vector containing only the text
text_all <- rev2_wts$Q15_MostImp
docs_all <- Corpus(VectorSource(text_all))

# Clean the text
docs_all <- docs_all %>%
  tm_map(removeNumbers) %>%
  tm_map(removePunctuation) %>%
  tm_map(stripWhitespace)
docs_all <- tm_map(docs_all, content_transformer(tolower))
docs_all <- tm_map(docs_all, removeWords, stopwords("english"))
docs_all <- tm_map(docs_all, removeWords, c("getting", "get","spending","something","able","just","one","day","enjoy","enjoying","fishing","experience","catch","say","immediate","must","factor")) 

# Create Document term matrix
dtm_all <- TermDocumentMatrix(docs_all) 
matrix_all <- as.matrix(dtm_all) 
words_all <- sort(rowSums(matrix_all),decreasing=TRUE) 
df_all <- data.frame(word = names(words_all),freq=words_all)
head(df_all, 3) # fish, time, catching, family, outdoors

# Comparison of top 20 words
tab_df(df_comparison)

COVID.word	COVID.freq	Not_COVID.word	Not_COVID.freq
fish	188	fish	92
time	183	time	81
family	183	catching	76
catching	135	family	69
relaxing	89	outdoors	48
outdoors	82	water	44
nature	79	relaxing	35
fun	59	friends	32
weather	53	nature	28
friends	51	weather	22
water	51	relax	21
good	44	good	20
outside	43	fun	20
kids	39	away	19
quiet	30	peace	16
relaxation	30	husband	16
peace	28	quiet	14
relax	25	kids	14
husband	21	outside	13
away	18	relaxation	12

# Generate wordcloud for all anglers
set.seed(122) # for reproducibility 
wordcloud(words = df_all$word, freq = df_all$freq, min.freq = 10,max.words=200, random.order=FALSE, rot.per=0.35,colors=brewer.pal(8, "Dark2"))

We compared the most frequent words used by new and retained anglers to describe the most important factors in enjoying their fishing experiences. The top 5 words used by new anglers (in order) were fish, time, family, catching, and relaxing. The top 5 words used by retained anglers (in order) were fish, time, catching, family, and outdoors. The top 20 words used by both new and retained anglers to describe their fishing experiences were exactly the same, although ordered slightly differently.

We created one word cloud for all anglers combined because there were not large differences between the groups.

Q16. The last time you went fishing, how long did it take you to get there from home? Choose only one answer.

# Order the factor
rev2_wts$Income <- relevel(rev2_wts$Income, ref = "Between $100,000 and $149,999")
rev2_wts$Q16_TravelTime <- as.factor(rev2_wts$Q16_TravelTime)
rev2_wts$Q16_TravelTime  <- ordered(rev2_wts$Q16_TravelTime, levels = c("Less than 30 minutes","Between 30 and 60 minutes", "Between 1 and 2 hours", "More than 2 hours"))
svy_desn_Q16 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts,na.rm = TRUE) # New survey design

# Ordinal regression
Q16_mod1 <- svyolr(Q16_TravelTime ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age + Coastal, design=svy_desn_Q16)
summary(Q16_mod1) # examine model results
VIF(Q16_mod1) # collinearity
Q16_mod2 <- svyolr(Q16_TravelTime ~ Coastal + Income + SW_FW2 + Race_Ethn, design=svy_desn_Q16)
summary(Q16_mod2) # examine model results
VIF(Q16_mod2) # collinearity

# Create summary table
Q16_coef_table <- coef(summary(Q16_mod2))
Q16_pval <- pnorm(abs(Q16_coef_table[, "t value"]),lower.tail = FALSE)* 2
Q16_AOR <- exp(coef(Q16_mod2))
Q16_summary_table <- cbind(Q16_coef_table, "p value" = round(Q16_pval,3),Odds_Ratio = Q16_AOR)

# Survey response proportions
Q16_prop <- rev2_wts %>%
  filter(!(is.na(Q16_TravelTime))) %>% # remove NA values from Gender & Age
  count(Q16_TravelTime) %>%
  mutate(freq = prop.table(n)) # get proportions

# Table of responses
Q16_prop %>% 
  kbl(caption = "Q16 How long did it take you to get to your fishing spot?",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q16 How long did it take you to get to your fishing spot?

Q16_TravelTime	n	freq
Less than 30 minutes	599	0.33
Between 30 and 60 minutes	521	0.29
Between 1 and 2 hours	395	0.22
More than 2 hours	283	0.16

# Summary table
Q16_summary_table %>% 
  kbl(caption = "Q16 Model Odds Ratios",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q16 Model Odds Ratios

	Value	Std. Error	t value	p value	Odds_Ratio
CoastalY	-0.52	0.16	-3.31	0.00	0.60
IncomeLess than $20,000	-1.02	0.37	-2.77	0.01	0.36
IncomeBetween $20,000 and $39,999	-0.86	0.27	-3.14	0.00	0.42
IncomeBetween $40,000 and $59,999	-0.07	0.21	-0.31	0.76	0.94
IncomeBetween $60,000 and $79,999	-0.38	0.20	-1.90	0.06	0.68
IncomeBetween $80,000 and $99,999	-0.13	0.22	-0.58	0.56	0.88
IncomeOver $150,000	0.09	0.20	0.47	0.64	1.10
SW_FW2SW	1.10	0.19	5.67	0.00	3.02
SW_FW2SW.FW	0.81	0.18	4.62	0.00	2.25
Race_EthnAsian or Pacific Islander	0.60	0.26	2.33	0.02	1.82
Race_EthnBlack or African American	0.08	0.20	0.40	0.69	1.09
Race_EthnHispanic	0.89	0.35	2.51	0.01	2.43
Race_EthnNative American or Alaskan Native	-0.48	0.59	-0.82	0.42	0.62
Race_EthnOther Hispanic	0.11	0.27	0.42	0.68	1.12
Race_EthnWhite Hispanic	-0.13	0.19	-0.68	0.50	0.88
Less than 30 minutes|Between 30 and 60 minutes	-0.52	0.17	-3.04	0.00	0.60
Between 30 and 60 minutes|Between 1 and 2 hours	0.76	0.17	4.40	0.00	2.13
Between 1 and 2 hours|More than 2 hours	1.95	0.18	10.77	0.00	7.00

# Histogram of responses
model.frame(Q16_mod2) %>% ggplot(aes(x = Q16_TravelTime, fill = Income)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q16. How long did it take you to get to your fishing spot?")

model.frame(Q16_mod2) %>% ggplot(aes(x = Q16_TravelTime, fill = SW_FW2)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q16. How long did it take you to get to your fishing spot?")

model.frame(Q16_mod2) %>% ggplot(aes(x = Q16_TravelTime, fill = Race_Ethn)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q16. How long did it take you to get to your fishing spot?")

The majority of anglers (62%) reported that they traveled less than 1 hour the last time they went fishing (less than 30 mins or between 30-60 mins), while 38% of anglers traveled more than 1 hour (between 1-2 hrs or 2+ hrs).

The odds of traveling further to fish were: 2.772 to 1.679 times less likely for anglers with an income <$20K-$39K, 1.098 times more likely for saltwater anglers, 3.016 times more likely for anglers that fish in both freshwater and saltwater, 2.249 times more likely for Asian or Pacific Islander anglers, 1.086 times more likely for Hispanic anglers, and 1.141 times less likely for anglers in coastal counties.

Q17. How far are you willing to travel to go fishing?

# Order the factor
rev2_wts$Q17_MaxTravel <- as.factor(rev2_wts$Q17_MaxTravel)
rev2_wts$Q17_MaxTravel  <- ordered(rev2_wts$Q17_MaxTravel, levels = c("Less than 30 minutes","Between 30 and 60 minutes", "Between 1 and 2 hours", "Between 2 and 4 hours","More than 4 hours"))
svy_desn_Q17 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts,na.rm = TRUE) # New survey design

# Ordinal regression
Q17_mod1 <- svyolr(Q17_MaxTravel ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age + Coastal, design=svy_desn_Q17)
summary(Q17_mod1) # examine model results
VIF(Q17_mod1) # collinearity
Q17_mod2 <- svyolr(Q17_MaxTravel ~ Income + COVID + SW_FW2 + Children + New.Old + Race_Ethn + Coastal, design=svy_desn_Q17)
summary(Q17_mod2) # examine model results
VIF(Q17_mod2) # collinearity
Q17_mod3 <- svyolr(Q17_MaxTravel ~ Income + COVID + SW_FW2 + New.Old + Coastal, design=svy_desn_Q17)
summary(Q17_mod3) # examine model results
VIF(Q17_mod3) # collinearity

# Create summary table
Q17_coef_table <- coef(summary(Q17_mod3))
Q17_pval <- pnorm(abs(Q17_coef_table[, "t value"]),lower.tail = FALSE)* 2
Q17_AOR <- exp(coef(Q17_mod3))
Q17_summary_table <- cbind(Q17_coef_table, "p value" = round(Q17_pval,3),Odds_ratio = Q17_AOR)

# Survey response proportions
Q17_prop <- rev2_wts %>%
  filter(!(is.na(Q17_MaxTravel))) %>% # remove NA values from Gender & Age
  count(Q17_MaxTravel) %>%
  mutate(freq = prop.table(n)) # get proportions

# Table of responses
Q17_prop %>% 
  kbl(caption = "Q17 How far are you willing to travel to go fishing?",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria")

Q17 How far are you willing to travel to go fishing?

Q17_MaxTravel	n	freq
Less than 30 minutes	68	0.04
Between 30 and 60 minutes	341	0.19
Between 1 and 2 hours	556	0.31
Between 2 and 4 hours	430	0.24
More than 4 hours	400	0.22

# Summary Table
Q17_summary_table %>% 
  kbl(caption = "Q17 Model Odds Ratios",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q17 Model Odds Ratios

	Value	Std. Error	t value	p value	Odds_ratio
IncomeLess than $20,000	-0.98	0.50	-1.98	0.05	0.38
IncomeBetween $20,000 and $39,999	-0.12	0.25	-0.49	0.62	0.88
IncomeBetween $40,000 and $59,999	0.07	0.21	0.32	0.75	1.07
IncomeBetween $60,000 and $79,999	0.06	0.22	0.29	0.77	1.07
IncomeBetween $80,000 and $99,999	-0.20	0.17	-1.14	0.26	0.82
IncomeOver $150,000	0.33	0.17	1.92	0.06	1.39
COVIDYes	-0.34	0.15	-2.19	0.03	0.72
SW_FW2SW	0.38	0.19	2.00	0.05	1.47
SW_FW2SW.FW	0.83	0.16	5.07	0.00	2.30
New.OldNew	-0.33	0.13	-2.45	0.01	0.72
CoastalY	-0.61	0.15	-4.02	0.00	0.55
Less than 30 minutes|Between 30 and 60 minutes	-3.79	0.26	-14.38	0.00	0.02
Between 30 and 60 minutes|Between 1 and 2 hours	-1.46	0.19	-7.63	0.00	0.23
Between 1 and 2 hours|Between 2 and 4 hours	-0.09	0.18	-0.52	0.61	0.91
Between 2 and 4 hours|More than 4 hours	1.03	0.19	5.55	0.00	2.79

# Histogram of responses
model.frame(Q17_mod3) %>% ggplot(aes(x = Q17_MaxTravel, fill = Income)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q17. How far are you willing to travel to go fishing?")

model.frame(Q17_mod3) %>% ggplot(aes(x = Q17_MaxTravel, fill = COVID)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q17. How far are you willing to travel to go fishing?")

model.frame(Q17_mod3) %>% ggplot(aes(x = Q17_MaxTravel, fill = New.Old)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q17. How far are you willing to travel to go fishing?")

model.frame(Q17_mod3) %>% ggplot(aes(x = Q17_MaxTravel, fill = SW_FW2)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q17. How far are you willing to travel to go fishing?")

model.frame(Q17_mod3) %>% ggplot(aes(x = Q17_MaxTravel, fill = Coastal)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q17. How far are you willing to travel to go fishing?")

Most anglers (77%) reported that they were willing to travel over an hour to go fishing with a travel time mode of between 1-2 hrs. Only 4% of anglers reported that they were only willing to drive up to 30 minutes to go fishing.

The odds of being willing to travel further to fish were: 2.666 times less likely for anglers with an income <$20K, 1.398 times less likely for anglers that were prompted by COVID-19 to start fishing, 1.468 times more likely for saltwater anglers, 2.298 times more likely for anglers that fish in both freshwater and saltwater, 1.391 times less likely for new anglers, 1.834 times less likely for anglers in coastal counties.

Q18 & 19. Please tell us how these features would affect your willingness to fish. For each row, indicate your level of agreement with the following statement: I am more likely to go fishing at a place where this feature, facility, or amenity is available:

# Clean Q18 and Q19 survey answers
rev2_wts <- rev2_wts %>%  
  mutate_at(vars(Q18A:Q18I), list(~recode_factor(.,"Strongly Disagree" = "-2","Disagree" = "-1","Neither Agree nor Disagree" = "0","Agree" = "1","Strongly Agree"= "2"))) %>%
  mutate(across(Q18A:Q18I, ~as.numeric(as.character(.)))) %>%
  mutate_at(vars(Q18A:Q18I), list(~.*weights)) %>%
  mutate(Motorboat_access = Q18A,
         Peaceful_areas = Q18B,
         Portable_restrooms = Q18C,
         Paddling_trail = Q18D,
         Fish_stockings = Q18E,
         Parking_area = Q18F,
         ADA_accessible = Q18G,
         Baby_changing_stations = Q18H,
         Shoreline_access = Q18I)

# Clean Q19 survey answers
rev2_wts <- rev2_wts %>%  
  mutate_at(vars(Q19A:Q19I), list(~recode_factor(.,"Strongly Disagree" = "-2","Disagree" = "-1","Neither Agree nor Disagree" = "0","Agree" = "1","Strongly Agree"= "2"))) %>%
  mutate(across(Q19A:Q19I, ~as.numeric(as.character(.)))) %>%
  mutate_at(vars(Q19A:Q19I), list(~.*weights)) %>%
  mutate(Canoe_access = Q19A,
         Space = Q19B,
         Brick_restrooms = Q19C,
         Lighting = Q19D,
         Picnic_areas = Q19E,
         Playground = Q19F,
         Shade = Q19G,
         Swimming = Q19H,
         Camping = Q19I)

# Remove rows with missing independent variable data
Q18_19_df1 <- rev2_wts %>%
  filter_at(vars(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2,Coastal),all_vars(!is.na(.))) 
# Select independent variables
Q18_19_indep1 <- Q18_19_df1 %>%
  dplyr::select(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2,Coastal)
# Select dependent variables
Q18_19_dep1 <- Q18_19_df1 %>%
  dplyr::select(Motorboat_access,Peaceful_areas,Portable_restrooms,Paddling_trail,Fish_stockings,Parking_area,ADA_accessible,Baby_changing_stations,Shoreline_access,Canoe_access,Space,Brick_restrooms,Lighting,Picnic_areas,Playground,Shade,Swimming,Camping)
# Transform dependent data
Q18_19_dep2 <- decostand(Q18_19_dep1, method = 'normalize') # normalize
# MCA, keep 2 axes, on independent data
Q18_19_indep2 <- dudi.mix(Q18_19_indep1, scannf=F, nf=2)$tab # MCA, keep 2 axes

# Run RDA
Q18_19_rda_mod1 <- rda(Q18_19_dep2 ~ ., data = Q18_19_indep2,scale = TRUE) 
# remove variables with species scores for RDA1-2 <0.1
summary(Q18_19_rda_mod1)
Q18_19_dep3 <- Q18_19_df1 %>%
  dplyr::select(Motorboat_access,Fish_stockings,Baby_changing_stations,Canoe_access,Picnic_areas,Playground,Shade)
Q18_19_dep4 <- decostand(Q18_19_dep3, method = 'normalize') # normalize
Q18_19_rda_mod1 <- rda(Q18_19_dep4 ~ ., data = Q18_19_indep2,scale = TRUE) #redo RDA w/fewer dep vars
envfit(Q18_19_rda_mod1,Q18_19_indep2) 
Q18_19_indep3 <- Q18_19_indep2 %>% # select variables at alpha = 0.01 for next model 
  dplyr::select(Incom.Between..40.000.and..59.999,Incom.Over..150.000,Age.45.54,Age.25.34,Age.35.44,Age.55.64,Age.65.or.older ,Child.Under.17,Child.18.or.over,Child.No.children,New.O.Old,New.O.New,Race_.White,Race_.White.Hispanic,SW_FW.FW,SW_FW.SW)
Q18_19_rda_mod2 <- rda(Q18_19_dep4 ~ ., data = Q18_19_indep3)
envfit(Q18_19_rda_mod2,Q18_19_indep3) # select variables at alpha = 0.001 for next model
Q18_19_indep4 <- Q18_19_indep2 %>%
  dplyr::select(Incom.Over..150.000,Age.25.34,Age.35.44,Age.55.64,Child.Under.17,New.O.New,New.O.Old,Race_.White,SW_FW.FW,SW_FW.SW)
Q18_19_rda_mod3 <- rda(Q18_19_dep4 ~ ., data = Q18_19_indep4)
envfit(Q18_19_rda_mod3,Q18_19_indep4) # select variables at alpha = 0.001 for next model
Q18_19_rda_mod3$call
RsquareAdj(Q18_19_rda_mod3) # 0.06 
sqrt(vif.cca(Q18_19_rda_mod3))

# Extract scores for visualization
# Basic plot for visualization
vare_spp_sco <- scores(Q18_19_rda_mod3, display = "species",scaling = "species")
vare_sam_sco <- scores(Q18_19_rda_mod3, display = "sites",scaling = "sites")
vare_env_sco <- scores(Q18_19_rda_mod3, display = "bp",scaling = "sites")
vare_spp_tbl <- as_tibble(vare_spp_sco)
vare_sam_tbl <- as_tibble(vare_sam_sco)
vare_env_tbl <- as_tibble(vare_env_sco)
vare_spp_tbl <- mutate(vare_spp_tbl, vgntxt=rownames(vare_spp_sco), ccatype = "species")
vare_sam_tbl <- mutate(vare_sam_tbl, vgntxt=rownames(vare_sam_sco), ccatype = "sites")
vare_env_tbl <- mutate(vare_env_tbl, vgntxt=rownames(vare_env_sco), ccatype = "bp")

vare_env_tbl <- mutate(vare_env_tbl, vgntxt=c("Income over $150K","Age 25-34","Age 35-44","Age 55-64","Children Under 17","New Angler","White","Freshwater","Saltwater"),ccatype = "bp")
rescaled <- vare_env_tbl %>% 
            dplyr::select(RDA1, RDA2) %>%
            as.matrix() *5
vare_tbl <- dplyr::select(vare_env_tbl, vgntxt, ccatype) %>%
            bind_cols(as_tibble(rescaled)) %>%
            bind_rows(vare_spp_tbl)

# Get mean, SD, and calculate cronbach's alpha
# Recode factors 1-5 for comparison to other questions in final table
Q18_19_summary1 <- rev1_wts %>%  
  mutate_at(vars(Q18A:Q19I), list(~recode_factor(.,"Strongly Disagree" = "1","Disagree" = "2","Neither Agree nor Disagree" = "3","Agree" = "4","Strongly Agree"= "5"))) %>%
  mutate(across(Q18A:Q19I, ~as.numeric(as.character(.)))) %>%
  mutate_at(vars(Q18A:Q19I), list(~.*weights)) %>%
  mutate(Motorboat_access = Q18A,
         Peaceful_areas = Q18B,
         Portable_restrooms = Q18C,
         Paddling_trail = Q18D,
         Fish_stockings = Q18E,
         Parking_area = Q18F,
         ADA_accessible = Q18G,
         Baby_changing_stations = Q18H,
         Shoreline_access = Q18I,
         Canoe_access = Q19A,
         Space = Q19B,
         Brick_restrooms = Q19C,
         Lighting = Q19D,
         Picnic_areas = Q19E,
         Playground = Q19F,
         Shade = Q19G,
         Swimming = Q19H,
         Camping = Q19I) %>%
  dplyr::select(Motorboat_access,Peaceful_areas,Portable_restrooms,Paddling_trail,Fish_stockings,Parking_area,ADA_accessible,Baby_changing_stations,Shoreline_access,Canoe_access,Space,Brick_restrooms,Lighting,Picnic_areas,Playground,Shade,Swimming,Camping)

Q18_19_summary2 <- describe(Q18_19_summary1,na.rm = TRUE,IQR = TRUE) 
Q18_19_summary2$Variables <- row.names(Q18_19_summary2) 
Q18_19_summary3 <- Q18_19_summary2[,c("Variables","median","IQR")]
Q18_19_summary3$Variables<-str_replace_all(Q18_19_summary3$Variables, c("_" = " ")) # replaces spaces in names with periods

# Calculate cronbach's alpha
Q18_19_cronbach.alpha <- cronbach.alpha(abs(Q18_19_summary1),na.rm = TRUE,standardized = TRUE)
Q18_19_cronbach.alpha$alpha

# Pull out significant variable terms for each axis
Q18_19_envfit1 <- envfit(Q18_19_rda_mod3,Q18_19_indep4) # select variables at alpha = 0.001 for next model
Q18_19_envfit2 <- data.frame(Q18_19_envfit1$vectors[1],Q18_19_envfit1$vectors[2],Q18_19_envfit1$vectors[4])
names(Q18_19_envfit2) <- c("RDA 1","RDA 2","r2","p")
Q18_19_envfit2$Variables <- row.names(Q18_19_envfit2) 
Q18_19_envfit2 <- Q18_19_envfit2 %>% relocate(Variables, .before = "RDA 1")

# Descriptive stats table
tab_df(Q18_19_summary3)

Variables	median	IQR
Motorboat access	2.58	2.12
Peaceful areas	3.36	2.54
Portable restrooms	3.09	2.17
Paddling trail	2.44	2.28
Fish stockings	3.15	2.25
Parking area	3.25	2.24
ADA accessible	2.26	1.72
Baby changing stations	1.91	1.88
Shoreline access	3.25	2.49
Canoe access	2.54	2.29
Space	3.34	2.49
Brick restrooms	2.69	2.29
Lighting	3.09	2.28
Picnic areas	2.69	2.16
Playground	2.14	1.99
Shade	3.15	2.25
Swimming	2.44	2.39
Camping	2.67	2.16

# Reorganize data to plot
Q18_a <- rev1_wts[c(which(colnames(rev1_wts)=="Q18A"):which(colnames(rev1_wts)=="Q18I"))]
colnames(Q18_a) <- c("Motorboat_access","Peaceful_areas","Portable_restrooms","Paddling_trail","Fish_stockings","Parking_area","ADA_accessible","Baby_changing_stations","Shoreline_access") 
Q18_b <- Q18_a %>% 
  pivot_longer(cols = Motorboat_access:Shoreline_access,names_to = "Q18_answers",values_to = "Agreement")
Q18_c <- Q18_b %>%
  group_by(Q18_answers,Agreement) %>%
  count() %>%
  na.omit()
Q18_c$Agreement <- ordered(Q18_c$Agreement,levels = c("Strongly Agree", "Agree", "Neither Agree nor Disagree","Disagree","Strongly Disagree") )

ggplot(Q18_c, aes(fill=Agreement, y=n, x=Q18_answers)) + geom_bar(position="fill", stat="identity") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8),axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + labs(y = "Percent", x = "Q18. Agreement on Fishing Statements") 

# Reorganize data to plot
Q19_a <- rev1_wts[c(which(colnames(rev1_wts)=="Q19A"):which(colnames(rev1_wts)=="Q19H"))]
colnames(Q19_a) <- c("Canoe_access","Space","Brick_restrooms","Lighting","Picnic_areas","Playground","Shade","Swimming") 
Q19_b <- Q19_a %>% 
  pivot_longer(cols = Canoe_access:Swimming,names_to = "Q19_answers",values_to = "Agreement")
Q19_c <- Q19_b %>%
  group_by(Q19_answers,Agreement) %>%
  count() %>%
  na.omit()
Q19_c$Agreement <- ordered(Q19_c$Agreement,levels = c("Strongly Agree", "Agree", "Neither Agree nor Disagree","Disagree","Strongly Disagree") )

ggplot(Q19_c, aes(fill=Agreement, y=n, x=Q19_answers)) + geom_bar(position="fill", stat="identity") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8),axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + labs(y = "Percent", x = "Q19. Agreement on Fishing Statements")

# PERMANOVA
anova.cca(Q18_19_rda_mod3, step = 1000, by = "axis",permutations = 999)

## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = Q18_19_dep4 ~ Incom.Over..150.000 + Age.25.34 + Age.35.44 + Age.55.64 + Child.Under.17 + New.O.New + New.O.Old + Race_.White + SW_FW.FW + SW_FW.SW, data = Q18_19_indep4)
##           Df Variance       F Pr(>F)    
## RDA1       1  0.02585 38.0469  0.001 ***
## RDA2       1  0.01068 15.7149  0.001 ***
## RDA3       1  0.00427  6.2805  0.020 *  
## RDA4       1  0.00157  2.3047  0.720    
## RDA5       1  0.00147  2.1696  0.656    
## RDA6       1  0.00034  0.5078  1.000    
## RDA7       1  0.00003  0.0444  1.000    
## Residual 916  0.62224                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_df(Q18_19_envfit2)

Variables	RDA.1	RDA.2	r2	p
Incom.Over..150.000	0.98	-0.22	0.03	0.00
Age.25.34	-1.00	0.03	0.03	0.00
Age.35.44	-0.88	-0.48	0.02	0.00
Age.55.64	0.93	-0.36	0.02	0.00
Child.Under.17	-0.77	-0.64	0.09	0.00
New.O.New	-0.79	0.61	0.04	0.00
New.O.Old	0.79	-0.61	0.04	0.00
Race_.White	1.00	0.00	0.04	0.00
SW_FW.FW	-0.69	0.73	0.02	0.00
SW_FW.SW	0.52	-0.86	0.04	0.00

# Basic plot for visualization
ggplot() +
  geom_point(aes(x=RDA1, y=RDA2), data=filter(vare_tbl, ccatype=="species"))  +
  geom_text_repel(aes(x=RDA1, y=RDA2, label=vgntxt, size=3,colour=ccatype),
                  data=vare_tbl, seed=123) + 
  geom_segment(aes(x=0, y=0, xend=RDA1, yend=RDA2), arrow=arrow(length = unit(0.2,"cm")),
               data=filter(vare_tbl, ccatype=="bp"), color="blue") +
  coord_fixed() +
  theme_classic() +
  scale_colour_manual(values = c("blue", "black")) +
  theme(legend.position="none")

Anglers agreed more strongly that they would fish in areas that had the following amenities: peaceful areas,space, parking areas, and shoreline access. Baby changing stations, playgrounds, ADA accessible facilities were least important overall.

Anglers between the age of 25-34, new anglers, and anglers that primarily fish in freshwater more strongly favored shade, canoe access, fish stockings, and picnic areas than typical anglers. Anglers aged 55-64 yrs, White anglers, saltwater anglers, and anglers with incomes over $150K were more likely to favor sites with motorboat access. Finally, anglers with children under 17 and anglers between the ages of 35-44 yrs were more likely to favor sites with baby changing stations and playgrounds.

Cronbach’s alpha for Q18-19 is 0.9937515

Q20. How likely are you to purchase another fishing license when your current license expires?

# Order the factor
rev2_wts$Q20_NewLicense <- recode_factor(rev2_wts$Q20_NewLicense,"Not likely at all" = "Not likely","Not very likely" = "Not likely")
rev2_wts$Q20_NewLicense  <- ordered(rev2_wts$Q20_NewLicense, levels = c("Not likely", "Somewhat likely", "Very likely"))
svy_desn_Q20 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts,na.rm = TRUE) # New survey design

# Ordinal regression
Q20_mod1 <- svyolr(Q20_NewLicense ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age + Coastal, design=svy_desn_Q20)
summary(Q20_mod1) # examine model results
VIF(Q20_mod1) # collinearity
Q20_mod2 <- svyolr(Q20_NewLicense ~ New.Old + COVID + SW_FW2, design=svy_desn_Q20)
summary(Q20_mod2) # examine model results
VIF(Q20_mod2) # collinearity

# Survey response proportions
Q20_prop <- rev2_wts %>%
  filter(!(is.na(Q20_NewLicense))) %>% # remove NA values from Gender & Age
  count(Q20_NewLicense) %>%
  mutate(freq = prop.table(n)) # get proportions

# Create summary table
Q20_coef_table <- coef(summary(Q20_mod2))
Q20_pval <- pnorm(abs(Q20_coef_table[, "t value"]),lower.tail = FALSE)* 2
Q20_AOR <- exp(coef(Q20_mod2))
Q20_summary_table <- cbind(Q20_coef_table, "p value" = round(Q20_pval,3),Odds_Ratio = Q20_AOR)

# Table of responses
Q20_prop %>% 
  kbl(caption = "Q20 How likely are you to purchase another fishing license?",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q20 How likely are you to purchase another fishing license?

Q20_NewLicense	n	freq
Not likely	18	0.01
Somewhat likely	139	0.07
Very likely	1798	0.92

# Summary table
Q20_summary_table %>% 
  kbl(caption = "Q20 Model Odds Ratios",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q20 Model Odds Ratios

	Value	Std. Error	t value	p value	Odds_Ratio
New.OldNew	-1.21	0.31	-3.93	0.00	0.30
COVIDYes	-1.04	0.25	-4.17	0.00	0.35
SW_FW2SW	-0.40	0.28	-1.43	0.15	0.67
SW_FW2SW.FW	0.63	0.29	2.16	0.03	1.88
Not likely|Somewhat likely	-6.18	0.51	-12.18	0.00	0.00
Somewhat likely|Very likely	-3.81	0.34	-11.09	0.00	0.02

# Histogram of responses
model.frame(Q20_mod2) %>% ggplot(aes(x = Q20_NewLicense, fill = New.Old)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q20. How likely are you to purchase another license?")

model.frame(Q20_mod2) %>% ggplot(aes(x = Q20_NewLicense, fill = COVID)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q20. How likely are you to purchase another license?")

model.frame(Q20_mod2) %>% ggplot(aes(x = Q20_NewLicense, fill = SW_FW2)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q20. How likely are you to purchase another license?")

Most anglers (99%) reported that they were either somewhat likely (7%) or very likely (92%) to purchase another fishing license.

The odds of being more likely to purchase another fishing license were: 3.366 times less likely for new anglers that purchased a license during the pandemic, 2.841 times less likely for anglers that were prompted by COVID-19 to purchase a license, 1.884 times more likely for anglers that fish in both freshwater and saltwater.

Q21. How likely would you be to enroll in a program where your license renews automatically every year?

# Order the factor
rev2_wts$Q21_AutoRenew <- as.factor(rev2_wts$Q21_AutoRenew)
rev2_wts$Q21_AutoRenew  <- ordered(rev2_wts$Q21_AutoRenew, levels = c("Extremely unlikely","Unlikely", "Not sure", "Likely", "Extremely likely"))
svy_desn_Q21 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts,na.rm = TRUE) # New survey design

# Ordinal regression
Q21_mod1 <- svyolr(Q21_AutoRenew ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age + Coastal, design=svy_desn_Q21)
summary(Q21_mod1) # examine model results
VIF(Q21_mod1) # collinearity
Q21_mod2 <- svyolr(Q21_AutoRenew ~ Income + Race_Ethn, design=svy_desn_Q21)
summary(Q21_mod2) # examine model results
VIF(Q21_mod2) # collinearity

# Survey response proportions
Q21_prop <- rev2_wts %>%
  filter(!(is.na(Q21_AutoRenew))) %>% # remove NA values from Gender & Age
  count(Q21_AutoRenew) %>%
  mutate(freq = prop.table(n)) # get proportions

# Create summary table
Q21_coef_table <- coef(summary(Q21_mod2))
Q21_pval <- pnorm(abs(Q21_coef_table[, "t value"]),lower.tail = FALSE)* 2
Q21_AOR <- exp(coef(Q21_mod2))
Q21_summary_table <- cbind(Q21_coef_table, "p value" = round(Q21_pval,3),Odds_ratio = Q21_AOR)

# Table of responses
Q21_prop %>% 
  kbl(caption = "Q21 Likeliness to enroll in auto-renew program",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q21 Likeliness to enroll in auto-renew program

Q21_AutoRenew	n	freq
Extremely unlikely	28	0.01
Unlikely	114	0.06
Not sure	471	0.24
Likely	546	0.28
Extremely likely	788	0.40

# Summary table
Q21_summary_table %>% 
  kbl(caption = "Q21 Model Odds Ratios",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q21 Model Odds Ratios

	Value	Std. Error	t value	p value	Odds_ratio
IncomeLess than $20,000	-0.68	0.39	-1.72	0.09	0.51
IncomeBetween $20,000 and $39,999	-0.76	0.27	-2.84	0.00	0.47
IncomeBetween $40,000 and $59,999	-0.39	0.20	-1.94	0.05	0.67
IncomeBetween $60,000 and $79,999	0.01	0.19	0.05	0.96	1.01
IncomeBetween $80,000 and $99,999	-0.07	0.21	-0.31	0.76	0.94
IncomeOver $150,000	0.45	0.18	2.49	0.01	1.57
Race_EthnAsian or Pacific Islander	0.06	0.28	0.21	0.83	1.06
Race_EthnBlack or African American	-0.67	0.32	-2.08	0.04	0.51
Race_EthnHispanic	0.05	0.36	0.15	0.88	1.06
Race_EthnNative American or Alaskan Native	-0.28	0.80	-0.35	0.72	0.75
Race_EthnOther Hispanic	0.91	0.33	2.74	0.01	2.49
Race_EthnWhite Hispanic	0.16	0.18	0.89	0.37	1.17
Extremely unlikely|Unlikely	-5.68	0.41	-13.85	0.00	0.00
Unlikely|Not sure	-2.99	0.19	-16.01	0.00	0.05
Not sure|Likely	-1.14	0.14	-8.37	0.00	0.32
Likely|Extremely likely	0.18	0.13	1.35	0.18	1.20

# Histogram of responses
model.frame(Q21_mod2) %>% ggplot(aes(x = Q21_AutoRenew, fill = Income)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q21. How likely are you to enroll in auto-renew license program?")

model.frame(Q21_mod2) %>% ggplot(aes(x = Q21_AutoRenew, fill = Race_Ethn)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q21. How likely are you to enroll in auto-renew license program?")

Most anglers (68%) reported that they were either likely (28%) or extremely likely (40%) to enroll in an auto-renew program, while only 7% were unlikely and 24% were unsure.

The odds of being likely to enroll in an autorenew program were: 1.482 to 2.131 times less likely for anglers with an income between $20K-$59K, 1.57 times more likely for anglers with an income over $150K, 1.964 times less likely for anglers that identified as Black or African American, and 2.487 times more likely for anglers that identified as Other Hispanic.

Q24. Which types of fishing information would you find useful? Please choose all that you would find useful.

# Evaluate frequencies by categorical responses
rev2_wts$Income <- factor(rev2_wts$Income,levels = c("Less than $20,000", "Between $20,000 and $39,999", "Between $40,000 and $59,999","Between $60,000 and $79,999","Between $80,000 and $99,999","Between $100,000 and $149,999","Over $150,000"))
Q24_df1 <- rev2_wts %>% 
  dplyr::select(c("Q24A":"Q24I")) %>%
  pivot_longer(cols = "Q24A":"Q24I",names_to = "Q24_answers",values_to = "Agreement") %>%
  group_by(Q24_answers,Agreement) %>%
  count() %>%
  na.omit()
freq = prop.table(Q24_df1$n)
Q24_df1 <- data.frame(cbind(Q24_df1,freq=freq))
Q24_df1$Agreement <- recode_factor(Q24_df1$Agreement,"I wouldn't be interested in receiving any information." = "Not interested","How to improve my fishing skills" = "Improve skills","What species to fish for" = "What species","What different kinds of bait or lures work best for certain fish" = "Baits and lures","How to clean fish" = "Clean fish","How to buy a fishing license" = "Buy license","What amenities are available at local fishing areas" = "Available amenities","How to understand fishing regulations better" = "Understand Regulations")

# Format data for analysis
Q24_df2 <- rev2_wts %>% 
    dplyr::select(Q24A,Q24B,Q24C,Q24D,Q24E,Q24F,Q24G,Q24H,Q24I,Gender,Age,Race_Ethn,Income,SW_FW2,New.Old,Children,COVID,Coastal,weights) %>%
  pivot_longer(cols = "Q24A":"Q24I",names_to = "Q24_answers",values_to = "Agreement") %>%
  group_by(Q24_answers,Agreement,Gender,Age,Race_Ethn,Income,SW_FW2,New.Old,Children,COVID,Coastal,weights) %>%
  filter(!is.na(Agreement)) %>%
  mutate(Agreement = as.factor(Agreement)) %>%
  data.frame()

Q24_df2$Agreement <- recode_factor(Q24_df2$Agreement,"I wouldn't be interested in receiving any information." = "Not interested","How to improve my fishing skills" = "Improve skills","What species to fish for" = "What species","What different kinds of bait or lures work best for certain fish" = "Baits and lures","How to clean fish" = "Clean fish","How to buy a fishing license" = "Buy license","What amenities are available at local fishing areas" = "Available amenities","How to understand fishing regulations better" = "Understand Regulations")

Q24_df2$Age <- ordered(Q24_df2$Age,levels = c("18-24", "25-34", "35-44","45-54","55-64","65 or older"))

Q24_df2$Income <- ordered(Q24_df2$Income,levels = c("Less than $20,000", "Between $20,000 and $39,999", "Between $40,000 and $59,999","Between $60,000 and $79,999","Between $80,000 and $99,999","Between $100,000 and $149,999","Over $150,000"))

# Chi-square tests
svy_desn_Q24 <- svydesign(ids = ~0, weights = ~weights, data = Q24_df2) # Survey design

Q24_Gender_chisq <- svychisq(~Gender + Agreement, svy_desn_Q24, statistic = "adjWald") # not significant
Q24_Age_chisq <- svychisq(~Age + Agreement, svy_desn_Q24, statistic = "adjWald") # significant
Q24_Race_Ethn_chisq <- svychisq(~Race_Ethn + Agreement, svy_desn_Q24, statistic = "adjWald") # not significant
Q24_Income_chisq <- svychisq(~Income + Agreement, svy_desn_Q24, statistic = "adjWald") # not significant
Q24_SW_FW2_chisq <- svychisq(~SW_FW2 + Agreement, svy_desn_Q24, statistic = "adjWald") # not significant
Q24_New.Old_chisq <- svychisq(~New.Old + Agreement, svy_desn_Q24, statistic = "adjWald") # not significant
Q24_Children_chisq <- svychisq(~Children + Agreement, svy_desn_Q24, statistic = "adjWald") # significant
Q24_COVID_chisq <- svychisq(~COVID + Agreement, svy_desn_Q24, statistic = "adjWald") # not significant
Q24_Coastal_chisq <- svychisq(~Coastal + Agreement, svy_desn_Q24, statistic = "adjWald") # not significant

# Contingency tables for significant factors
Q24_Age_xtab <- tab_xtab(Q24_df2$Age,Q24_df2$Agreement,weight.by = Q24_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)
Q24_Children_xtab <- tab_xtab(Q24_df2$Children,Q24_df2$Agreement,weight.by = Q24_df2$weights,show.cell.prc = TRUE,show.summary = FALSE)

# Correlation plots of residuals to examine results
Q24_corr_Age <- corrplot(Q24_Age_chisq$residuals, is.cor = FALSE,method = 'number')

Q24_corr_Children <- corrplot(Q24_Children_chisq$residuals, is.cor = FALSE,method = 'number')

# Combine chi-squared test results into dataframe
Q24_Age <- data.frame(cbind("Age",Q24_Age_chisq[[c(2,1)]],Q24_Age_chisq$statistic,Q24_Age_chisq$p.value))
Q24_Children <- data.frame(cbind("Children",Q24_Children_chisq[[c(2,1)]],Q24_Children_chisq$statistic,Q24_Children_chisq$p.value))

Q24_chi_sq <- data.frame(rbind(Q24_Age,Q24_Children))
names(Q24_chi_sq) <- c("Variable", "df","F","p")
Q24_chi_sq$F <- as.numeric(Q24_chi_sq$F)
Q24_chi_sq$p <- as.numeric(Q24_chi_sq$p)

# Response percentages graph
ggplot(Q24_df1,aes(x=Agreement,y=freq)) + geom_bar(stat="identity") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Q24. Which types of fishing information would you find useful?") 

tab_df(Q24_chi_sq,digits = 3,title = "Q24. Which types of fishing information would you find useful?")

Q24. Which types of fishing information would you find useful?

Variable	df	F	p
Age	40	2.092	0.000
Children	16	2.478	0.001

# Chi-square test results
Q24_Age_xtab

Age	Agreement									Total
	Not interested	Improve skills	What species	Baits and lures	Clean fish	Buy license	Available amenities	Understand Regulations	Where to fish
18-24	29 0.4 %	248 3.8 %	172 2.6 %	252 3.8 %	147 2.2 %	38 0.6 %	144 2.2 %	100 1.5 %	224 3.4 %	1354 20.5 %
25-34	24 0.4 %	259 3.9 %	209 3.2 %	294 4.5 %	183 2.8 %	30 0.5 %	186 2.8 %	130 2 %	270 4.1 %	1585 24.2 %
35-44	50 0.8 %	239 3.6 %	188 2.9 %	292 4.4 %	142 2.2 %	55 0.8 %	192 2.9 %	122 1.8 %	267 4 %	1547 23.4 %
45-54	44 0.7 %	158 2.4 %	125 1.9 %	216 3.3 %	91 1.4 %	28 0.4 %	143 2.2 %	97 1.5 %	199 3 %	1101 16.8 %
55-64	29 0.4 %	96 1.5 %	73 1.1 %	135 2 %	53 0.8 %	22 0.3 %	105 1.6 %	59 0.9 %	117 1.8 %	689 10.4 %
65 or older	19 0.3 %	42 0.6 %	32 0.5 %	70 1.1 %	19 0.3 %	8 0.1 %	45 0.7 %	32 0.5 %	53 0.8 %	320 4.9 %
Total	195 3 %	1042 15.8 %	799 12.1 %	1259 19.1 %	635 9.6 %	181 2.7 %	815 12.4 %	540 8.2 %	1130 17.1 %	6596 100 %

Q24_Age_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Age + Agreement, svy_desn_Q24, statistic = "adjWald")
## F = 2.0924, ndf = 40, ddf = 6369, p-value = 6.889e-05

Q24_Age_xtab

Age	Agreement									Total
	Not interested	Improve skills	What species	Baits and lures	Clean fish	Buy license	Available amenities	Understand Regulations	Where to fish
18-24	29 0.4 %	248 3.8 %	172 2.6 %	252 3.8 %	147 2.2 %	38 0.6 %	144 2.2 %	100 1.5 %	224 3.4 %	1354 20.5 %
25-34	24 0.4 %	259 3.9 %	209 3.2 %	294 4.5 %	183 2.8 %	30 0.5 %	186 2.8 %	130 2 %	270 4.1 %	1585 24.2 %
35-44	50 0.8 %	239 3.6 %	188 2.9 %	292 4.4 %	142 2.2 %	55 0.8 %	192 2.9 %	122 1.8 %	267 4 %	1547 23.4 %
45-54	44 0.7 %	158 2.4 %	125 1.9 %	216 3.3 %	91 1.4 %	28 0.4 %	143 2.2 %	97 1.5 %	199 3 %	1101 16.8 %
55-64	29 0.4 %	96 1.5 %	73 1.1 %	135 2 %	53 0.8 %	22 0.3 %	105 1.6 %	59 0.9 %	117 1.8 %	689 10.4 %
65 or older	19 0.3 %	42 0.6 %	32 0.5 %	70 1.1 %	19 0.3 %	8 0.1 %	45 0.7 %	32 0.5 %	53 0.8 %	320 4.9 %
Total	195 3 %	1042 15.8 %	799 12.1 %	1259 19.1 %	635 9.6 %	181 2.7 %	815 12.4 %	540 8.2 %	1130 17.1 %	6596 100 %

Q24_Children_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Children + Agreement, svy_desn_Q24, statistic = "adjWald")
## F = 2.4778, ndf = 16, ddf = 6393, p-value = 0.0009007

Anglers were most interested in information about different kinds of bait or lures, followed by information about where to fish, and how to improve their fishing skills. Anglers were least interested in information about how to buy a fishing license. Very few anglers reported that they wouldn’t be interested in receiving any information.

The youngest anglers (ages 18-24) were most interested in receiving information about how to improve their fishing skills, anglers 25-34 wanted information about how to clean fish, anglers 35-44 were interested in information about how to buy a license, anglers 45+ were largely not interested in receiving any information. Anglers with young children wanted information about how to clean fish, anglers with adult children were not interested in receiving information, and those with no children were interested in improving their skills.

Q25. Did the survey respondent purchase another license in 2021 and/or 2022?

# Combine purchase years into one factor
LCYR_purchases <- rev2_wts %>%
  filter(!is.na(LCYR_2020)) %>% # not all RespIDs could be matched to the license database, exclude missing ones
  mutate(LC = paste(LCYR_2021,LCYR_2022,sep = "_"),
         LC21_22 = as.factor(case_when(grepl("2021",LC) & grepl("2022", LC) ~ "Both_Yrs",
                                                     grepl("2021", LC) ~ "One_Yr",
                                                     grepl("2022", LC) ~ "One_Yr",
                                                     grepl("NA_NA", LC) ~ "None")))

LCYR_purchases$LC21_22 <- as.factor(LCYR_purchases$LC21_22)
LCYR_purchases$LC21_22  <- ordered(LCYR_purchases$LC21_22, levels = c("None","One_Yr", "Both_Yrs"))
LCYR_purchases$Q20_NewLicense <- factor(LCYR_purchases$Q20_NewLicense, ordered = FALSE)
LCYR_purchases$Q20_NewLicense <- factor(LCYR_purchases$Q20_NewLicense,levels = c("Somewhat likely", "Not likely","Very likely"))
svy_desn_Q25 <- svydesign(ids = ~0, weights = ~weights, data = LCYR_purchases,na.rm = TRUE) # New survey design

# Ordinal regression
Q25_mod1 <- svyolr(LC21_22 ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age + Coastal+ Q20_NewLicense, design=svy_desn_Q25)
summary(Q25_mod1) # examine model results
VIF(Q25_mod1) # collinearity
Q25_mod2 <- svyolr(LC21_22 ~ Gender + New.Old + COVID + Age + Q20_NewLicense, design=svy_desn_Q25)
summary(Q25_mod2) # examine model results
VIF(Q25_mod2) # collinearity

# Create summary table
Q25_coef_table <- coef(summary(Q25_mod2))
Q25_pval <- pnorm(abs(Q25_coef_table[, "t value"]),lower.tail = FALSE)* 2
Q25_AOR <- exp(coef(Q25_mod2))
Q25_summary_table <- cbind(Q25_coef_table, "p value" = round(Q25_pval,3),Odds_Ratio = Q25_AOR)

# Survey response proportions New/Old
Q25_prop_New_Old_2021 <- rev2_wts %>%
  filter(!is.na(LCYR_2020)) %>%
  group_by(New.Old) %>%
  count(LCYR_2021)

Q25_prop_New_Old_2022 <- rev2_wts %>%
  filter(!is.na(LCYR_2020)) %>%
  group_by(New.Old) %>%
  count(LCYR_2022)

# Survey response proportions COVID-prompted
Q25_prop_COVID_2021 <- rev2_wts %>%
  filter(!is.na(LCYR_2020)) %>%
  group_by(COVID) %>%
  count(LCYR_2021)

Q25_prop_COVID_2022 <- rev2_wts %>%
  filter(!is.na(LCYR_2020)) %>%
  group_by(COVID) %>%
  count(LCYR_2022)

# Survey response proportions intention to purchase a license
Q25_prop_Q20_2021 <- rev2_wts %>%
  filter(!is.na(LCYR_2020)) %>%
  group_by(Q20_NewLicense) %>%
  count(LCYR_2021)

Q25_prop_Q20_2022 <- rev2_wts %>%
  filter(!is.na(LCYR_2020)) %>%
  group_by(Q20_NewLicense) %>%
  count(LCYR_2022)

# Summary table
Q25_summary_table %>% 
  kbl(caption = "Q25 Model Odds Ratios",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Q25 Model Odds Ratios

	Value	Std. Error	t value	p value	Odds_Ratio
GenderFemale	-0.37	0.12	-3.04	0.00	0.69
New.OldNew	-0.58	0.13	-4.60	0.00	0.56
COVIDYes	-0.31	0.16	-1.96	0.05	0.73
Age18-24	0.01	0.22	0.05	0.96	1.01
Age25-34	0.09	0.17	0.54	0.59	1.09
Age35-44	-0.01	0.13	-0.09	0.93	0.99
Age55-64	-0.01	0.15	-0.05	0.96	0.99
Age65 or older	1.20	0.20	6.15	0.00	3.33
Q20_NewLicenseNot likely	0.54	0.81	0.67	0.50	1.72
Q20_NewLicenseVery likely	1.44	0.27	5.34	0.00	4.23
None|One_Yr	-0.69	0.30	-2.27	0.02	0.50
One_Yr|Both_Yrs	1.36	0.31	4.42	0.00	3.88

# Histogram of responses
model.frame(Q25_mod2) %>% ggplot(aes(x = LC21_22, fill = Gender)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Did Respondents Purchase Another Fishing License in 2021 or 2022?")

model.frame(Q25_mod2) %>% ggplot(aes(x = LC21_22, fill = New.Old)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Did Respondents Purchase Another Fishing License in 2021 or 2022?")

model.frame(Q25_mod2) %>% ggplot(aes(x = LC21_22, fill = COVID)) + geom_bar(position = "fill") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + labs(y = "Percentage", x = "Did Respondents Purchase Another Fishing License in 2021 or 2022?")

Approximately 58% of new anglers and 54% of COVID-prompted anglers purchased a license in 2021, compared to 69% of retained anglers 64% of anglers not COVID-prompted. In 2022, 51% of new anglers and COVID-prompted anglers purchased a license, while 69% of retained anglers and 61% of anglers not prompted by COVID purchased a license.

Anglers that said they were not likely to purchase another fishing license had a 61% churn rate overall, those that were somewhat likely had a 68% churn rate, and those that were very likely to purchase another fishing license had a 41% churn rate.

The odds of purchasing a license in one year (2021 or 2022) or both years (2021 and 2022) were: 1.451 times less likely for female anglers, 1.782 times less likely for new anglers, 1.363 times less likely for COVID-prompted anglers, 3.334 times more likely for anglers 65 years or older, and 4.23 times more likely for anglers that reported they were very likely to purchase another fishing license compared to anglers that reported they were only somewhat likely to purchase another license.

Tables for publication

tab_xtab(rev2_wts$Income,rev2_wts$New.Old,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE,file = "Income.doc")

tab_xtab(rev2_wts$Children,rev2_wts$New.Old,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE,file = "Children.doc")

tab_xtab(rev2_wts$Gender,rev2_wts$New.Old,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE,file = "Gender.doc")

tab_xtab(rev2_wts$Age,rev2_wts$New.Old,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE,file = "Age.doc")

tab_xtab(rev2_wts$Race_Ethn,rev2_wts$New.Old,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE,file = "Race_Ethnicity.doc")

tab_df(weights2,file = "Age and Gender Weights.doc") # Numeric weights and frequencies of age and gender

tab_dfs(list(Q5_summary3,Q11_summary3,Q18_19_summary3),file = "Likert_summary_tables.doc") # Summary responses to likert scale questions

tab_dfs(list(Q5_envfit2,Q11_envfit2,Q18_19_envfit2),file = "RDA_summary_tables.doc") # summary of significant RDA variables

tab_dfs(list(Q7_chi_sq_SW_FW2,Q10_chi_sq,Q14_chi_sq,Q24_chi_sq),file = "Chi_sq_results.doc") # summary of significant chi-squared test results

# Corrected variable names for binomial regression table
var_names <- c("Intercept","Age:25-34 yrs","Age:35-44 yrs","Age:45-54 yrs","Age:55-64 yrs","Age:65 yrs or older","Children:No Children","Children:Under 17","Gender:Male","New/Retained:New","Race/Ethnicity:Black or African American","Race/Ethnicity:Hispanic","Race/Ethnicity:Native American or Alaskan Native","Race/Ethnicity:Other Hispanic","Race/Ethnicity:White","Race/Ethnicity:White Hispanic","SW/FW:SW","SW/FW:Both","Income:$20K-$39K","Income:$40K-$59K","Income:$60K-$79K","Income:$80K-$99K","Income:$100K-$149K","Income:>$150K","COVID:Yes")

tab_model(Q4_mod3,Q0_mod2,Q3_mod3,Q9_mod2,Q8_mod2,wrap.labels = FALSE,show.r2 = FALSE,order.terms = c(1,25,10,9,7,8,17,18,2,3,4,5,6,11,12,13,14,15,16,19,20,21,22,23,24),pred.labels = var_names,dv.labels = c("Prompted by COVID-19 to Fish","New Angler (i.e. Not Retained)","Fished in Last 6 months","Fishing was the Main Trip Purpose","Fished from a Boat"),file = "Binomial_logistic_reg.doc") # binomial logistic regression results

# Corrected variable names for ordinal regression table
var_names2 <- c("SW/FW:SW","SW/FW:Both","Gender:Male","Age:25-34 yrs","Age:35-44 yrs","Age:45-54 yrs","Age:55-64 yrs","Age:65 yrs or older","Days Fished:1 day|2-3 days","Days Fished:2-3 days|4-10 days","Days Fished:4-10 days|>10days","COVID:Yes","Race/Ethnicity:Black or African American","Race/Ethnicity:Hispanic","Race/Ethnicity:Native American or Alaskan Native","Race/Ethnicity:Other Hispanic","Race/Ethnicity:White","Race/Ethnicity:White Hispanic","Catch/Keep Fish:Never|Rarely","Catch/Keep Fish:Rarely|Sometimes","Catch/Keep Fish:Sometimes|Most of the time","Catch/Keep Fish:Most of the time|Always","New/Retained:New", "Income:$20K-$39K","Income:$40K-$59K","Income:$60K-$79K","Income:$80K-$99K","Income:$100K-$149K","Income:>$150K","Travel Time:<30 mins|30-60 mins","Travel Time:30-60 mins|1-2 hrs","Travel Time:1-2 hrs|>2 hrs","Travel Time:2-4 hrs|>4 hrs","Children:No Children","Children:Under 17","Travel Time:1-2hrs|2-4 hrs","Purchase License:Not likely|Somewhat likely","Purchase License:Somewhat likely|Very likely","Auto-renew:Extremely unlikely|Unlikely","Auto-renew:Unlikely|Not sure","Auto-renew:Not sure|Likely","Auto-renew:Likely|Extremely likely")

tab_model(Q6_mod3,Q12_mod3,Q13_mod2,Q16_mod2,Q17_mod3,Q20_mod2,Q21_mod2,Q25_mod2,pred.labels = var_names2,order.terms = c(12,23,3,34,35,1,2,4,5,6,7,8,13,14,15,16,17,18,24,25,26,27,28,29,9,10,11,19,20,21,22,30,31,32,33,36,37,38,39,40,41,42), show.r2 = FALSE,dv.labels = c("How many days did you fish in the last 6 months?","How often do you catch fish?","How often do you keep fish?","How far did you travel to fish?","How far would you travel to fish?","How likely are you to purchase another license?","How likely would you be to enroll in an auto-renew license program?","Did they purchase a license in 2021 or 2022?"),file = "Ordinal_logistic_reg.doc") # ordinal logistic regression results

# Save supplemental contingency tables 
# Q7
tab_xtab(rev2_wts$New.Old,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.exp = TRUE,show.summary = FALSE,tdcol.expected = "#993333",file = "Q7_New.Old.doc")
tab_xtab(rev2_wts$Gender,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q7_Gender.doc")
tab_xtab(rev2_wts$Age,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q7_Age.doc")
tab_xtab(rev2_wts$Children,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q7_Children.doc")
tab_xtab(rev2_wts$Race_Ethn,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q7_Race_Ethn.doc")
tab_xtab(rev2_wts$Income,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q7_Income.doc")
tab_xtab(rev2_wts$Coastal,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q7_Coastal.doc")

# Q10
tab_xtab(Q10_df2$Gender,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q10_Gender.doc")
tab_xtab(Q10_df2$Age,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q10_Age.doc")
tab_xtab(Q10_df2$Race_Ethn,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q10_Race_Ethn.doc")
tab_xtab(Q10_df2$Income,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q10_Income.doc")
tab_xtab(Q10_df2$SW_FW2,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q10_SW_FW.doc")
tab_xtab(Q10_df2$COVID,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q10_COVID.doc")
tab_xtab(Q10_df2$Children,Q10_df2$Agreement,weight.by = Q10_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q10_Children.doc")

# Q14
tab_xtab(Q14_df2$Gender,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q14_Gender.doc")
tab_xtab(Q14_df2$Age,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q14_Age.doc")
tab_xtab(Q14_df2$Race_Ethn,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q14_Race_Ethn.doc")
tab_xtab(Q14_df2$SW_FW2,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q14_SW_FW.doc")
tab_xtab(Q14_df2$Children,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q14_Children.doc")
tab_xtab(Q14_df2$COVID,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q14_COVID.doc")
tab_xtab(Q14_df2$New.Old,Q14_df2$Agreement,weight.by = Q14_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q14_New.Old.doc")

# Q24
tab_xtab(Q24_df2$Age,Q24_df2$Agreement,weight.by = Q24_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q24_Age.doc")
tab_xtab(Q24_df2$Children,Q24_df2$Agreement,weight.by = Q24_df2$weights,show.exp = TRUE,tdcol.expected = "#993333",show.summary = FALSE,file = "Q24_Children.doc")

tab_df(weights2) # Numeric weights and frequencies of age and gender

Age	Survey_Female	Database_Female	Weights_Female	Survey_Male	Database_Male	Weights_Male
18-24	2.07	5.29	2.55	2.71	13.66	5.04
25-34	7.33	8.07	1.10	6.59	14.21	2.16
35-44	9.09	8.52	0.94	12.59	14.85	1.18
45-54	8.55	6.46	0.75	16.21	11.57	0.71
55-64	8.13	3.97	0.49	13.50	7.33	0.54
65 or older	4.09	1.98	0.48	9.14	4.10	0.45

tab_dfs(list(Q5_summary3,Q11_summary3,Q18_19_summary3)) # Summary responses to likert scale questions

Variables	median	IQR
Food	1.24	1.48
Sport	1.86	1.87
Outdoors	2.67	2.02
Nature	2.58	2.02
Friends	2.54	1.82
Children1	2.04	2.15
Stress	2.51	1.90
Relax	2.58	1.85
Get away	1.91	2.13

 

Variables	median	IQR
Enjoyed fishing	3.25	1.98
Too expensive	1.67	1.41
Fishing rules confusing	1.63	1.23
Enjoyed surroundings	3.15	1.99

 

Variables	median	IQR
Motorboat access	2.58	2.12
Peaceful areas	3.36	2.54
Portable restrooms	3.09	2.17
Paddling trail	2.44	2.28
Fish stockings	3.15	2.25
Parking area	3.25	2.24
ADA accessible	2.26	1.72
Baby changing stations	1.91	1.88
Shoreline access	3.25	2.49
Canoe access	2.54	2.29
Space	3.34	2.49
Brick restrooms	2.69	2.29
Lighting	3.09	2.28
Picnic areas	2.69	2.16
Playground	2.14	1.99
Shade	3.15	2.25
Swimming	2.44	2.39
Camping	2.67	2.16

tab_dfs(list(Q7_chi_sq_SW_FW2,Q10_chi_sq,Q14_chi_sq,Q24_chi_sq)) # summary of significant chi-squared test results

Variable	df	F	p
SW/FW vs New/Retained	2	10.58	0.00
SW/FW vs Age	10	4.09	0.00
SW/FW vs Gender	2	4.43	0.01
SW/FW vs Children	4	5.69	0.00
SW/FW vs Race/Ethnicity	12	3.76	0.00
SW/FW vs Coastal	2	135.15	0.00

 

Variable	df	F	p
Age	25	21.67	0.00
Race/Ethnicity	30	1.97	0.00
Income	30	1.51	0.04
SW/FW	10	3.38	0.00
Gender	5	7.21	0.00
Children	10	37.24	0.00

 

Variable	df	F	p
COVID-Prompted	5	13.98	0.00
Age	25	4.91	0.00
Race/Ethnicity	30	2.03	0.00
New/Retained	5	3.04	0.01
SW/FW	10	3.66	0.00
Gender	5	5.86	0.00
Children	10	10.56	0.00

 

Variable	df	F	p
Age	40	2.09	0.00
Children	16	2.48	0.00

tab_dfs(list(Q5_envfit2,Q11_envfit2,Q18_19_envfit2)) # summary of significant RDA variables

Variables	RDA 1	RDA 2	r2	p
Gende.Female	-0.44	0.90	0.09	0.00
Gende.Male	0.44	-0.90	0.09	0.00
Age.18.24	-0.96	0.27	0.02	0.00
Age.25.34	-0.98	0.21	0.05	0.00
Age.65.or.older	0.81	-0.58	0.01	0.00
Child.18.or.over	1.00	-0.10	0.14	0.00
Child.No.children	-1.00	0.09	0.26	0.00
Child.Under.17	1.00	-0.05	0.02	0.00
COVID.Yes	-0.89	0.46	0.03	0.00
New.O.Old	0.67	-0.74	0.03	0.00

 

Variables	RDA 1	RDA 2	r2	p
Age.18.24	1.00	0.01	0.54	0.00
Age.25.34	1.00	0.01	0.08	0.00
Age.55.64	-1.00	-0.01	0.09	0.00
Age.65.or.older	-1.00	-0.01	0.05	0.00
Child.18.or.over	-1.00	-0.01	0.21	0.00
Child.No.children	1.00	0.02	0.12	0.00

 

Variables	RDA 1	RDA 2	r2	p
Incom.Over..150.000	0.98	-0.22	0.03	0.00
Age.25.34	-1.00	0.03	0.03	0.00
Age.35.44	-0.88	-0.48	0.02	0.00
Age.55.64	0.93	-0.36	0.02	0.00
Child.Under.17	-0.77	-0.64	0.09	0.00
New.O.New	-0.79	0.61	0.04	0.00
New.O.Old	0.79	-0.61	0.04	0.00
Race_.White	1.00	0.00	0.04	0.00
SW_FW.FW	-0.69	0.73	0.02	0.00
SW_FW.SW	0.52	-0.86	0.04	0.00

tab_model(Q4_mod3,Q0_mod2,Q3_mod3,Q9_mod2,Q8_mod2,wrap.labels = FALSE,show.r2 = FALSE,order.terms = c(1,25,10,9,7,8,17,18,2,3,4,5,6,11,12,13,14,15,16,19,20,21,22,23,24),pred.labels = var_names,dv.labels = c("Prompted by COVID-19 to Fish","New Angler (i.e. Not Retained)","Fished in Last 6 months","Fishing was the Main Trip Purpose","Fished from a Boat")) # binomial logistic regression results

	Prompted by COVID-19 to Fish			New Angler (i.e. Not Retained)			Fished in Last 6 months			Fishing was the Main Trip Purpose			Fished from a Boat
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	0.09	0.06 – 0.13	<0.001	1.65	1.23 – 2.22	0.001	7.37	5.53 – 9.81	<0.001	2.79	1.73 – 4.49	<0.001	2.06	1.35 – 3.14	0.001
ChildrenNo children	2.14	1.50 – 3.06	<0.001
ChildrenUnder 17	1.66	1.05 – 2.64	0.031
GenderFemale	1.89	1.33 – 2.69	<0.001	1.90	1.45 – 2.51	<0.001				0.73	0.53 – 1.01	0.055
New.OldNew	1.49	0.99 – 2.27	0.059							0.84	0.60 – 1.19	0.327	0.68	0.50 – 0.94	0.018
Race_EthnAsian or Pacific Islander				2.81	1.37 – 5.74	0.005				1.13	0.48 – 2.66	0.785	0.59	0.27 – 1.28	0.182
Race_EthnBlack or African American				1.86	1.00 – 3.49	0.052				2.25	1.00 – 5.08	0.050	0.20	0.10 – 0.39	<0.001
Race_EthnHispanic				0.89	0.28 – 2.88	0.850				2.59	0.65 – 10.29	0.178	0.49	0.20 – 1.19	0.115
Race_EthnNative American or Alaskan Native				1.28	0.39 – 4.23	0.688				4.28	0.93 – 19.64	0.061	1.98	0.78 – 5.01	0.151
Race_EthnOther Hispanic				0.71	0.35 – 1.43	0.339				3.47	1.39 – 8.67	0.008	0.28	0.14 – 0.58	<0.001
Race_EthnWhite Hispanic				1.20	0.80 – 1.81	0.382				1.40	0.86 – 2.30	0.178	0.61	0.41 – 0.92	0.019
SW_FW2SW				0.51	0.35 – 0.73	<0.001							1.91	1.32 – 2.77	0.001
SW_FW2SW.FW				0.71	0.51 – 0.98	0.038							2.20	1.55 – 3.11	<0.001
Age18-24				2.72	1.44 – 5.14	0.002	1.40	0.60 – 3.29	0.441	2.35	0.93 – 5.94	0.072
Age25-34				1.68	1.13 – 2.48	0.010	2.64	1.36 – 5.10	0.004	0.83	0.52 – 1.34	0.446
Age35-44				1.41	1.02 – 1.95	0.039	1.63	1.02 – 2.59	0.041	0.86	0.57 – 1.29	0.466
Age55-64				1.16	0.84 – 1.60	0.375	0.98	0.65 – 1.47	0.906	0.91	0.59 – 1.40	0.672
Age65 or older				1.08	0.73 – 1.60	0.696	0.40	0.26 – 0.59	<0.001	0.53	0.31 – 0.88	0.015
IncomeLess than $20,000										0.56	0.24 – 1.30	0.177	0.31	0.15 – 0.65	0.002
IncomeBetween $20,000 and $39,999										3.40	1.66 – 6.97	0.001	0.37	0.20 – 0.67	0.001
IncomeBetween $40,000 and $59,999										1.73	0.98 – 3.04	0.059	0.54	0.32 – 0.91	0.020
IncomeBetween $60,000 and $79,999										1.48	0.87 – 2.51	0.150	0.63	0.39 – 1.01	0.054
IncomeBetween $80,000 and $99,999										1.11	0.67 – 1.84	0.692	0.84	0.51 – 1.36	0.470
IncomeOver $150,000										0.76	0.48 – 1.21	0.249	1.04	0.66 – 1.65	0.864
CoastalY										1.51	1.09 – 2.09	0.014
COVIDYes													0.69	0.48 – 0.97	0.034
Observations	1314			1553			1901			1190			1238

tab_model(Q6_mod3,Q12_mod3,Q13_mod2,Q16_mod2,Q17_mod3,Q20_mod2,Q21_mod2,Q25_mod2,show.r2 = FALSE,pred.labels = var_names2,order.terms = c(12,23,3,34,35,1,2,4,5,6,7,8,13,14,15,16,17,18,24,25,26,27,28,29,9,10,11,19,20,21,22,30,31,32,33,36,37,38,39,40,41,42),dv.labels = c("How many days did you fish in the last 6 months?","How often do you catch fish?","How often do you keep fish?","How far did you travel to fish?","How far would you travel to fish?","How likely are you to purchase another license?","How likely would you be to enroll in an auto-renew license program?","Did they purchase a license in 2021 or 2022?")) # ordinal logistic regression results

	How many days did you fish in the last 6 months?			How often do you catch fish?			How often do you keep fish?			How far did you travel to fish?			How far would you travel to fish?			How likely are you to purchase another license?			How likely would you be to enroll in an auto-renew license program?			Did they purchase a license in 2021 or 2022?
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
SW_FW2SW	0.46	0.34 – 0.63	<0.001	1.19	0.87 – 1.63	0.271	3.45	2.54 – 4.68	<0.001	3.02	2.06 – 4.42	<0.001	1.47	1.01 – 2.14	0.046	0.67	0.38 – 1.16	0.154
SW_FW2SW.FW	1.64	1.25 – 2.15	<0.001	1.36	1.03 – 1.81	0.032	2.20	1.67 – 2.90	<0.001	2.25	1.59 – 3.17	<0.001	2.30	1.67 – 3.17	<0.001	1.88	1.06 – 3.35	0.031
GenderFemale	0.73	0.58 – 0.91	0.007																			0.69	0.54 – 0.88	0.002
Age18-24	1.36	0.83 – 2.22	0.226				0.84	0.51 – 1.39	0.507													1.01	0.66 – 1.56	0.956
Age25-34	1.24	0.89 – 1.71	0.201				0.94	0.69 – 1.27	0.676													1.09	0.79 – 1.51	0.586
Age35-44	0.98	0.74 – 1.29	0.868				1.19	0.91 – 1.57	0.210													0.99	0.76 – 1.28	0.931
Age55-64	0.89	0.68 – 1.18	0.433				1.42	1.07 – 1.89	0.016													0.99	0.75 – 1.32	0.963
Age65 or older	0.65	0.46 – 0.93	0.018				2.24	1.57 – 3.18	<0.001													3.33	2.27 – 4.90	<0.001
Race_EthnAsian or Pacific Islander	0.88	0.53 – 1.46	0.611	0.68	0.41 – 1.14	0.143	1.56	0.95 – 2.56	0.082	1.82	1.10 – 3.01	0.020							1.06	0.62 – 1.82	0.832
Race_EthnBlack or African American	0.52	0.31 – 0.88	0.014	0.49	0.30 – 0.78	0.003	2.52	1.09 – 5.84	0.031	1.09	0.73 – 1.62	0.688							0.51	0.27 – 0.96	0.037
Race_EthnHispanic	0.42	0.13 – 1.33	0.141	0.27	0.07 – 0.99	0.049	0.99	0.58 – 1.67	0.957	2.43	1.21 – 4.87	0.012							1.06	0.53 – 2.12	0.880
Race_EthnNative American or Alaskan Native	2.06	0.73 – 5.76	0.169	1.59	1.16 – 2.17	0.004	1.45	0.73 – 2.88	0.283	0.62	0.20 – 1.96	0.415							0.75	0.16 – 3.62	0.723
Race_EthnOther Hispanic	1.28	0.65 – 2.50	0.478	0.94	0.51 – 1.75	0.852	0.85	0.50 – 1.44	0.537	1.12	0.66 – 1.91	0.675							2.49	1.30 – 4.77	0.006
Race_EthnWhite Hispanic	0.87	0.63 – 1.21	0.408	0.46	0.33 – 0.65	<0.001	0.93	0.68 – 1.26	0.623	0.88	0.60 – 1.28	0.495							1.17	0.82 – 1.67	0.374
Once (1 day)|2–3 days	0.04	0.02 – 0.05	<0.001
2–3 days|4–10 days	0.29	0.22 – 0.39	<0.001
4–10 days|More than 10 days	1.27	0.97 – 1.66	0.084
COVIDYes				0.56	0.41 – 0.77	<0.001	0.63	0.48 – 0.84	0.002				0.72	0.53 – 0.97	0.029	0.35	0.22 – 0.58	<0.001				0.73	0.54 – 1.00	0.050
Never|Rarely				0.02	0.01 – 0.04	<0.001	0.43	0.32 – 0.58	<0.001
Rarely|Sometimes				0.10	0.08 – 0.13	<0.001	1.43	1.06 – 1.92	0.018
Sometimes|Most of the time				0.59	0.48 – 0.74	<0.001	6.93	4.98 – 9.66	<0.001
Most of the time|Always				6.03	4.64 – 7.83	<0.001	38.59	26.06 – 57.13	<0.001
New.OldNew							0.69	0.54 – 0.88	0.003				0.72	0.55 – 0.94	0.015	0.30	0.16 – 0.54	<0.001				0.56	0.44 – 0.72	<0.001
CoastalY										0.60	0.44 – 0.81	0.001	0.55	0.41 – 0.73	<0.001
IncomeLess than $20,000										0.36	0.18 – 0.74	0.006	0.38	0.14 – 0.99	0.048				0.51	0.23 – 1.10	0.085
IncomeBetween $20,000 and $39,999										0.42	0.25 – 0.72	0.002	0.88	0.54 – 1.45	0.626				0.47	0.28 – 0.79	0.005
IncomeBetween $40,000 and $59,999										0.94	0.62 – 1.42	0.756	1.07	0.71 – 1.60	0.750				0.67	0.45 – 1.01	0.053
IncomeBetween $60,000 and $79,999										0.68	0.46 – 1.01	0.058	1.07	0.69 – 1.65	0.774				1.01	0.70 – 1.46	0.964
IncomeBetween $80,000 and $99,999										0.88	0.57 – 1.35	0.564	0.82	0.59 – 1.15	0.256				0.94	0.62 – 1.42	0.757
IncomeOver $150,000										1.10	0.74 – 1.62	0.638	1.39	0.99 – 1.96	0.056				1.57	1.10 – 2.24	0.013
Less than 30 minutes|Between 30 and 60 minutes										0.60	0.43 – 0.83	0.002	0.02	0.01 – 0.04	<0.001
Between 30 and 60 minutes|Between 1 and 2 hours										2.13	1.52 – 2.99	<0.001	0.23	0.16 – 0.34	<0.001
Between 1 and 2 hours|More than 2 hours										7.00	4.91 – 9.98	<0.001
Between 1 and 2 hours|Between 2 and 4 hours													0.91	0.63 – 1.30	0.606
Between 2 and 4 hours|More than 4 hours													2.79	1.94 – 4.02	<0.001
Not likely|Somewhat likely																0.00	0.00 – 0.01	<0.001
Somewhat likely|Very likely																0.02	0.01 – 0.04	<0.001
Extremely unlikely|Unlikely																			0.00	0.00 – 0.01	<0.001
Unlikely|Not sure																			0.05	0.04 – 0.07	<0.001
Not sure|Likely																			0.32	0.24 – 0.42	<0.001
Likely|Extremely likely																			1.20	0.92 – 1.55	0.178
Q20_NewLicenseNot likely																						1.72	0.35 – 8.40	0.500
Q20_NewLicenseVery likely																						4.23	2.49 – 7.18	<0.001
None|One_Yr																						0.50	0.28 – 0.91	0.024
One_Yr|Both_Yrs																						3.88	2.13 – 7.09	<0.001
Observations	1553			1565			1527			1190			1248			1717			1386			1608