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

This is an analysis of a survey of 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)

# 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)

# Clean independent variables
rev1 <- raw3 %>%
  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),
         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 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

## Adjust responses using weights
rev1_wts <- rev1 %>%
  left_join(weights1, by = c("Gender","Age")) %>%
  replace_na(list(weights = 1)) %>%
  mutate(Gender = as.factor(Gender))

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.

Other demographics of survey respondents

# Don't include Age & Gender factor b/c this is given above
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 <- rev2_wts %>%
  filter(!(is.na(Income))) %>% 
  count(Income) %>%
  mutate(freq = prop.table(n)) # get proportions

Children <- rev2_wts %>%
  filter(!(is.na(Children))) %>% 
  count(Children) %>%
  mutate(freq = prop.table(n)) # get proportions

Race_Ethn <- rev2_wts %>%
  filter(!(is.na(Race_Ethn))) %>% 
  count(Race_Ethn) %>%
  mutate(freq = prop.table(n)) # get proportions

New.Old <- rev2_wts %>%
  filter(!(is.na(New.Old))) %>% 
  count(New.Old) %>%
  mutate(freq = prop.table(n)) # get proportions

COVID <- rev2_wts %>%
  filter(!(is.na(COVID))) %>% 
  count(COVID) %>%
  mutate(freq = prop.table(n)) # get proportions

Children %>% 
  kbl(caption = "Children",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Children

Children	n	freq
18 or over	852	0.56
No children	501	0.33
Under 17	182	0.12

Income %>% 
  kbl(caption = "Income",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Income

Income	n	freq
Less than $20,000	70	0.05
Between $20,000 and $39,999	142	0.10
Between $40,000 and $59,999	186	0.13
Between $60,000 and $79,999	226	0.15
Between $80,000 and $99,999	211	0.14
Between $100,000 and $149,999	332	0.23
Over $150,000	298	0.20

Race_Ethn %>% 
  kbl(caption = "Race and Ethnicity",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Race and Ethnicity

Race_Ethn	n	freq
Asian or Pacific Islander	64	0.04
Black or African American	84	0.05
Hispanic	41	0.02
Native American or Alaskan Native	40	0.02
Other Hispanic	74	0.04
White	1220	0.69
White Hispanic	242	0.14

New.Old %>% 
  kbl(caption = "Retained or New Angler during pandemic",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

Retained or New Angler during pandemic

New.Old	n	freq
New	1517	0.68
Old	721	0.32

COVID %>% 
  kbl(caption = "COVID-Prompted Angler",digits = 2) %>%
  kable_classic(full_width = F, html_font = "Cambria") 

COVID-Prompted Angler

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

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). Most anglers (81%) were already planning on going fishing and were not prompted by COVID-19 to go fishing. Approximately half of the anglers surveyed (56%) had adult children, 12% had children under 17 yrs old, and 32% had no children. There were approximately twice as many freshwater anglers (42%) compared to saltwater anglers (24%), while 34% of anglers reported that they fish in both freshwater and saltwater. The survey respondents were primarily composed of anglers that identified as White (69%), followed by White Hispanic (14%), Black or African American (5%), Other Hispanic (4%), Asian or Pacific Islander (4%), Hispanic (2%), and Native American or Alaskan Native (2%). Finally, the majority of survey respondents had incomes over $80K (57%), while 38% of respondents had incomes between $20K-$79K, and only 5% had incomes less than $20K.

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

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 <- 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"))
rev2_wts$Income <- relevel(rev2_wts$Income, ref = "Less than $20,000")
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, 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)) 

# Reformat model data for figure
Q0_df1 <- model.frame(Q0_mod2) %>% 
  group_by(New.Old,Age,Gender,Race_Ethn,SW_FW2) %>%
  count() %>%
  na.omit()
Q0_freq = prop.table(Q0_df1$n)
Q0_df1 <- data.frame(cbind(Q0_df1,freq=Q0_freq))

# 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

# Barplot of responses
ggplot(Q0_df1,aes(x=New.Old,y=freq,fill = Age)) + 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 = "Q0. Fished due to COVID") + facet_wrap(~Gender, ncol = 2)

ggplot(Q0_df1,aes(x=New.Old,y=freq,fill = Race_Ethn)) + 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 = "Q0. Fished due to COVID") + facet_wrap(~SW_FW2, ncol = 3)

# 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)	3.233	0.494	6.538	0.000	25.352	9.612	66.868
GenderMale	-0.654	0.139	-4.701	0.000	0.520	0.396	0.683
Race_EthnBlack or African American	-0.441	0.473	-0.932	0.351	0.644	0.255	1.627
Race_EthnHispanic	-1.107	0.663	-1.669	0.095	0.331	0.090	1.214
Race_EthnNative American or Alaskan Native	-0.891	0.710	-1.254	0.210	0.410	0.102	1.653
Race_EthnOther Hispanic	-1.381	0.490	-2.818	0.005	0.251	0.096	0.657
Race_EthnWhite	-1.067	0.365	-2.925	0.003	0.344	0.168	0.704
Race_EthnWhite Hispanic	-0.868	0.396	-2.194	0.028	0.420	0.193	0.912
SW_FW2SW	-0.673	0.183	-3.669	0.000	0.510	0.356	0.731
SW_FW2SW.FW	-0.336	0.166	-2.029	0.043	0.714	0.516	0.989
Age25-34	-0.537	0.348	-1.544	0.123	0.585	0.296	1.156
Age35-44	-0.681	0.329	-2.072	0.038	0.506	0.265	0.964
Age45-54	-1.035	0.325	-3.188	0.001	0.355	0.188	0.672
Age55-64	-0.882	0.328	-2.686	0.007	0.414	0.217	0.788
Age65 or older	-0.960	0.345	-2.783	0.005	0.383	0.195	0.753

# Effect plot
plot(Effect(focal.predictors = c("Age","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.

New anglers recruited during the pandemic were: 1.923 times more likely to be female, 1.960 times more likely to fish in freshwater compared to saltwater, 1.400 times more likely to fish in freshwater compared to fishing in both saltwater and freshwater, 2.382-3.981 times more likely to identify as Asian or Pacific Islander compared to White Hispanic, White, or Other Hispanic, and 1.977 - 2.612 times more likely to be between the ages of 18-24 compared to 35 years or older.

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?

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"))
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, 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

# Reformat model data for figure
Q3_df1 <- model.frame(Q3_mod3) %>% 
  group_by(Q3_Fished,Age) %>%
  count() %>%
  na.omit()
Q3_freq = prop.table(Q3_df1$n)
Q3_df1 <- data.frame(cbind(Q3_df1,freq=Q3_freq))

# 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

# Histogram of responses
ggplot(Q3_df1,aes(x=Q3_Fished,y=freq,fill = Age)) + 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 = "Q3. Fished in last 6 months")

# 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)	2.380	0.407	5.845	0.000	10.800	4.860	23.996
Age25-34	0.584	0.506	1.155	0.248	1.793	0.665	4.831
Age35-44	0.085	0.447	0.189	0.850	1.088	0.453	2.616
Age45-54	-0.357	0.432	-0.826	0.409	0.700	0.300	1.634
Age55-64	-0.407	0.433	-0.940	0.347	0.665	0.284	1.557
Age65 or older	-1.298	0.432	-3.001	0.003	0.273	0.117	0.638

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

Most respondents (88%) went fishing between mid-March 2020 and 2021. Younger individuals (18-24 yrs) were 3.66 times more likely than older individuals (65+) to have fished in the last 6 months.

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

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"))
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, 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 ~ Age + 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"))

# Reformat model data for figure
Q4_df1 <- model.frame(Q4_mod3) %>% 
  group_by(Q4_COVID_Fished,Age,Gender,New.Old,Children) %>%
  count() %>%
  na.omit()
Q4_freq = prop.table(Q4_df1$n)
Q4_df1 <- data.frame(cbind(Q4_df1,freq=Q4_freq))

# 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

Q4_AOR_Male <- 1/0.52 # 1.92 x

# 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

# Barplot of responses
ggplot(Q4_df1,aes(x=Q4_COVID_Fished,y=freq,fill = Age)) + 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 = "Q4. Fished due to COVID") + facet_wrap(~Gender, ncol = 2)

ggplot(Q4_df1,aes(x=Q4_COVID_Fished,y=freq,fill = Children)) + 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 = "Q4. Fished due to COVID") + facet_wrap(~New.Old, ncol = 2)

# 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.212	0.383	-5.775	0.000	0.109	0.052	0.232
Age25-34	0.721	0.334	2.161	0.031	2.057	1.069	3.960
Age35-44	0.562	0.332	1.691	0.091	1.753	0.914	3.364
Age45-54	0.469	0.346	1.354	0.176	1.599	0.810	3.154
Age55-64	0.492	0.341	1.445	0.149	1.636	0.839	3.193
Age65 or older	-0.044	0.399	-0.111	0.911	0.957	0.437	2.092
ChildrenNo children	0.753	0.189	3.990	0.000	2.123	1.466	3.074
ChildrenUnder 17	0.343	0.268	1.283	0.200	1.410	0.834	2.384
GenderMale	-0.604	0.178	-3.400	0.001	0.547	0.386	0.775
New.OldNew	0.424	0.213	1.988	0.047	1.528	1.006	2.322

# Effect plot
plot(Effect(focal.predictors = c("Children","Age","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.02 times more likely for young anglers (25-34) compared to other ages, 2.15 times more likely for anglers with no children compared to those with adult children, 1.92 times more likely for female anglers compared to males, 1.53 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$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"))
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),all_vars(!is.na(.))) 
Q5_df2 <- Q5_df1 %>%
  dplyr::select(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2,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)
# 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.No..I.was.already.planning.to.go.fishing.,COVID.Yes,New.O.New,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=c("Female","Age 18-24","Age 25-34","Age 65 over","Children 18 over","No Children","COVID-No","New"),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
# Recode factors 1-5 for comparison to other questions in final table
Q5_summary1 <- rev2_wts %>%
  dplyr::select(Food,Sport,Outdoors,Nature,Friends,Children1,Stress,Relax,Get_away)

tab_itemscale(Q5_summary1)

Q5_summary2 <- describe(Q5_summary1,na.rm = TRUE) 
Q5_summary2$Variables <- row.names(Q5_summary2) 
Q5_summary3 <- Q5_summary2[,c("Variables","mean","sd")]
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	mean	sd
Food	1.82	1.91
Sport	2.46	2.43
Outdoors	3.42	2.82
Nature	3.41	2.85
Friends	3.22	2.68
Children1	2.52	2.08
Stress	3.16	2.77
Relax	3.40	2.86
Get away	2.55	2.56

# 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)"))
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.No..I.was.already.planning.to.go.fishing. + COVID.Yes + New.O.New + New.O.Old, data = Q5_indep3)
##           Df Variance       F Pr(>F)    
## RDA1       1 0.005123 62.9800  0.001 ***
## RDA2       1 0.001495 18.3818  0.001 ***
## RDA3       1 0.000435  5.3462  0.006 ** 
## RDA4       1 0.000222  2.7263  0.299    
## RDA5       1 0.000114  1.3991  0.867    
## RDA6       1 0.000033  0.4019  1.000    
## RDA7       1 0.000015  0.1790  1.000    
## RDA8       1 0.000001  0.0075  1.000    
## Residual 950 0.077277                   
## ---
## 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.43	0.90	0.09	0.00
Gende.Male	0.43	-0.90	0.09	0.00
Age.18.24	-0.98	0.22	0.02	0.00
Age.25.34	-0.98	0.22	0.05	0.00
Age.65.or.older	0.81	-0.59	0.02	0.00
Child.18.or.over	1.00	-0.10	0.13	0.00
Child.No.children	-1.00	0.08	0.26	0.00
Child.Under.17	1.00	-0.03	0.03	0.00
COVID.No..I.was.already.planning.to.go.fishing.	0.89	-0.45	0.02	0.00
COVID.Yes	-0.89	0.45	0.02	0.00
New.O.New	-0.65	0.76	0.03	0.00
New.O.Old	0.65	-0.76	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 anglers not prompted by COVID to go fishing (COVID-No) were more likely to fish for food or sport. Anglers between the ages of 18-34 were likely to go fishing to enjoy nature, the outdoors, and to relieve stress. While anglers that recruited during COVID and female anglers were more likely to fish to spend time with friends. 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.9688616

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

# Order the factor
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"))
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, 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
Q6_mod4 <- svyolr(Q6_DaysFished ~ SW_FW2 + Gender + Age, design=svy_desn_Q6)
summary(Q6_mod4) # examine model results
VIF(Q6_mod4) # collinearity
FW_AOR <- 1/0.441 # 2.268x
Angler_65over <- 1/0.498 # 2.008x

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

# Reformat model data for figure
Q6_df1 <- model.frame(Q6_mod3) %>% 
  group_by(Q6_DaysFished,Age,Gender,SW_FW2) %>%
  count() %>%
  na.omit()
Q6_freq = prop.table(Q6_df1$n)
Q6_df1 <- data.frame(cbind(Q6_df1,freq=Q6_freq))

# 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

# Histogram of responses
ggplot(Q6_df1,aes(x=Q6_DaysFished,y=freq,fill = Age)) + 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 = "Q6. How many days have you fished?") + facet_wrap(~Gender, ncol = 2)

ggplot(Q6_df1,aes(x=Q6_DaysFished,y=freq,fill = SW_FW2)) + 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 = "Q6. How many days have you fished?")

# 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.82	0.15	-5.45	0.00	0.44
SW_FW2SW.FW	0.52	0.13	3.94	0.00	1.68
GenderMale	0.39	0.11	3.40	0.00	1.47
Age25-34	-0.13	0.25	-0.52	0.60	0.88
Age35-44	-0.33	0.23	-1.43	0.15	0.72
Age45-54	-0.34	0.23	-1.49	0.14	0.71
Age55-64	-0.45	0.24	-1.91	0.06	0.64
Age65 or older	-0.70	0.26	-2.69	0.01	0.50
Once (1 day)|2–3 days	-3.22	0.29	-11.15	0.00	0.04
2–3 days|4–10 days	-1.15	0.24	-4.71	0.00	0.32
4–10 days|More than 10 days	0.32	0.24	1.34	0.18	1.38

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

The odds of fishing the most often (>10 days) were: 1.47 times more likely for males compared to females, 1.68 times more likely for anglers that fish in both freshwater and saltwater compared to just freshwater anglers, 2.27 times more likely for freshwater anglers compared to saltwater anglers, and 2.01 times more likely for young anglers (18-24 yrs) compared to the oldest anglers (65 yrs+).

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

# 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)

# 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)

# 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_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))
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	193 10.1 %	156 8.1 %	230 12 %	579 30.2 %
New	611 31.9 %	248 12.9 %	480 25 %	1339 69.8 %
Total	804 41.9 %	404 21.1 %	710 37 %	1918 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.305, ndf = 2, ddf = 2236, p-value = 3.508e-05

Q7_Gender_xtab

Gender	SW_FW2			Total
	FW	SW	SW.FW
Female	264 15.4 %	138 8.1 %	188 11 %	590 34.5 %
Male	464 27.1 %	210 12.3 %	447 26.1 %	1121 65.5 %
Total	728 42.5 %	348 20.3 %	635 37.1 %	1711 100 %

Q7_Gender_SW_FW2_chisq

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

Q7_Age_xtab

Age	SW_FW2			Total
	FW	SW	SW.FW
18-24	137 8 %	46 2.7 %	140 8.2 %	323 18.9 %
25-34	167 9.8 %	64 3.7 %	168 9.8 %	399 23.3 %
35-44	169 9.9 %	85 5 %	154 9 %	408 23.9 %
45-54	124 7.2 %	73 4.3 %	104 6.1 %	301 17.6 %
55-64	85 5 %	56 3.3 %	54 3.2 %	195 11.5 %
65 or older	45 2.6 %	21 1.2 %	19 1.1 %	85 4.9 %
Total	727 42.5 %	345 20.2 %	639 37.3 %	1711 100 %

Q7_Age_SW_FW2_chisq

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

Q7_Children_xtab

Children	SW_FW2			Total
	FW	SW	SW.FW
18 or over	211 15.4 %	136 9.9 %	150 10.9 %	497 36.2 %
No children	262 19.1 %	104 7.6 %	278 20.2 %	644 46.9 %
Under 17	105 7.6 %	45 3.3 %	83 6 %	233 16.9 %
Total	578 42.1 %	285 20.7 %	511 37.2 %	1374 100 %

Q7_Children_SW_FW2_chisq

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

Q7_Race_Ethn_xtab

Race_Ethn	SW_FW2			Total
	FW	SW	SW.FW
Asian or Pacific Islander	25 1.6 %	11 0.7 %	25 1.6 %	61 3.9 %
Black or African American	33 2.1 %	7 0.4 %	27 1.7 %	67 4.2 %
Hispanic	6 0.4 %	13 0.8 %	20 1.2 %	39 2.4 %
Native American or Alaskan Native	16 1 %	10 0.6 %	18 1.1 %	44 2.7 %
Other Hispanic	37 2.3 %	22 1.4 %	49 3.1 %	108 6.8 %
White	480 30 %	180 11.2 %	329 20.5 %	989 61.7 %
White Hispanic	86 5.4 %	87 5.4 %	120 7.5 %	293 18.3 %
Total	683 42.7 %	330 20.6 %	588 36.7 %	1601 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.3973, ndf = 12, ddf = 2226, p-value = 5.964e-05

tab_df(Q7_chi_sq_SW_FW2)

Variable	df	F	p
SW/FW vs New/Retained	2	10.31	0.00
SW/FW vs Age	10	4.17	0.00
SW/FW vs Gender	2	4.24	0.01
SW/FW vs Children	4	5.42	0.00
SW/FW vs Race/Ethnicity	12	3.40	0.00

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

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"))
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, design=svy_desn_Q8, family=quasibinomial)
summary(Q8_mod1) # examine model results
varImp(Q8_mod1) # Variable importance
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)
varImp(Q8_mod2) # Variable importance
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)) 

# Reformat model data for figure
Q8_df1 <- model.frame(Q8_mod2) %>% 
  group_by(Boat_Y_N,Income,Race_Ethn,SW_FW2,COVID,New.Old) %>%
  count() %>%
  na.omit()
Q8_freq = prop.table(Q8_df1$n)
Q8_df1 <- data.frame(cbind(Q8_df1,freq=Q8_freq))

# 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

# Barplot of responses
ggplot(Q8_df1,aes(x=Boat_Y_N,y=freq,fill = Income)) + 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 = "Q8. Did you fish from a boat?") 

ggplot(Q8_df1,aes(x=Boat_Y_N,y=freq,fill = Race_Ethn)) + 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 = "Q8. Did you fish from a boat?") 

ggplot(Q8_df1,aes(x=Boat_Y_N,y=freq,fill = COVID)) + 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 = "Q8. Did you fish from a boat?") + facet_wrap(~SW_FW2,ncol = 3)

# 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)	-1.020	0.506	-2.017	0.044	0.360	0.134	0.972
IncomeBetween $20,000 and $39,999	0.188	0.429	0.439	0.661	1.207	0.520	2.800
IncomeBetween $40,000 and $59,999	0.559	0.398	1.403	0.161	1.749	0.800	3.820
IncomeBetween $60,000 and $79,999	0.695	0.380	1.829	0.068	2.005	0.951	4.227
IncomeBetween $80,000 and $99,999	0.940	0.386	2.434	0.015	2.559	1.200	5.458
IncomeBetween $100,000 and $149,999	1.141	0.373	3.059	0.002	3.131	1.506	6.512
IncomeOver $150,000	1.187	0.382	3.103	0.002	3.276	1.547	6.937
Race_EthnBlack or African American	-0.918	0.508	-1.806	0.071	0.399	0.147	1.082
Race_EthnHispanic	-0.112	0.578	-0.193	0.847	0.894	0.288	2.782
Race_EthnNative American or Alaskan Native	1.275	0.593	2.150	0.032	3.580	1.118	11.462
Race_EthnOther Hispanic	-0.729	0.514	-1.418	0.156	0.482	0.176	1.323
Race_EthnWhite	0.615	0.395	1.557	0.120	1.849	0.852	4.013
Race_EthnWhite Hispanic	0.087	0.425	0.206	0.837	1.091	0.475	2.510
SW_FW2SW	0.663	0.190	3.480	0.001	1.940	1.335	2.819
SW_FW2SW.FW	0.781	0.176	4.433	0.000	2.183	1.545	3.084
COVIDYes	-0.370	0.176	-2.101	0.036	0.691	0.489	0.976
New.OldNew	-0.378	0.162	-2.325	0.020	0.685	0.498	0.943

# 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: 2.56 - 3.28 times more likely for anglers with an income between $80-$150K compared to anglers with incomes <$20K, 3.58 times more likely for Native American or Alaskan Native anglers compared to Asian or Pacific Islander anglers, 1.94 times more likely for anglers that fish in saltwater compared to those that fish in freshwater, 2.18 times more likely for anglers that fish in both saltwater and freshwater compared to those that fish in freshwater, 1.45 times more likely for anglers that were not prompted by COVID-19 to start fishing, and 1.46 times more likely for retained anglers compared to new 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$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"))
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, 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, design=svy_desn_Q9, family=quasibinomial)
summary(Q9_mod2)
varImp(Q9_mod2) # Variable importance
VIF(Q9_mod2) # collinearity
car::Anova(Q9_mod2, type = 3, test.statistic = "F",error.df = degf(Q9_mod2$survey.design)) # Evaluate variable significance

# Check model accuracy
logitgof(obs = Q9_mod2$y, fitted(Q9_mod2)) #run hoslem test to test goodness of fit
glm.probs <- predict(Q9_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 == Q9_mod2$y)) # Correct assignment 75% of the time 
Q9_age_xtab <- xtabs(~Age + Gender+Q9_Purpose, data = rev2_wts)
Q9_age_xtab <- Q9_age_xtab%>%ftable(row.vars=c("Gender","Age"))

# Estimate and odds ratio table
Q9_coef_table <- coef(summary(Q9_mod2))
Q9_AOR <- standardize_parameters(Q9_mod2, 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

# Reformat model data for figure
Q9_df1 <- model.frame(Q9_mod2) %>% 
  group_by(Q9_Purpose,Age,Gender,Income) %>%
  count() %>%
  na.omit()
Q9_freq = prop.table(Q9_df1$n)
Q9_df1 <- data.frame(cbind(Q9_df1,freq=Q9_freq))

Q9_AOR_65 <- 1/0.21 # 4.76 x
Q9_AOR_45 <- 1/0.41 # 2.44 x

# 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

# Histogram of responses
ggplot(Q9_df1,aes(x=Q9_Purpose,y=freq,fill = Age)) + 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 = "Q9. Was Fishing the Main Purpose?") + facet_wrap(~Gender, ncol = 2)

ggplot(Q9_df1,aes(x=Q9_Purpose,y=freq,fill = Income)) + 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 = "Q9. Was Fishing the Main Purpose?") 

# 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)	1.14	0.57	2.00	0.05	3.14	1.02	9.63
GenderMale	0.42	0.16	2.72	0.01	1.53	1.12	2.07
Age25-34	-0.96	0.44	-2.20	0.03	0.38	0.16	0.90
Age35-44	-0.99	0.43	-2.30	0.02	0.37	0.16	0.86
Age45-54	-0.89	0.44	-2.03	0.04	0.41	0.17	0.97
Age55-64	-0.95	0.44	-2.16	0.03	0.39	0.16	0.92
Age65 or older	-1.57	0.46	-3.45	0.00	0.21	0.08	0.51
IncomeBetween $20,000 and $39,999	1.78	0.49	3.64	0.00	5.92	2.27	15.42
IncomeBetween $40,000 and $59,999	1.13	0.42	2.66	0.01	3.09	1.34	7.10
IncomeBetween $60,000 and $79,999	1.03	0.42	2.45	0.01	2.81	1.23	6.41
IncomeBetween $80,000 and $99,999	0.65	0.41	1.59	0.11	1.91	0.86	4.26
IncomeBetween $100,000 and $149,999	0.63	0.40	1.58	0.11	1.89	0.86	4.15
IncomeOver $150,000	0.25	0.42	0.61	0.54	1.29	0.57	2.91

# Effect plot
plot(Effect(focal.predictors = c("Gender","Age"),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.53 times more likely for male anglers compared to females, 2.44 - 4.76 times more likely for anglers ages 18-24 compared to older anglers, 2.81 - 5.92 times more likely for anglers with incomes between $20-$79K compared to anglers with incomes <$20K.

Q10. Who do you typically fish with?

# 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"))
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,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,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") # 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_COVID_xtab <- tab_xtab(Q10_df2$COVID,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)

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

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

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

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

Q10_corr_COVID <- corrplot(Q10_COVID_chisq$residuals, is.cor = FALSE)

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

# 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_COVID <- data.frame(cbind("COVID-Prompted",Q10_COVID_chisq[[c(2,1)]],Q10_COVID_chisq$statistic,Q10_COVID_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_COVID,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
Female	38 1.5 %	159 6.1 %	239 9.1 %	281 10.7 %	99 3.8 %	98 3.7 %	914 34.9 %
Male	184 7 %	265 10.1 %	350 13.4 %	590 22.6 %	177 6.8 %	135 5.2 %	1701 65.1 %
Total	222 8.5 %	424 16.2 %	589 22.5 %	871 33.3 %	276 10.6 %	233 8.9 %	2615 100 %

Q10_Gender_chisq

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

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	51 1.9 %	10 0.4 %	120 4.6 %	231 8.8 %	10 0.4 %	21 0.8 %	443 16.9 %
25-34	57 2.2 %	125 4.8 %	143 5.5 %	239 9.1 %	24 0.9 %	5 0.2 %	593 22.7 %
35-44	39 1.5 %	185 7.1 %	143 5.5 %	169 6.5 %	120 4.6 %	29 1.1 %	685 26.3 %
45-54	37 1.4 %	69 2.6 %	88 3.4 %	132 5 %	92 3.5 %	67 2.6 %	485 18.5 %
55-64	30 1.1 %	24 0.9 %	66 2.5 %	70 2.7 %	23 0.9 %	69 2.6 %	282 10.7 %
65 or older	10 0.4 %	13 0.5 %	32 1.2 %	30 1.1 %	6 0.2 %	37 1.4 %	128 4.8 %
Total	224 8.6 %	426 16.3 %	592 22.6 %	871 33.3 %	275 10.5 %	228 8.7 %	2616 100 %

Q10_Age_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Age + Agreement, svy_desn_Q10, statistic = "adjWald")
## F = 22.169, 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
Asian or Pacific Islander	9 0.4 %	11 0.5 %	20 0.8 %	35 1.4 %	6 0.2 %	2 0.1 %	83 3.4 %
Black or African American	8 0.3 %	18 0.7 %	19 0.8 %	33 1.4 %	4 0.2 %	4 0.2 %	86 3.6 %
Hispanic	1 0 %	14 0.6 %	22 0.9 %	21 0.9 %	9 0.4 %	6 0.2 %	73 3 %
Native American or Alaskan Native	9 0.4 %	14 0.6 %	12 0.5 %	24 1 %	10 0.4 %	7 0.3 %	76 3.2 %
Other Hispanic	15 0.6 %	33 1.4 %	29 1.2 %	60 2.5 %	25 1 %	20 0.8 %	182 7.5 %
White	135 5.6 %	236 9.7 %	342 14.1 %	474 19.5 %	160 6.6 %	141 5.8 %	1488 61.3 %
White Hispanic	32 1.3 %	60 2.5 %	104 4.3 %	171 7 %	40 1.6 %	36 1.5 %	443 18.2 %
Total	209 8.6 %	386 15.9 %	548 22.5 %	818 33.6 %	254 10.4 %	216 8.9 %	2431 100 %

Q10_Race_Ethn_chisq

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

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	21 1 %	12 0.6 %	39 1.9 %	36 1.8 %	7 0.3 %	4 0.2 %	119 5.8 %
Between $20,000 and $39,999	27 1.3 %	42 2.1 %	59 2.9 %	94 4.6 %	14 0.7 %	14 0.7 %	250 12.3 %
Between $40,000 and $59,999	35 1.7 %	42 2.1 %	62 3 %	88 4.3 %	22 1.1 %	17 0.8 %	266 13 %
Between $60,000 and $79,999	31 1.5 %	63 3.1 %	68 3.3 %	117 5.7 %	29 1.4 %	27 1.3 %	335 16.3 %
Between $80,000 and $99,999	25 1.2 %	52 2.6 %	59 2.9 %	85 4.2 %	36 1.8 %	29 1.4 %	286 14.1 %
Between $100,000 and $149,999	28 1.4 %	69 3.4 %	91 4.5 %	137 6.7 %	59 2.9 %	43 2.1 %	427 21 %
Over $150,000	27 1.3 %	63 3.1 %	71 3.5 %	105 5.2 %	49 2.4 %	38 1.9 %	353 17.4 %
Total	194 9.5 %	343 16.8 %	449 22.1 %	662 32.5 %	216 10.6 %	172 8.4 %	2036 100 %

Q10_Income_chisq

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

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	131 4.6 %	218 7.6 %	253 8.8 %	345 12.1 %	107 3.7 %	92 3.2 %	1146 40 %
SW	27 0.9 %	84 2.9 %	137 4.8 %	209 7.3 %	68 2.4 %	68 2.4 %	593 20.7 %
SW.FW	77 2.7 %	172 6 %	251 8.8 %	392 13.7 %	127 4.4 %	101 3.5 %	1120 39.1 %
Total	235 8.2 %	474 16.6 %	641 22.4 %	946 33.1 %	302 10.6 %	261 9.1 %	2859 100 %

Q10_SW_FW2_chisq

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

Q10_COVID_xtab 

COVID	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
No, I was already planning to go fishing.	181 6.3 %	387 13.5 %	532 18.5 %	765 26.6 %	253 8.8 %	230 8 %	2348 81.7 %
Yes	54 1.9 %	88 3.1 %	116 4 %	186 6.5 %	49 1.7 %	33 1.1 %	526 18.3 %
Total	235 8.2 %	475 16.5 %	648 22.5 %	951 33.1 %	302 10.5 %	263 9.2 %	2874 100 %

Q10_COVID_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~COVID + Agreement, svy_desn_Q10, statistic = "adjWald")
## F = 2.3068, ndf = 5, ddf = 2854, p-value = 0.04205

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 %	83 4.1 %	182 9 %	187 9.3 %	132 6.5 %	179 8.9 %	810 40.1 %
No children	126 6.2 %	13 0.6 %	207 10.3 %	425 21.1 %	18 0.9 %	30 1.5 %	819 40.6 %
Under 17	24 1.2 %	118 5.9 %	85 4.2 %	119 5.9 %	28 1.4 %	14 0.7 %	388 19.3 %
Total	197 9.8 %	214 10.6 %	474 23.5 %	731 36.2 %	178 8.8 %	223 11.1 %	2017 100 %

Q10_Children_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Children + Agreement, svy_desn_Q10, statistic = "adjWald")
## F = 39.008, 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 likley 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$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"))
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),all_vars(!is.na(.))) 
# Select independent variables
Q11_indep1 <- Q11_df1 %>%
  dplyr::select(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2)
# 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, 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, SW_FW.SW,New.O.New, New.O.Old)
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,New.O.New, New.O.Old)
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=c("Age 18-24","Age 25-34","Age 55-64","Age 65 over","Children 18 over","No Children","New Angler"),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) 
Q11_summary2$Variables <- row.names(Q11_summary2) 
Q11_summary3 <- Q11_summary2[,c("Variables","mean","sd")]
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	mean	sd
Enjoyed fishing	3.91	3.12
Too expensive	2.13	1.96
Fishing rules confusing	1.98	1.64
Enjoyed surroundings	3.82	3.02

# 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 + New.O.New + New.O.Old, data = Q11_indep4)
##           Df Variance        F Pr(>F)    
## RDA1       1  1.91742 880.6484  0.001 ***
## RDA2       1  0.00362   1.6626  0.959    
## RDA3       1  0.00130   0.5950  0.998    
## RDA4       1  0.00053   0.2452  1.000    
## Residual 954  2.07713                    
## ---
## 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.02	0.58	0.00
Age.25.34	1.00	-0.03	0.07	0.00
Age.55.64	-1.00	0.02	0.10	0.00
Age.65.or.older	-1.00	0.02	0.05	0.00
Child.18.or.over	-1.00	0.02	0.20	0.00
Child.No.children	1.00	-0.02	0.12	0.00
New.O.New	1.00	0.02	0.01	0.00
New.O.Old	-1.00	-0.02	0.01	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")) +
  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 that recruited during the COVID-19 pandemic, those without children and between the ages of 18-34 were more likely to enjoy their surroundings and enjoy fishing, while those with children and 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.9658277

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

# Order the factor
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"))
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, 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

# Reformat model data for figure
Q12_df1 <- model.frame(Q12_mod3) %>% 
  group_by(Q12_CatchOften,SW_FW2,COVID,Race_Ethn) %>%
  count() %>%
  na.omit()
Q12_freq = prop.table(Q12_df1$n)
Q12_df1 <- data.frame(cbind(Q12_df1,freq=Q12_freq))

Q12_COVID_AOR <- 1/0.564 # 1.773x

# 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

# Histogram of responses
ggplot(Q12_df1,aes(x=Q12_CatchOften,y=freq,fill = COVID)) + 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 = "Q12. How often do you catch fish?") + facet_wrap(~SW_FW2, ncol = 3)

ggplot(Q12_df1,aes(x=Q12_CatchOften,y=freq,fill = Race_Ethn)) + 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 = "Q12. How often do you catch fish?") 

# 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.20	0.16	1.22	0.22	1.22
SW_FW2SW.FW	0.32	0.14	2.19	0.03	1.37
COVIDYes	-0.57	0.16	-3.58	0.00	0.56
Race_EthnBlack or African American	-0.30	0.33	-0.91	0.36	0.74
Race_EthnHispanic	-0.82	0.68	-1.20	0.23	0.44
Race_EthnNative American or Alaskan Native	0.88	0.29	3.08	0.00	2.42
Race_EthnOther Hispanic	0.35	0.41	0.85	0.39	1.42
Race_EthnWhite	0.40	0.26	1.53	0.12	1.49
Race_EthnWhite Hispanic	-0.37	0.29	-1.29	0.20	0.69
Never|Rarely	-3.45	0.36	-9.70	0.00	0.03
Rarely|Sometimes	-1.88	0.28	-6.78	0.00	0.15
Sometimes|Most of the time	-0.10	0.27	-0.38	0.70	0.90
Most of the time|Always	2.19	0.29	7.68	0.00	8.93

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.37 times more likely for for anglers that fish in both freshwater and saltwater compared to just freshwater anglers, 1.77 times more likely for anglers that were not prompted by COVID-19 to fish compared to anglers that were prompted by COVID-19 to fish, and 2.42 times more likely for Native American or Alaskan Native anglers.

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

# Order the factor
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"))
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, 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 + Income, design=svy_desn_Q13)
summary(Q13_mod2) # examine model results
VIF(Q13_mod2) # collinearity
Q13_mod3 <- svyolr(Q13_KeepOften ~ New.Old + SW_FW2 + Age + COVID, design=svy_desn_Q13)
summary(Q13_mod3) # examine model results
VIF(Q13_mod3) # collinearity

# Create odds ratio summary table
Q13_coef_table <- coef(summary(Q13_mod3))
Q13_pval <- pnorm(abs(Q13_coef_table[, "t value"]),lower.tail = FALSE)* 2
Q13_AOR <- exp(coef(Q13_mod3))
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
Q13_prop

# Reformat model data for figure
Q13_df1 <- model.frame(Q13_mod3) %>% 
  group_by(Q13_KeepOften,SW_FW2,New.Old,Age,COVID) %>%
  count() %>%
  na.omit()
Q13_freq = prop.table(Q13_df1$n)
Q13_df1 <- data.frame(cbind(Q13_df1,freq=Q13_freq))

# 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

# Histogram of responses
ggplot(Q13_df1,aes(x=Q13_KeepOften,y=freq,fill = New.Old)) + 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 = "Q13. How often do you keep fish?") + facet_wrap(~COVID, ncol = 2)

ggplot(Q13_df1,aes(x=Q13_KeepOften,y=freq,fill = Age)) + 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 = "Q13. How often do you keep fish?") + facet_wrap(~SW_FW2, ncol = 3)

# 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.40	0.12	-3.39	0.00	0.67
SW_FW2SW	1.22	0.15	8.38	0.00	3.40
SW_FW2SW.FW	0.78	0.14	5.80	0.00	2.19
Age25-34	0.19	0.23	0.81	0.42	1.21
Age35-44	0.35	0.22	1.55	0.12	1.42
Age45-54	0.18	0.22	0.81	0.42	1.20
Age55-64	0.57	0.23	2.49	0.01	1.77
Age65 or older	1.01	0.26	3.95	0.00	2.74
COVIDYes	-0.36	0.14	-2.61	0.01	0.70
Never|Rarely	-0.67	0.25	-2.69	0.01	0.51
Rarely|Sometimes	0.50	0.25	1.98	0.05	1.65
Sometimes|Most of the time	2.08	0.26	7.94	0.00	8.00
Most of the time|Always	3.72	0.28	13.37	0.00	41.27

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.49 times more likely for for anglers that had licenses prior to the pandemic (retained anglers), 1.43 times more likely for anglers that weren’t prompted by COVID-19 to start fishing, 3.40 times more likely for saltwater anglers compared to freshwater anglers, 2.19 times more likely for anglers that fish in both freshwater and saltwater compared to freshwater anglers, and 1.77-2.74 times more likely for anglers 55+.

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,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,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

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

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

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

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

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

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

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

# 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	14.078	0.000
Age	25	4.910	0.000
Race/Ethnicity	30	1.929	0.002
New/Retained	5	3.296	0.006
SW/FW	10	3.680	0.000
Gender	5	5.738	0.000
Children	10	10.450	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
Female	285 13.7 %	125 6 %	27 1.3 %	102 4.9 %	32 1.5 %	128 6.1 %	699 33.5 %
Male	433 20.7 %	367 17.6 %	114 5.5 %	174 8.3 %	60 2.9 %	240 11.5 %	1388 66.5 %
Total	718 34.4 %	492 23.6 %	141 6.8 %	276 13.2 %	92 4.4 %	368 17.6 %	2087 100 %

Q14_Gender_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Gender + Agreement, svy_desn_Q14, statistic = "adjWald")
## F = 5.7381, ndf = 5, ddf = 2086, p-value = 2.886e-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	89 4.3 %	144 6.9 %	41 2 %	63 3 %	5 0.2 %	101 4.8 %	443 21.2 %
25-34	146 7 %	131 6.3 %	34 1.6 %	98 4.7 %	33 1.6 %	79 3.8 %	521 25 %
35-44	175 8.4 %	101 4.8 %	36 1.7 %	56 2.7 %	37 1.8 %	89 4.3 %	494 23.7 %
45-54	153 7.3 %	68 3.3 %	18 0.9 %	33 1.6 %	11 0.5 %	52 2.5 %	335 16.1 %
55-64	110 5.3 %	33 1.6 %	9 0.4 %	19 0.9 %	4 0.2 %	28 1.3 %	203 9.7 %
65 or older	48 2.3 %	15 0.7 %	3 0.1 %	7 0.3 %	3 0.1 %	15 0.7 %	91 4.2 %
Total	721 34.5 %	492 23.6 %	141 6.8 %	276 13.2 %	93 4.5 %	364 17.4 %	2087 100 %

Q14_Age_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Age + Agreement, svy_desn_Q14, statistic = "adjWald")
## F = 4.9101, ndf = 25, ddf = 2066, p-value = 2.302e-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
Asian or Pacific Islander	19 1 %	24 1.2 %	9 0.5 %	22 1.1 %	2 0.1 %	14 0.7 %	90 4.6 %
Black or African American	20 1 %	27 1.4 %	13 0.7 %	17 0.9 %	7 0.4 %	23 1.2 %	107 5.6 %
Hispanic	14 0.7 %	15 0.8 %	3 0.2 %	12 0.6 %	5 0.3 %	10 0.5 %	59 3.1 %
Native American or Alaskan Native	22 1.1 %	12 0.6 %	7 0.4 %	12 0.6 %	3 0.2 %	9 0.5 %	65 3.4 %
Other Hispanic	33 1.7 %	44 2.2 %	6 0.3 %	16 0.8 %	4 0.2 %	9 0.5 %	112 5.7 %
White	453 23.1 %	239 12.2 %	71 3.6 %	127 6.5 %	58 3 %	197 10.1 %	1145 58.5 %
White Hispanic	112 5.7 %	91 4.6 %	31 1.6 %	67 3.4 %	8 0.4 %	73 3.7 %	382 19.4 %
Total	673 34.3 %	452 23.1 %	140 7.1 %	273 13.9 %	87 4.4 %	335 17.1 %	1960 100 %

Q14_Race_Ethn_chisq

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

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	309 13.8 %	230 10.3 %	56 2.5 %	137 6.1 %	44 2 %	186 8.3 %	962 43 %
SW	215 9.6 %	68 3 %	22 1 %	51 2.3 %	15 0.7 %	59 2.6 %	430 19.2 %
SW.FW	261 11.7 %	214 9.6 %	68 3 %	111 5 %	42 1.9 %	152 6.8 %	848 38 %
Total	785 35 %	512 22.9 %	146 6.5 %	299 13.3 %	101 4.5 %	397 17.7 %	2240 100 %

Q14_SW_FW2_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~SW_FW2 + Agreement, svy_desn_Q14, statistic = "adjWald")
## F = 3.6805, ndf = 10, ddf = 2081, p-value = 6.761e-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.	684 30.4 %	390 17.3 %	114 5.1 %	180 8 %	64 2.8 %	282 12.5 %	1714 76.1 %
Yes	103 4.6 %	128 5.7 %	32 1.4 %	120 5.3 %	37 1.6 %	116 5.2 %	536 23.8 %
Total	787 35 %	518 23 %	146 6.5 %	300 13.3 %	101 4.5 %	398 17.7 %	2250 100 %

Q14_COVID_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~COVID + Agreement, svy_desn_Q14, statistic = "adjWald")
## F = 14.078, ndf = 5, ddf = 2086, p-value = 1.447e-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	270 16 %	97 5.7 %	31 1.8 %	45 2.7 %	17 1 %	92 5.4 %	552 32.6 %
No children	221 13.1 %	233 13.8 %	74 4.4 %	152 9 %	5 0.3 %	151 8.9 %	836 49.5 %
Under 17	81 4.8 %	81 4.8 %	19 1.1 %	32 1.9 %	31 1.8 %	57 3.4 %	301 17.8 %
Total	572 33.9 %	411 24.3 %	124 7.3 %	229 13.6 %	53 3.1 %	300 17.8 %	1689 100 %

Q14_Children_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Children + Agreement, svy_desn_Q14, statistic = "adjWald")
## F = 10.45, 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$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")) 

# 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, family

# Find associations
findAssocs(dtm, terms = "time", corlimit = 0.1)
# top 5 words associated with time
# family quality   spend   husband   son 
#   0.37    0.21    0.17   0.15      0.15

findAssocs(dtm, terms = "fish", corlimit = 0.1)
# top 5 words associated with fish
# catching    size          big    efficient     trophy
# 0.67         0.16         0.13         0.12        0.12
      

findAssocs(dtm, terms = "family", corlimit = 0.1)    
# top 5 words associated with family
# time    friends   memories  longer   walk
# 0.37       0.32       0.15    0.13   0.13


# COVID anglers
#Create a vector containing only the text
text_COVID <- rev2_wts %>%
  filter(COVID == "Yes") 
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")) 

# 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, 3) # fish, time, family

# Find associations
findAssocs(dtm_COVID, terms = "time", corlimit = 0.2)
# top 5 words associated with time
# kids     family      loved ones    quality      covid 
# 0.31       0.30       0.28          0.28        0.27 

findAssocs(dtm_COVID, terms = "fish", corlimit = 0.2)
# words associated with fish
# catching     many          big        allot        banks 
# 0.64         0.27         0.27         0.24         0.24 

findAssocs(dtm_COVID, terms = "family", corlimit = 0.1)    
# words associated with family
# friends   time    outdoor    quality   covid 
# 0.33       0.30   0.19      0.16        0.13 

# 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 largely the same, only differences are the words "stress" and "people" in the COVID group and "quiet" and "husband" in the non-COVID group

findAssocs(dtm_COVID, terms = "stress", corlimit = 0.1)    
# free  thinking   reduce   relief    anything  release    relax 
# 0.43  0.37       0.37      0.3        0.26     0.21     0.11 

findAssocs(dtm_COVID, terms = "people", corlimit = 0.1)    
# lot      without        allot        banks         bass 
# 0.53         0.50         0.37         0.37         0.37 

findAssocs(dtm, terms = "quiet", corlimit = 0.1)    
# peace         shady          connect       come        husband  family 
# 0.40          0.15          0.15          0.15          0.15 

findAssocs(dtm, terms = "husband", corlimit = 0.1)    
# wade      trash       pick      adult      
# 0.16       0.16       0.16       0.16 

# Comparison of top 20 words
tab_df(df_comparison)

COVID.word	COVID.freq	Not_COVID.word	Not_COVID.freq
fish	64	fish	280
time	51	time	264
family	49	family	252
catching	37	catching	211
outdoors	26	outdoors	130
nature	21	relaxing	124
relaxing	18	nature	107
fun	18	water	95
water	17	friends	83
weather	15	fun	79
good	14	weather	75
outside	13	good	64
friends	12	outside	56
kids	10	kids	53
relax	9	relax	46
relaxation	9	quiet	44
peace	7	peace	44
stress	7	relaxation	42
people	7	husband	37
away	6	away	37

# Generate wordcloud
set.seed(122) # for reproducibility 
wordcloud(words = df$word, freq = df$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 anglers that were prompted or not prompted by COVID-19 to start fishing (COVID-prompted or non-COVID-prompted). The top 5 most frequent words were fish, time, family, catching, and outdoors for both groups. There was an overlap of 18 of the top 20 words associated with each group. The words that did not overlap in the top 20 word list were “stress” and “people” for the COVID-prompted group and “quiet” and “husband” for the non-COVID-prompted angling group. The word “people” was associated with the words: lot, without, and banks. The word “stress” was associated with the words: free, thinking, and reduce. The word “husband” was associated with the words: trash, wade, and pick. The word “quiet” was associated with the words: peace, shady, and connect.

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 <- 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"))
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, design=svy_desn_Q16)
summary(Q16_mod1) # examine model results
VIF(Q16_mod1) # collinearity
Q16_mod2 <- svyolr(Q16_TravelTime ~ Income + SW_FW2 + Age + 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

# Reformat model data for figure
Q16_df1 <- model.frame(Q16_mod2) %>% 
  group_by(Q16_TravelTime,Income,SW_FW2,Age,Race_Ethn) %>%
  count() %>%
  na.omit()
Q16_freq = prop.table(Q16_df1$n)
Q16_df1 <- data.frame(cbind(Q16_df1,freq=Q16_freq))

Q16_AOR_I20K <- 1/0.47 # 2.128
Q16_AOR_RW <- 1/0.54 # 1.85
Q16_AOR_RB <- 1/0.58 # 1.72

# 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

# Histogram of responses
ggplot(Q16_df1,aes(x=Q16_TravelTime,y=freq,fill = Age)) + 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 = "Q16. How long did it take you to get to your fishing spot?") + facet_wrap(~SW_FW2, ncol = 3)

ggplot(Q16_df1,aes(x=Q16_TravelTime,y=freq,fill = Income)) + 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 = "Q16. How long did it take you to get to your fishing spot?") + facet_wrap(~SW_FW2, ncol = 3)

ggplot(Q16_df1,aes(x=Q16_TravelTime,y=freq,fill = Race_Ethn)) + 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 = "Q16. How long did it take you to get to your fishing spot?") + facet_wrap(~SW_FW2, ncol = 3)

# 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
IncomeBetween $20,000 and $39,999	0.14	0.41	0.35	0.73	1.16
IncomeBetween $40,000 and $59,999	0.70	0.37	1.92	0.06	2.01
IncomeBetween $60,000 and $79,999	0.49	0.37	1.33	0.18	1.63
IncomeBetween $80,000 and $99,999	0.69	0.37	1.85	0.06	1.99
IncomeBetween $100,000 and $149,999	0.75	0.35	2.13	0.03	2.11
IncomeOver $150,000	0.96	0.36	2.64	0.01	2.61
SW_FW2SW	0.79	0.17	4.75	0.00	2.20
SW_FW2SW.FW	0.65	0.15	4.27	0.00	1.92
Age25-34	0.36	0.29	1.22	0.22	1.43
Age35-44	0.33	0.28	1.19	0.23	1.40
Age45-54	0.42	0.28	1.52	0.13	1.53
Age55-64	0.53	0.29	1.81	0.07	1.69
Age65 or older	0.70	0.31	2.24	0.03	2.01
Race_EthnBlack or African American	-0.61	0.31	-1.97	0.05	0.54
Race_EthnHispanic	0.05	0.45	0.11	0.91	1.05
Race_EthnNative American or Alaskan Native	-0.84	0.56	-1.49	0.14	0.43
Race_EthnOther Hispanic	-0.37	0.38	-0.99	0.32	0.69
Race_EthnWhite	-0.55	0.27	-2.03	0.04	0.58
Race_EthnWhite Hispanic	-0.52	0.30	-1.74	0.08	0.59
Less than 30 minutes|Between 30 and 60 minutes	0.25	0.45	0.56	0.58	1.28
Between 30 and 60 minutes|Between 1 and 2 hours	1.52	0.45	3.38	0.00	4.57
Between 1 and 2 hours|More than 2 hours	2.68	0.45	5.90	0.00	14.53

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.11-2.61 times more likely for anglers with an income >$100K compared to anglers with the lowest incomes (<$20K), 2.20 times more likely for saltwater anglers compared to freshwater anglers, 1.92 times more likely for anglers that fish in both freshwater and saltwater compared to freshwater anglers, 1.69 - 2.01 times more likely for anglers 55+ compared to younger anglers, and 1.72 - 1.85 times more likely for Asian or Pacific Islander anglers compared to white and Black or African American anglers

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

# Order the factor
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"))
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, 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, 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, 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

# Reformat model data for figure
Q17_df1 <- model.frame(Q17_mod3) %>% 
  group_by(Q17_MaxTravel,Income,COVID,SW_FW2) %>%
  count() %>%
  na.omit()
Q17_freq = prop.table(Q17_df1$n)
Q17_df1 <- data.frame(cbind(Q17_df1,freq=Q17_freq))

Q17_AOR_I20K <- 1/0.39 # 2.56
Q16_AOR_Covid <- 1/0.66 # 1.52
Q16_AOR_RB <- 1/0.58 # 1.72

# 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

# Histogram of responses
ggplot(Q17_df1,aes(x=Q17_MaxTravel,y=freq,fill = Income)) + 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 = "Q17. How far are you willing to travel to go fishing?") + facet_wrap(~SW_FW2, ncol = 3)

ggplot(Q17_df1,aes(x=Q17_MaxTravel,y=freq,fill = COVID)) + 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 = "Q17. How far are you willing to travel to go fishing?")

# 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
IncomeBetween $20,000 and $39,999	0.94	0.49	1.91	0.06	2.57
IncomeBetween $40,000 and $59,999	1.13	0.47	2.38	0.02	3.08
IncomeBetween $60,000 and $79,999	1.09	0.48	2.25	0.03	2.97
IncomeBetween $80,000 and $99,999	0.79	0.46	1.73	0.08	2.20
IncomeBetween $100,000 and $149,999	0.95	0.45	2.12	0.03	2.60
IncomeOver $150,000	1.44	0.46	3.11	0.00	4.24
COVIDYes	-0.38	0.15	-2.43	0.01	0.69
SW_FW2SW	0.07	0.17	0.41	0.68	1.07
SW_FW2SW.FW	0.65	0.15	4.34	0.00	1.92
New.OldNew	-0.34	0.13	-2.56	0.01	0.71
Less than 30 minutes|Between 30 and 60 minutes	-2.62	0.49	-5.40	0.00	0.07
Between 30 and 60 minutes|Between 1 and 2 hours	-0.33	0.45	-0.72	0.47	0.72
Between 1 and 2 hours|Between 2 and 4 hours	1.00	0.46	2.19	0.03	2.73
Between 2 and 4 hours|More than 4 hours	2.11	0.47	4.48	0.00	8.27

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 desired fishing spots within a 30 minute drive.

The odds of being willing to travel further to fish were: 2.20-4.24 times more likely for anglers with an income >$40K compared to anglers with the lowest incomes (<$20K), 1.41 times more likely for anglers that had licenses prior to the pandemic (retained anglers), 1.45 times more likely for anglers that were not prompted by COVID-19 to start fishing, and 1.92 times more likely for anglers that fish in both freshwater and saltwater compared to freshwater anglers.

Q18. 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 survey answers
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"))
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)
  
# Remove rows with missing independent variable data
Q18_df1 <- rev2_wts %>%
  filter_at(vars(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2),all_vars(!is.na(.))) 
# Select independent variables
Q18_indep1 <- Q18_df1 %>%
  dplyr::select(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2)
# Select dependent variables
Q18_dep1 <- Q18_df1 %>%
  dplyr::select(Motorboat_access,Peaceful_areas,Portable_restrooms,Paddling_trail,Fish_stockings,Parking_area,ADA_accessible,Baby_changing_stations,Shoreline_access)
# Transform dependent data
Q18_dep2 <- decostand(Q18_dep1, method = 'normalize') # normalize
# MCA, keep 2 axes, on independent data
Q18_indep2 <- dudi.mix(Q18_indep1, scannf=F, nf=2)$tab # MCA, keep 2 axes

# Run RDA
Q18_rda_mod1 <- rda(Q18_dep2 ~ ., data = Q18_indep2,scale = TRUE) 
envfit(Q18_rda_mod1,Q18_indep2) # select variables at alpha = 0.001 for next model 
Q18_indep3 <- Q18_indep2 %>%
  dplyr::select(Incom.Between..100.000.and..149.999,Incom.Over..150.000,Age.25.34, Age.35.44, Age.45.54, Age.65.or.older,Child.No.children,Child.Under.17,New.O.New,New.O.Old, Race_.Black.or.African.American,Race_.White, Race_.White.Hispanic,SW_FW.SW,SW_FW.FW,COVID.No..I.was.already.planning.to.go.fishing.,COVID.Yes)
Q18_rda_mod2 <- rda(Q18_dep2 ~ ., data = Q18_indep3)
envfit(Q18_rda_mod2,Q18_indep3) # select variables at alpha = 0.001 for next model
Q18_indep4 <- Q18_indep2 %>%
  dplyr::select(Incom.Over..150.000,Age.25.34,Age.65.or.older,Child.No.children,Child.Under.17,New.O.New,New.O.Old, Race_.Black.or.African.American,Race_.White, Race_.White.Hispanic, SW_FW.SW,SW_FW.FW,COVID.No..I.was.already.planning.to.go.fishing.,COVID.Yes)
Q18_rda_mod3 <- rda(Q18_dep2 ~ ., data = Q18_indep4)
envfit(Q18_rda_mod3,Q18_indep4) # select variables at alpha = 0.001 for next model
Q18_rda_mod3$call
RsquareAdj(Q18_rda_mod3) # 0.04 
sqrt(vif.cca(Q18_rda_mod3))

# Extract scores for visualization
# Basic plot for visualization
vare_spp_sco <- scores(Q18_rda_mod3, display = "species",scaling = "species")
vare_sam_sco <- scores(Q18_rda_mod3, display = "sites",scaling = "sites")
vare_env_sco <- scores(Q18_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 65 over","No children", "Children Under 17","New Angler","Black or African American","White","White Hispanic","Saltwater","Freshwater","COVID-No"),ccatype = "bp")
rescaled <- vare_env_tbl %>% 
            dplyr::select(RDA1, RDA2) %>%
            as.matrix() *4.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_summary1 <- rev1_wts %>%  
  mutate_at(vars(Q18A:Q18I), list(~recode_factor(.,"Strongly Disagree" = "1","Disagree" = "2","Neither Agree nor Disagree" = "3","Agree" = "4","Strongly Agree"= "5"))) %>%
  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) %>%
  dplyr::select(Motorboat_access,Peaceful_areas,Portable_restrooms,Paddling_trail,Fish_stockings,Parking_area,ADA_accessible,Baby_changing_stations,Shoreline_access)

Q18_summary2 <- describe(Q18_summary1,na.rm = TRUE) 
Q18_summary2$Variables <- row.names(Q18_summary2) 
Q18_summary3 <- Q18_summary2[,c("Variables","mean","sd")]
Q18_summary3$Variables<-str_replace_all(Q18_summary3$Variables, c("_" = " ")) # replaces spaces in names with periods

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

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

# Descriptive stats table
tab_df(Q18_summary3)

Variables	mean	sd
Motorboat access	3.52	3.13
Peaceful areas	4.46	3.70
Portable restrooms	3.89	3.28
Paddling trail	3.19	2.99
Fish stockings	3.96	3.55
Parking area	4.27	3.65
ADA accessible	2.97	2.75
Baby changing stations	2.68	2.68
Shoreline access	4.36	3.80

# 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") 

# PERMANOVA
anova.cca(Q18_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_dep2 ~ Incom.Over..150.000 + Age.25.34 + Age.65.or.older + Child.No.children + Child.Under.17 + New.O.New + New.O.Old + Race_.Black.or.African.American + Race_.White + Race_.White.Hispanic + SW_FW.SW + SW_FW.FW + COVID.No..I.was.already.planning.to.go.fishing. + COVID.Yes, data = Q18_indep4)
##           Df Variance       F Pr(>F)    
## RDA1       1  0.01245 26.3930  0.001 ***
## RDA2       1  0.00516 10.9349  0.001 ***
## RDA3       1  0.00265  5.6118  0.114    
## RDA4       1  0.00252  5.3508  0.087 .  
## RDA5       1  0.00100  2.1188  0.947    
## RDA6       1  0.00047  0.9921  1.000    
## RDA7       1  0.00021  0.4432  1.000    
## RDA8       1  0.00015  0.3190  1.000    
## RDA9       1  0.00011  0.2251  1.000    
## Residual 949  0.44780                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_df(Q18_envfit2)

Variables	RDA.1	RDA.2	r2	p
Incom.Over..150.000	0.93	0.36	0.02	0.00
Age.25.34	-1.00	-0.02	0.03	0.00
Age.65.or.older	1.00	0.05	0.02	0.00
Child.No.children	-0.51	-0.86	0.02	0.00
Child.Under.17	-0.85	0.53	0.05	0.00
New.O.New	-0.80	-0.60	0.03	0.00
New.O.Old	0.80	0.60	0.03	0.00
Race_.Black.or.African.American	-0.61	-0.79	0.01	0.00
Race_.White	0.99	0.11	0.05	0.00
Race_.White.Hispanic	-0.99	0.16	0.02	0.00
SW_FW.SW	0.30	0.96	0.04	0.00
SW_FW.FW	-0.26	-0.97	0.02	0.00
COVID.No..I.was.already.planning.to.go.fishing.	0.72	0.70	0.02	0.00
COVID.Yes	-0.72	-0.70	0.02	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: parking areas, peaceful areas, portable restrooms, and shoreline access. Fish stockings, motorboat access, and paddling trails were less important, while having baby changing stations and ADA accessible facilities were least important overall.

Anglers between the ages of 25-34 and those that identified as White Hispanic favored sites with paddling trails. Anglers that recruited during the pandemic, Black or African American anglers, anglers without children, and freshwater anglers favored fish stockings and peaceful areas. Anglers age 65 and over, anglers that were not prompted by the pandemic to begin fishing, 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 were more likely to favor sites with baby changing stations and portable restrooms.

Cronbach’s alpha for Q18 is 0.9868816

Q19. 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:

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"))

# Clean Q19 survey answers
rev2_wts <- rev2_wts %>%  
  mutate_at(vars(Q19A:Q19H), list(~recode_factor(.,"Strongly Disagree" = "-2","Disagree" = "-1","Neither Agree nor Disagree" = "0","Agree" = "1","Strongly Agree"= "2"))) %>%
  mutate(across(Q19A:Q19H, ~as.numeric(as.character(.)))) %>%
  mutate_at(vars(Q19A:Q19H), list(~.*weights)) %>%
  mutate(Canoe_access = Q19A,
         Space = Q19B,
         Brick_restrooms = Q19C,
         Lighting = Q19D,
         Picnic_areas = Q19E,
         Playground = Q19F,
         Shade = Q19G,
         Swimming = Q19H)
  
# Remove rows with missing independent variable data
Q19_df1 <- rev2_wts %>%
  filter_at(vars(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2),all_vars(!is.na(.))) 
# Select independent variables
Q19_indep1 <- Q19_df1 %>%
  dplyr::select(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2)
# Select dependent variables
Q19_dep1 <- Q19_df1 %>%
  dplyr::select(Canoe_access,Space,Brick_restrooms,Lighting,Picnic_areas,Playground,Shade,Swimming)
# Transform dependent data
Q19_dep2 <- decostand(Q19_dep1, method = 'normalize')
# MCA, keep 2 axes, on independent data
Q19_indep2 <- dudi.mix(Q19_indep1, scannf=F, nf=2)$tab # MCA, keep 2 axes

# Run RDA
Q19_rda_mod1 <- rda(Q19_dep2 ~ ., data = Q19_indep2) 
envfit(Q19_rda_mod1,Q19_indep2) # select variables at alpha = 0.001 for next model 
Q19_indep3 <- Q19_indep2 %>%
  dplyr::select(Age.25.34,Age.65.or.older,Child.18.or.over,Child.No.children,Child.Under.17,New.O.New,New.O.Old, Race_.White, Race_.White.Hispanic,SW_FW.SW)
Q19_rda_mod2 <- rda(Q19_dep2 ~ ., data = Q19_indep3)
envfit(Q19_rda_mod2,Q19_indep3) # select variables at alpha = 0.001 for next model
Q19_rda_mod2$call
RsquareAdj(Q19_rda_mod2) 
sqrt(vif.cca(Q19_rda_mod2))

# Extract scores for visualization
# Basic plot for visualization
vare_spp_sco <- scores(Q19_rda_mod2, display = "species")
vare_sam_sco <- scores(Q19_rda_mod2, display = "sites")
vare_env_sco <- scores(Q19_rda_mod2, 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=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=c("Age 25-34","Age 65 over","Children over 18","No Children","New Angler","White","White Hispanic","Saltwater"),ccatype = "bp")
rescaled <- vare_env_tbl %>% 
            dplyr::select(RDA1, RDA2) %>%
            as.matrix() * 0.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
Q19_summary1 <- rev1_wts %>%  
  mutate_at(vars(Q19A:Q19H), list(~recode_factor(.,"Strongly Disagree" = "1","Disagree" = "2","Neither Agree nor Disagree" = "3","Agree" = "4","Strongly Agree"= "5"))) %>%
  mutate(across(Q19A:Q19H, ~as.numeric(as.character(.)))) %>%
  mutate_at(vars(Q19A:Q19H), list(~.*weights)) %>%
  mutate(Canoe_access = Q19A,
         Space = Q19B,
         Brick_restrooms = Q19C,
         Lighting = Q19D,
         Picnic_areas = Q19E,
         Playground = Q19F,
         Shade = Q19G,
         Swimming = Q19H) %>%
  dplyr::select(Canoe_access,Space,Brick_restrooms,Lighting,Picnic_areas,Playground,Shade,Swimming)

Q19_summary2 <- describe(Q19_summary1,na.rm = TRUE) 
Q19_summary2$Variables <- row.names(Q19_summary2) 
Q19_summary3 <- Q19_summary2[,c("Variables","mean","sd")]
Q19_summary3$Variables<-str_replace_all(Q19_summary3$Variables, c("_" = " ")) # replaces spaces in names with periods

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

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

# 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(Q19_rda_mod2, 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 = Q19_dep2 ~ Age.25.34 + Age.65.or.older + Child.18.or.over + Child.No.children + Child.Under.17 + New.O.New + New.O.Old + Race_.White + Race_.White.Hispanic + SW_FW.SW, data = Q19_indep3)
##           Df Variance       F Pr(>F)    
## RDA1       1  0.01291 22.3734  0.001 ***
## RDA2       1  0.00637 11.0451  0.001 ***
## RDA3       1  0.00183  3.1669  0.368    
## RDA4       1  0.00136  2.3558  0.546    
## RDA5       1  0.00112  1.9393  0.568    
## RDA6       1  0.00040  0.6848  0.990    
## RDA7       1  0.00003  0.0520  1.000    
## RDA8       1  0.00001  0.0260  1.000    
## Residual 950  0.54800                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_df(Q19_envfit2)

Variables	RDA.1	RDA.2	r2	p
Age.25.34	-0.45	0.89	0.02	0.00
Age.65.or.older	0.59	-0.81	0.01	0.01
Child.18.or.over	0.43	-0.90	0.03	0.00
Child.No.children	0.25	0.97	0.06	0.00
Child.Under.17	-0.81	-0.58	0.06	0.00
New.O.New	-0.27	0.96	0.04	0.00
New.O.Old	0.27	-0.96	0.04	0.00
Race_.White	0.87	0.49	0.03	0.00
Race_.White.Hispanic	-0.96	-0.29	0.02	0.00
SW_FW.SW	0.02	-1.00	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=2,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")

Between 60-95% of anglers agreed that all of the following features would increase their willingness to fish at a site: space, shade, lighting, picnic areas and brick restrooms. Anglers were less willing to fish at sites with playgrounds, swimming, and were approximately neutral about canoe accessibility.

Anglers that recruited during the pandemic and anglers between the ages of 25-34 were most willing to fish at sites with shade and canoe access, while anglers over 65, with children over 18, and saltwater anglers were less likely to report those feature as important. White anglers were more willing to fish at sites with more space. Finally, white hispanic anglers were more willing to fish at sites with playgrounds, swimming, picnic areas, lighting, and brick restrooms.

Cronbach’s alpha for Q19 is 0.9857154

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

# Order the factor
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"))
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, 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

# Reformat model data for figure
Q20_df1 <- model.frame(Q20_mod2) %>% 
  group_by(Q20_NewLicense,New.Old,COVID,SW_FW2) %>%
  count() %>%
  na.omit()
Q20_freq = prop.table(Q20_df1$n)
Q20_df1 <- data.frame(cbind(Q20_df1,freq=Q20_freq))

# 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

# Histogram of responses
ggplot(Q20_df1,aes(x=Q20_NewLicense,y=freq,fill = New.Old)) + 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 = "Q20. How likely are you to purchase another license?") + facet_wrap(~COVID, ncol = 2)

ggplot(Q20_df1,aes(x=Q20_NewLicense,y=freq,fill = SW_FW2)) + 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 = "Q20. How likely are you to purchase another license?") 

ggplot(Q20_df1,aes(x=Q20_NewLicense,y=freq,fill = New.Old)) + 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 = "Q20. How likely are you to purchase another license?") 

ggplot(Q20_df1,aes(x=Q20_NewLicense,y=freq,fill = COVID)) + 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 = "Q20. How likely are you to purchase another license?") 

# 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.19	0.31	-3.86	0.00	0.31
COVIDYes	-1.02	0.25	-4.06	0.00	0.36
SW_FW2SW	-0.43	0.29	-1.49	0.14	0.65
SW_FW2SW.FW	0.63	0.29	2.19	0.03	1.88
Not likely|Somewhat likely	-6.12	0.50	-12.13	0.00	0.00
Somewhat likely|Very likely	-3.77	0.34	-10.98	0.00	0.02

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

The odds of being likely to purchase another fishing license were: 3.27 times more likely for anglers that purchased a license prior to the pandemic, 2.28 times more likely for anglers that reported that were not prompted by the pandemic to fish, and 1.88 times more likely for anglers that fish in both freshwater and saltwater compared to anglers that only fish in freshwater.

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

# Order the factor
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"))
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, 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

# Reformat model data for figure
Q21_df1 <- model.frame(Q21_mod2) %>% 
  group_by(Q21_AutoRenew,Income,Race_Ethn) %>%
  count() %>%
  na.omit()
Q21_freq = prop.table(Q21_df1$n)
Q21_df1 <- data.frame(cbind(Q21_df1,freq=Q21_freq))

# 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

# Histogram of responses
ggplot(Q21_df1,aes(x=Q21_AutoRenew,y=freq,fill = Income)) + 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 = "Q21. How likely are you to enroll in auto-renew license program?") 

ggplot(Q21_df1,aes(x=Q21_AutoRenew,y=freq,fill = Race_Ethn)) + 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 = "Q21. How likely are you to enroll in auto-renew license program?") 

# 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
IncomeBetween $20,000 and $39,999	-0.06	0.43	-0.14	0.89	0.94
IncomeBetween $40,000 and $59,999	0.32	0.41	0.78	0.44	1.37
IncomeBetween $60,000 and $79,999	0.70	0.39	1.77	0.08	2.01
IncomeBetween $80,000 and $99,999	0.63	0.41	1.55	0.12	1.88
IncomeBetween $100,000 and $149,999	0.69	0.39	1.77	0.08	2.00
IncomeOver $150,000	1.16	0.40	2.92	0.00	3.19
Race_EthnBlack or African American	-0.69	0.40	-1.73	0.08	0.50
Race_EthnHispanic	0.01	0.43	0.03	0.98	1.01
Race_EthnNative American or Alaskan Native	-0.51	0.85	-0.60	0.55	0.60
Race_EthnOther Hispanic	0.92	0.41	2.23	0.03	2.52
Race_EthnWhite	-0.03	0.27	-0.12	0.90	0.97
Race_EthnWhite Hispanic	0.14	0.31	0.44	0.66	1.15
Extremely unlikely|Unlikely	-5.02	0.55	-9.08	0.00	0.01
Unlikely|Not sure	-2.35	0.43	-5.40	0.00	0.10
Not sure|Likely	-0.48	0.43	-1.10	0.27	0.62
Likely|Extremely likely	0.84	0.44	1.93	0.05	2.32

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: 3.19 times more likely for anglers with an income >$150K compared to anglers with an income <$20K, 1.96 times more likely for anglers that identified as an Other Hispanic race compared to Asian and Pacific Islander anglers.

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,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,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

# 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)

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

# 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.074	0.000
Children	16	2.458	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	30 0.5 %	249 3.8 %	171 2.6 %	254 3.8 %	145 2.2 %	39 0.6 %	145 2.2 %	100 1.5 %	224 3.4 %	1357 20.6 %
25-34	25 0.4 %	255 3.9 %	209 3.2 %	290 4.4 %	185 2.8 %	30 0.5 %	188 2.8 %	131 2 %	270 4.1 %	1583 24.1 %
35-44	49 0.7 %	239 3.6 %	190 2.9 %	293 4.4 %	142 2.2 %	55 0.8 %	193 2.9 %	121 1.8 %	266 4 %	1548 23.3 %
45-54	44 0.7 %	160 2.4 %	126 1.9 %	216 3.3 %	91 1.4 %	28 0.4 %	144 2.2 %	97 1.5 %	201 3 %	1107 16.8 %
55-64	29 0.4 %	95 1.4 %	72 1.1 %	134 2 %	52 0.8 %	22 0.3 %	104 1.6 %	58 0.9 %	116 1.8 %	682 10.3 %
65 or older	19 0.3 %	42 0.6 %	33 0.5 %	71 1.1 %	18 0.3 %	8 0.1 %	46 0.7 %	32 0.5 %	53 0.8 %	322 4.9 %
Total	196 3 %	1040 15.8 %	801 12.1 %	1258 19.1 %	633 9.6 %	182 2.8 %	820 12.4 %	539 8.2 %	1130 17.1 %	6599 100 %

Q24_Age_chisq

## 
##  Design-based Wald test of association
## 
## data:  svychisq(~Age + Agreement, svy_desn_Q24, statistic = "adjWald")
## F = 2.0735, ndf = 40, ddf = 6369, p-value = 8.487e-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	30 0.5 %	249 3.8 %	171 2.6 %	254 3.8 %	145 2.2 %	39 0.6 %	145 2.2 %	100 1.5 %	224 3.4 %	1357 20.6 %
25-34	25 0.4 %	255 3.9 %	209 3.2 %	290 4.4 %	185 2.8 %	30 0.5 %	188 2.8 %	131 2 %	270 4.1 %	1583 24.1 %
35-44	49 0.7 %	239 3.6 %	190 2.9 %	293 4.4 %	142 2.2 %	55 0.8 %	193 2.9 %	121 1.8 %	266 4 %	1548 23.3 %
45-54	44 0.7 %	160 2.4 %	126 1.9 %	216 3.3 %	91 1.4 %	28 0.4 %	144 2.2 %	97 1.5 %	201 3 %	1107 16.8 %
55-64	29 0.4 %	95 1.4 %	72 1.1 %	134 2 %	52 0.8 %	22 0.3 %	104 1.6 %	58 0.9 %	116 1.8 %	682 10.3 %
65 or older	19 0.3 %	42 0.6 %	33 0.5 %	71 1.1 %	18 0.3 %	8 0.1 %	46 0.7 %	32 0.5 %	53 0.8 %	322 4.9 %
Total	196 3 %	1040 15.8 %	801 12.1 %	1258 19.1 %	633 9.6 %	182 2.8 %	820 12.4 %	539 8.2 %	1130 17.1 %	6599 100 %

Q24_Children_chisq

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

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.

Tables for publication

tab_dfs(list(COVID,New.Old,Children,SW_FW2,Race_Ethn,Income),file = "Demographics.doc") # Demographics summary tables

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_summary3,Q19_summary3),file = "Likert_summary_tables.doc") # Summary responses to likert scale questions

tab_dfs(list(Q5_envfit2,Q11_envfit2,Q18_envfit2,Q19_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_mod4,Q12_mod3,Q13_mod3,Q16_mod2,Q17_mod2,Q20_mod2,Q21_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?"),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.cell.prc = TRUE,show.summary = FALSE,file = "Q7_New.Old.doc")
tab_xtab(rev2_wts$Gender,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE,file = "Q7_Gender.doc")
tab_xtab(rev2_wts$Age,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE,file = "Q7_Age.doc")
tab_xtab(rev2_wts$Children,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE,file = "Q7_Children.doc")
tab_xtab(rev2_wts$Race_Ethn,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE,file = "Q7_Race_Ethn.doc")
tab_xtab(rev2_wts$Income,rev2_wts$SW_FW2,weight.by = rev2_wts$weights,show.cell.prc = TRUE,show.summary = FALSE,file = "Q7_Income.doc")


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

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

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

tab_dfs(list(COVID,New.Old,Children,SW_FW2,Race_Ethn,Income)) # Demographics summary tables

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

 

New.Old	n	freq
New	1517	0.68
Old	721	0.32

 

Children	n	freq
18 or over	852	0.56
No children	501	0.33
Under 17	182	0.12

 

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

 

Race_Ethn	n	freq
Asian or Pacific Islander	64	0.04
Black or African American	84	0.05
Hispanic	41	0.02
Native American or Alaskan Native	40	0.02
Other Hispanic	74	0.04
White	1220	0.69
White Hispanic	242	0.14

 

Income	n	freq
Less than $20,000	70	0.05
Between $20,000 and $39,999	142	0.10
Between $40,000 and $59,999	186	0.13
Between $60,000 and $79,999	226	0.15
Between $80,000 and $99,999	211	0.14
Between $100,000 and $149,999	332	0.23
Over $150,000	298	0.20

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_summary3,Q19_summary3)) # Summary responses to likert scale questions

Variables	mean	sd
Food	1.82	1.91
Sport	2.46	2.43
Outdoors	3.42	2.82
Nature	3.41	2.85
Friends	3.22	2.68
Children1	2.52	2.08
Stress	3.16	2.77
Relax	3.40	2.86
Get away	2.55	2.56

 

Variables	mean	sd
Enjoyed fishing	3.91	3.12
Too expensive	2.13	1.96
Fishing rules confusing	1.98	1.64
Enjoyed surroundings	3.82	3.02

 

Variables	mean	sd
Motorboat access	3.52	3.13
Peaceful areas	4.46	3.70
Portable restrooms	3.89	3.28
Paddling trail	3.19	2.99
Fish stockings	3.96	3.55
Parking area	4.27	3.65
ADA accessible	2.97	2.75
Baby changing stations	2.68	2.68
Shoreline access	4.36	3.80

 

Variables	mean	sd
Canoe access	3.62	3.55
Space	4.41	3.75
Brick restrooms	3.69	3.22
Lighting	3.95	3.58
Picnic areas	3.67	3.26
Playground	3.02	2.82
Shade	4.08	3.63
Swimming	3.18	3.09

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.31	0.00
SW/FW vs Age	10	4.17	0.00
SW/FW vs Gender	2	4.24	0.01
SW/FW vs Children	4	5.42	0.00
SW/FW vs Race/Ethnicity	12	3.40	0.00

 

Variable	df	F	p
COVID-Prompted	5	2.31	0.04
Age	25	22.17	0.00
Race/Ethnicity	30	1.89	0.00
Income	30	1.60	0.02
SW/FW	10	3.48	0.00
Gender	5	7.34	0.00
Children	10	39.01	0.00

 

Variable	df	F	p
COVID-Prompted	5	14.08	0.00
Age	25	4.91	0.00
Race/Ethnicity	30	1.93	0.00
New/Retained	5	3.30	0.01
SW/FW	10	3.68	0.00
Gender	5	5.74	0.00
Children	10	10.45	0.00

 

Variable	df	F	p
Age	40	2.07	0.00
Children	16	2.46	0.00

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

Variables	RDA 1	RDA 2	r2	p
Gende.Female	-0.43	0.90	0.09	0.00
Gende.Male	0.43	-0.90	0.09	0.00
Age.18.24	-0.98	0.22	0.02	0.00
Age.25.34	-0.98	0.22	0.05	0.00
Age.65.or.older	0.81	-0.59	0.02	0.00
Child.18.or.over	1.00	-0.10	0.13	0.00
Child.No.children	-1.00	0.08	0.26	0.00
Child.Under.17	1.00	-0.03	0.03	0.00
COVID.No..I.was.already.planning.to.go.fishing.	0.89	-0.45	0.02	0.00
COVID.Yes	-0.89	0.45	0.02	0.00
New.O.New	-0.65	0.76	0.03	0.00
New.O.Old	0.65	-0.76	0.03	0.00

 

Variables	RDA 1	RDA 2	r2	p
Age.18.24	1.00	-0.02	0.58	0.00
Age.25.34	1.00	-0.03	0.07	0.00
Age.55.64	-1.00	0.02	0.10	0.00
Age.65.or.older	-1.00	0.02	0.05	0.00
Child.18.or.over	-1.00	0.02	0.20	0.00
Child.No.children	1.00	-0.02	0.12	0.00
New.O.New	1.00	0.02	0.01	0.00
New.O.Old	-1.00	-0.02	0.01	0.00

 

Variables	RDA 1	RDA 2	r2	p
Incom.Over..150.000	0.93	0.36	0.02	0.00
Age.25.34	-1.00	-0.02	0.03	0.00
Age.65.or.older	1.00	0.05	0.02	0.00
Child.No.children	-0.51	-0.86	0.02	0.00
Child.Under.17	-0.85	0.53	0.05	0.00
New.O.New	-0.80	-0.60	0.03	0.00
New.O.Old	0.80	0.60	0.03	0.00
Race_.Black.or.African.American	-0.61	-0.79	0.01	0.00
Race_.White	0.99	0.11	0.05	0.00
Race_.White.Hispanic	-0.99	0.16	0.02	0.00
SW_FW.SW	0.30	0.96	0.04	0.00
SW_FW.FW	-0.26	-0.97	0.02	0.00
COVID.No..I.was.already.planning.to.go.fishing.	0.72	0.70	0.02	0.00
COVID.Yes	-0.72	-0.70	0.02	0.00

 

Variables	RDA 1	RDA 2	r2	p
Age.25.34	-0.45	0.89	0.02	0.00
Age.65.or.older	0.59	-0.81	0.01	0.01
Child.18.or.over	0.43	-0.90	0.03	0.00
Child.No.children	0.25	0.97	0.06	0.00
Child.Under.17	-0.81	-0.58	0.06	0.00
New.O.New	-0.27	0.96	0.04	0.00
New.O.Old	0.27	-0.96	0.04	0.00
Race_.White	0.87	0.49	0.03	0.00
Race_.White.Hispanic	-0.96	-0.29	0.02	0.00
SW_FW.SW	0.02	-1.00	0.03	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.11	0.05 – 0.23	<0.001	25.35	9.61 – 66.87	<0.001	10.80	4.86 – 24.00	<0.001	3.14	1.02 – 9.63	0.046	0.36	0.13 – 0.97	0.044
COVID:Yes													0.69	0.49 – 0.98	0.036
New/Retained:New	1.53	1.01 – 2.32	0.047										0.69	0.50 – 0.94	0.020
Gender:Male	0.55	0.39 – 0.77	0.001	0.52	0.40 – 0.68	<0.001				1.53	1.12 – 2.07	0.007
Children:No Children	2.12	1.47 – 3.07	<0.001
Children:Under 17	1.41	0.83 – 2.38	0.200
SW/FW:SW				0.51	0.36 – 0.73	<0.001							1.94	1.34 – 2.82	0.001
SW/FW:Both				0.71	0.52 – 0.99	0.043							2.18	1.55 – 3.08	<0.001
Age:25-34 yrs	2.06	1.07 – 3.96	0.031	0.58	0.30 – 1.16	0.123	1.79	0.67 – 4.83	0.248	0.38	0.16 – 0.90	0.028
Age:35-44 yrs	1.75	0.91 – 3.36	0.091	0.51	0.27 – 0.96	0.038	1.09	0.45 – 2.62	0.850	0.37	0.16 – 0.86	0.022
Age:45-54 yrs	1.60	0.81 – 3.15	0.176	0.36	0.19 – 0.67	0.001	0.70	0.30 – 1.63	0.409	0.41	0.17 – 0.97	0.043
Age:55-64 yrs	1.64	0.84 – 3.19	0.149	0.41	0.22 – 0.79	0.007	0.67	0.28 – 1.56	0.347	0.39	0.16 – 0.92	0.031
Age:65 yrs or older	0.96	0.44 – 2.09	0.911	0.38	0.19 – 0.75	0.005	0.27	0.12 – 0.64	0.003	0.21	0.08 – 0.51	0.001
Race/Ethnicity:Black or African American				0.64	0.25 – 1.63	0.351							0.40	0.15 – 1.08	0.071
Race/Ethnicity:Hispanic				0.33	0.09 – 1.21	0.095							0.89	0.29 – 2.78	0.847
Race/Ethnicity:Native American or Alaskan Native				0.41	0.10 – 1.65	0.210							3.58	1.12 – 11.46	0.032
Race/Ethnicity:Other Hispanic				0.25	0.10 – 0.66	0.005							0.48	0.18 – 1.32	0.156
Race/Ethnicity:White				0.34	0.17 – 0.70	0.003							1.85	0.85 – 4.01	0.120
Race/Ethnicity:White Hispanic				0.42	0.19 – 0.91	0.028							1.09	0.47 – 2.51	0.837
Income:$20K-$39K										5.92	2.27 – 15.42	<0.001	1.21	0.52 – 2.80	0.661
Income:$40K-$59K										3.09	1.34 – 7.10	0.008	1.75	0.80 – 3.82	0.161
Income:$60K-$79K										2.81	1.23 – 6.41	0.014	2.00	0.95 – 4.23	0.068
Income:$80K-$99K										1.91	0.86 – 4.26	0.113	2.56	1.20 – 5.46	0.015
Income:$100K-$149K										1.89	0.86 – 4.15	0.115	3.13	1.51 – 6.51	0.002
Income:>$150K										1.29	0.57 – 2.91	0.545	3.28	1.55 – 6.94	0.002
Observations	1303			1553			1901			1293			1238

tab_model(Q6_mod4,Q12_mod3,Q13_mod3,Q16_mod2,Q17_mod2,Q20_mod2,Q21_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?")) # 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?
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
COVID:Yes				0.56	0.41 – 0.77	<0.001	0.70	0.54 – 0.92	0.009				0.65	0.46 – 0.93	0.019	0.36	0.22 – 0.59	<0.001
New/Retained:New							0.67	0.53 – 0.85	0.001				0.78	0.58 – 1.05	0.101	0.31	0.17 – 0.56	<0.001
Gender:Male	1.47	1.18 – 1.84	0.001
Children:No Children													0.85	0.58 – 1.26	0.427
Children:Under 17													3.96	1.35 – 11.60	0.012
SW/FW:SW	0.44	0.33 – 0.59	<0.001	1.22	0.89 – 1.67	0.221	3.40	2.55 – 4.53	<0.001	2.20	1.59 – 3.04	<0.001	1.08	0.74 – 1.58	0.695	0.65	0.37 – 1.14	0.135
SW/FW:Both	1.68	1.30 – 2.18	<0.001	1.37	1.03 – 1.82	0.028	2.19	1.68 – 2.85	<0.001	1.92	1.42 – 2.60	<0.001	1.82	1.29 – 2.58	0.001	1.88	1.07 – 3.32	0.029
Age:25-34 yrs	0.88	0.54 – 1.43	0.602				1.21	0.76 – 1.91	0.420	1.43	0.81 – 2.54	0.222
Age:35-44 yrs	0.72	0.45 – 1.13	0.153				1.42	0.91 – 2.20	0.121	1.40	0.81 – 2.42	0.233
Age:45-54 yrs	0.71	0.45 – 1.12	0.138				1.20	0.77 – 1.86	0.416	1.53	0.88 – 2.64	0.130
Age:55-64 yrs	0.64	0.40 – 1.01	0.057				1.77	1.13 – 2.76	0.013	1.69	0.96 – 3.00	0.071
Age:65 yrs or older	0.50	0.30 – 0.83	0.007				2.74	1.66 – 4.52	<0.001	2.01	1.09 – 3.70	0.025
Race/Ethnicity:Black or African American				0.74	0.39 – 1.41	0.362				0.54	0.29 – 1.00	0.050							0.50	0.23 – 1.10	0.084
Race/Ethnicity:Hispanic				0.44	0.12 – 1.68	0.231				1.05	0.44 – 2.54	0.909							1.01	0.43 – 2.37	0.978
Race/Ethnicity:Native American or Alaskan Native				2.42	1.38 – 4.25	0.002				0.43	0.14 – 1.30	0.137							0.60	0.11 – 3.20	0.550
Race/Ethnicity:Other Hispanic				1.42	0.63 – 3.21	0.393				0.69	0.33 – 1.45	0.324							2.52	1.12 – 5.67	0.026
Race/Ethnicity:White				1.49	0.89 – 2.48	0.126				0.58	0.34 – 0.98	0.043							0.97	0.57 – 1.65	0.902
Race/Ethnicity:White Hispanic				0.69	0.39 – 1.22	0.199				0.59	0.33 – 1.07	0.082							1.15	0.62 – 2.11	0.660
Income:$20K-$39K										1.16	0.51 – 2.60	0.726	3.33	1.06 – 10.46	0.039				0.94	0.40 – 2.21	0.891
Income:$40K-$59K										2.01	0.98 – 4.13	0.056	4.34	1.41 – 13.36	0.011				1.37	0.62 – 3.04	0.436
Income:$60K-$79K										1.63	0.79 – 3.34	0.183	3.28	1.07 – 10.09	0.038				2.01	0.93 – 4.34	0.077
Income:$80K-$99K										1.99	0.96 – 4.11	0.065	2.80	0.97 – 8.11	0.058				1.88	0.85 – 4.19	0.121
Income:$100K-$149K										2.11	1.06 – 4.21	0.033	3.44	1.21 – 9.84	0.021				2.00	0.93 – 4.30	0.077
Income:>$150K										2.61	1.28 – 5.32	0.008	6.28	2.09 – 18.84	0.001				3.19	1.47 – 6.96	0.004
Days Fished:1 day|2-3 days	0.04	0.02 – 0.07	<0.001
Days Fished:2-3 days|4-10 days	0.32	0.20 – 0.51	<0.001
Days Fished:4-10 days|>10days	1.38	0.86 – 2.21	0.180
Catch/Keep Fish:Never|Rarely				0.03	0.02 – 0.06	<0.001	0.51	0.31 – 0.83	0.007
Catch/Keep Fish:Rarely|Sometimes				0.15	0.09 – 0.26	<0.001	1.65	1.01 – 2.69	0.048
Catch/Keep Fish:Sometimes|Most of the time				0.90	0.53 – 1.53	0.704	8.00	4.79 – 13.37	<0.001
Catch/Keep Fish:Most of the time|Always				8.93	5.10 – 15.62	<0.001	41.27	23.91 – 71.21	<0.001
Travel Time:<30 mins|30-60 mins										1.28	0.53 – 3.10	0.579	0.10	0.03 – 0.31	<0.001
Travel Time:30-60 mins|1-2 hrs										4.57	1.89 – 11.04	0.001	1.07	0.37 – 3.07	0.900
Travel Time:1-2 hrs|>2 hrs										14.53	5.96 – 35.41	<0.001
Travel Time:2-4 hrs|>4 hrs													1.37	1.01 – 1.85	0.040
Travel Time:1-2hrs|2-4 hrs													12.41	4.10 – 37.58	<0.001
Purchase License:Not likely|Somewhat likely																0.00	0.00 – 0.01	<0.001
Purchase License:Somewhat likely|Very likely																0.02	0.01 – 0.05	<0.001
Auto-renew:Extremely unlikely|Unlikely																			0.01	0.00 – 0.02	<0.001
Auto-renew:Unlikely|Not sure																			0.10	0.04 – 0.22	<0.001
Auto-renew:Not sure|Likely																			0.62	0.27 – 1.45	0.270
Auto-renew:Likely|Extremely likely																			2.32	0.98 – 5.48	0.054
Observations	1660			1565			1643			1237			1011			1717			1386