COVID Anglers Survey Analysis
Jen Granneman
5-16-2023
This is an analysis of a survey of anglers that recruited during the COVID-19 Pandemic. The data set includes a control group of anglers that recruited during the 3 years prior to the Pandemic. The purpose of this analysis is to identify the defining characteristics and preferences of these anglers to retain them.
library(dplyr)
library(pander)
library(knitr)
library(tidyverse)
library(readxl)
library(naniar)
library(vegan)
library(ggplot2)
library(rockchalk)
library(effectsize)
library(survey)
library(car)
library(ggeffects)
library(caret)
library(regclass)
library(generalhoslem)
library(effects)
library(faraway)
library(ecodist)
library(ape)
library(lmPerm)
library(factoextra)
library(kableExtra)
library(ade4)
library(BiodiversityR)
library(ggrepel)
library(corrplot)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)
## Adjust responses using weights
rev1_wts <- rev1 %>%
left_join(weights1, by = c("Gender","Age")) %>%
replace_na(list(weights = 1)) %>%
mutate(Gender = as.factor(Gender))
kable(weights1, caption = "Weights applied to the survey data", digits = 3)| Gender | Age | n_survey | freq_survey | n_random | freq_random | weights |
|---|---|---|---|---|---|---|
| Female | 18-24 | 39 | 0.021 | 1022 | 0.053 | 2.553 |
| Female | 25-34 | 138 | 0.073 | 1559 | 0.081 | 1.101 |
| Female | 35-44 | 171 | 0.091 | 1645 | 0.085 | 0.937 |
| Female | 45-54 | 161 | 0.086 | 1247 | 0.065 | 0.755 |
| Female | 55-64 | 153 | 0.081 | 767 | 0.040 | 0.488 |
| Female | 65 or older | 77 | 0.041 | 383 | 0.020 | 0.485 |
| Male | 18-24 | 51 | 0.027 | 2639 | 0.137 | 5.042 |
| Male | 25-34 | 124 | 0.066 | 2744 | 0.142 | 2.156 |
| Male | 35-44 | 237 | 0.126 | 2868 | 0.148 | 1.179 |
| Male | 45-54 | 305 | 0.162 | 2234 | 0.116 | 0.714 |
| Male | 55-64 | 254 | 0.135 | 1415 | 0.073 | 0.543 |
| Male | 65 or older | 172 | 0.091 | 792 | 0.041 | 0.449 |
Demographics of survey respondents
# Use survey design to incorporate non-response bias weights in age and gender
rev2_wts <- rev1_wts
svy_desn <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts) # Survey design
Gender <- svytable(~Gender, design = svy_desn)
Age <- svytable(~Age, design = svy_desn)
# Calculate proportions
Gender <- Gender %>%
data.frame() %>%
mutate(n = Freq) %>%
mutate(freq = prop.table(n)) %>%
dplyr::select(-Freq)
Age <- Age %>%
data.frame() %>%
mutate(n = Freq) %>%
mutate(freq = prop.table(n)) %>%
dplyr::select(-Freq)
COVID <- rev2_wts %>%
filter(!(is.na(COVID))) %>%
count(COVID) %>%
mutate(freq = prop.table(n))
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
SW_FW2 <- rev2_wts %>%
filter(!(is.na(SW_FW2))) %>%
count(SW_FW2) %>%
mutate(freq = prop.table(n)) # get proportionsAge %>%
kbl(caption = "Age",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | Age | n | freq |
|---|---|---|
| 18-24 | 357.72 | 0.19 |
| 25-34 | 420.27 | 0.22 |
| 35-44 | 442.73 | 0.23 |
| 45-54 | 343.18 | 0.18 |
| 55-64 | 222.61 | 0.12 |
| 65 or older | 114.49 | 0.06 |
Gender %>%
kbl(caption = "Gender",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | Gender | n | freq |
|---|---|---|
| Female | 650.33 | 0.34 |
| Male | 1247.67 | 0.66 |
Children %>%
kbl(caption = "Children",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | 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 | n | freq |
|---|---|---|
| Between $100,000 and $149,999 | 332 | 0.23 |
| 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 |
| Less than $20,000 | 70 | 0.05 |
| Over $150,000 | 298 | 0.20 |
Race_Ethn %>%
kbl(caption = "Race and Ethnicity",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | 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 |
COVID %>%
kbl(caption = "Prompted to fish by COVID",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | COVID | n | freq |
|---|---|---|
| No, I was already planning to go fishing. | 1583 | 0.81 |
| Yes | 381 | 0.19 |
New.Old %>%
kbl(caption = "Retained or New Angler during pandemic",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | New.Old | n | freq |
|---|---|---|
| New | 1517 | 0.68 |
| Old | 721 | 0.32 |
SW_FW2 %>%
kbl(caption = "Saltwater or Freshwater Angler",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | SW_FW2 | n | freq |
|---|---|---|
| FW | 802 | 0.43 |
| SW | 448 | 0.24 |
| SW.FW | 631 | 0.34 |
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 monthsResponses are summarized into likely or unlikely. Most anglers (99%) are likely to go fishing in the next 6 months, further analysis is not needed.
Q2. Why are you unlikely to go fishing in the next 6 months?
Q2_responses <- rev1_wts %>%
mutate(Q2_all = as.factor(case_when(Q2B == "It’s too hard to find a place to fish." ~ "Hard to find fishing place",
Q2C == "No one in my family or friend circle wants to go with me." ~ "No one to fish with",
Q2D == "I don’t know what to do with fish when I catch them." ~ "Don't know what to do with fish",
Q2E == "I decided I just don’t like fishing." ~ "Don't like fishing",
str_length(Q2Other)>0 ~ "Other"))) %>%
filter(!is.na(Q2_all))
ggplot(Q2_responses,aes(Q2_all)) + geom_histogram(stat = "count") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8),axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + labs(y = "Counts", x = "Q2. Reasons Unlikely to Go Fishing") # histogram counts of levels
Only 19 responses, not enough data for robust analysis
Q3. Did you go fishing between mid-March 2020 and today?
svy_desn_Q3 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts) # Survey design
Q3_mod1 <- svyglm(Q3_Fished ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age, 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_mod2$survey.design)) # Evaluate variable significance
# Estimate and odds ratio table
Q3_coef_table <- coef(summary(Q3_mod3))
Q3_AOR <- standardize_parameters(Q3_mod3, exp = TRUE) # Odds_ratio and 95% CI
Q3_summary_table <- cbind(Q3_coef_table,Odds_Ratio = Q3_AOR$Std_Odds_Ratio,CI_low = Q3_AOR$CI_low,CI_high = Q3_AOR$CI_high)
# Check model accuracy
logitgof(obs = Q3_mod3$y, fitted(Q3_mod3)) #run hoslem test to test goodness of fit
glm.probs <- predict(Q3_mod3,type = "response") #Predict() assigns fitted model values for each participants based on the model
glm.pred <- as.factor(ifelse(glm.probs > 0.5, 1, 0)) #If prediction is above 50% assign 1 and assign 0 if prediction is below 50%
prop.table(table(glm.pred == Q3_mod3$y)) # Correct assignment 88% of the time
Q3_age_xtab <- xtabs(~Q3_Fished + Age, data = rev2_wts)
Q3_AOR_18 <- 1/0.27
# Survey response proportions
Q3_prop <- rev2_wts %>%
filter(!(is.na(Q3_Fished))) %>% # remove NA values from Gender & Age
count(Q3_Fished) %>%
mutate(freq = prop.table(n)) # get proportions# Table of responses
Q3_prop %>%
kbl(caption = "Q3 Survey Responses",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") # 88% went fishing in the last 6 months| Q3_Fished | n | freq |
|---|---|---|
| No | 256 | 0.12 |
| Yes | 1965 | 0.88 |
# Histogram of responses
ggplot(model.frame(Q3_mod3), aes(Q3_Fished, fill = Age)) + geom_histogram(stat = "count") + theme_classic() + 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 = "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")| 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.70 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?
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, design=svy_desn_Q4, family=quasibinomial)
summary(Q4_mod2)
Q4_mod3 <- svyglm(Q4_COVID_Fished ~ Age + Children + Gender, 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_mod2$survey.design)) # Evaluate variable significance
# Estimate and odds ratio table
Q4_coef_table <- coef(summary(Q4_mod3))
Q4_AOR <- standardize_parameters(Q4_mod3, exp = TRUE) # Odds_ratio and 95% CI
Q4_summary_table <- cbind(Q4_coef_table,Odds_Ratio = Q4_AOR$Std_Odds_Ratio,CI_low = Q4_AOR$CI_low,CI_high = Q4_AOR$CI_high)
# Check model accuracy
logitgof(obs = Q4_mod3$y, fitted(Q4_mod3)) #run hoslem test to test goodness of fit
glm.probs <- predict(Q4_mod3,type = "response") #Predict() assigns fitted model values for each participants based on the model
glm.pred <- as.factor(ifelse(glm.probs > 0.5, 1, 0)) #If prediction is above 50% assign 1 and assign 0 if prediction is below 50%
prop.table(table(glm.pred == Q4_mod3$y))
Q4_xtab <- xtabs(~Q4_COVID_Fished + Age, data = rev2_wts)
Q4_xtab <- Q4_xtab%>%ftable(row.vars=c("Age"))
# Survey response proportions
Q4_prop <- rev2_wts %>%
filter(!(is.na(Q4_COVID_Fished))) %>% # remove NA values from Gender & Age
count(Q4_COVID_Fished) %>%
mutate(freq = prop.table(n)) # get proportions
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_COVID_Fished | n | freq |
|---|---|---|
| No, I was already planning to go fishing. | 1583 | 0.81 |
| Yes | 381 | 0.19 |
# Histogram of responses
ggplot(model.frame(Q4_mod3), aes(Q4_COVID_Fished, fill = Age)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q4. Fished due to COVID") + facet_wrap(~Gender, ncol = 2)# 
ggplot(model.frame(Q4_mod3), aes(Q4_COVID_Fished, fill = Children)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q4. Fished due to COVID")
# Model odds ratio results
Q4_summary_table %>%
kbl(caption = "Q4 Model Odds Ratios",digits = 3) %>%
kable_classic(full_width = F, html_font = "Cambria")| Estimate | Std. Error | t value | Pr(>|t|) | Odds_Ratio | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| (Intercept) | -1.837 | 0.334 | -5.509 | 0.000 | 0.159 | 0.083 | 0.306 |
| Age25-34 | 0.701 | 0.333 | 2.107 | 0.035 | 2.016 | 1.050 | 3.871 |
| Age35-44 | 0.522 | 0.331 | 1.578 | 0.115 | 1.685 | 0.881 | 3.223 |
| Age45-54 | 0.399 | 0.344 | 1.159 | 0.247 | 1.490 | 0.759 | 2.925 |
| Age55-64 | 0.428 | 0.339 | 1.262 | 0.207 | 1.534 | 0.789 | 2.982 |
| Age65 or older | -0.110 | 0.398 | -0.277 | 0.782 | 0.896 | 0.410 | 1.955 |
| ChildrenNo children | 0.767 | 0.187 | 4.100 | 0.000 | 2.154 | 1.492 | 3.110 |
| ChildrenUnder 17 | 0.325 | 0.267 | 1.218 | 0.224 | 1.384 | 0.820 | 2.336 |
| GenderMale | -0.654 | 0.177 | -3.692 | 0.000 | 0.520 | 0.367 | 0.736 |
# Effect plot
plot(Effect(focal.predictors = c("Children","Age","Gender"),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.01 times less likely for young anglers (25-34) compared to other ages, 2.15 times less likely for anglers with no children compared to those with adult children, 1.92 times less likely for male anglers compared to females.
Q5. How important were each of the following reasons why you went fishing? Please select no more than 2 of these for the “Very Important” column.
# Clean Q5 survey answers
rev2_wts <- rev2_wts %>%
mutate_at(vars(Q5A:Q5I), list(~recode_factor(.,"Not Important" = "1","Somewhat Important" = "2","Important" = "3","Very Important (please limit to 2 reasons)" = "4"))) %>%
mutate(across(Q5A:Q5I, ~as.numeric(as.character(.)))) %>%
mutate_at(vars(Q5A:Q5I), list(~.*weights)) %>%
mutate(Food = Q5A,
Sport = Q5B,
Outdoors = Q5C,
Nature = Q5D,
Friends = Q5E,
Children1 = Q5F,
Stress = Q5G,
Relax = Q5H,
Get_away = Q5I)
# Remove rows with missing independent variable data
Q5_df1 <- rev2_wts %>%
filter_at(vars(Gender,Income,Age,Children,COVID,New.Old,Race_Ethn,SW_FW2),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 CCA
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","pre-COVID","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)# 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.005 **
## RDA4 1 0.000222 2.7263 0.320
## RDA5 1 0.000114 1.3991 0.879
## 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 ' ' 1envfit(Q5_rda_mod2,Q5_indep3) # select variables at alpha = 0.001 for next model##
## ***VECTORS
##
## RDA1 RDA2 r2 Pr(>r)
## Gende.Female -0.43263 0.90157 0.0911 0.001
## Gende.Male 0.43263 -0.90157 0.0911 0.001
## Age.18.24 -0.97516 0.22150 0.0201 0.001
## Age.25.34 -0.97536 0.22062 0.0491 0.001
## Age.65.or.older 0.81081 -0.58531 0.0181 0.001
## Child.18.or.over 0.99540 -0.09579 0.1322 0.001
## Child.No.children -0.99675 0.08056 0.2559 0.001
## Child.Under.17 0.99952 -0.03097 0.0263 0.001
## COVID.No..I.was.already.planning.to.go.fishing. 0.89293 -0.45019 0.0242 0.001
## COVID.Yes -0.89293 0.45019 0.0242 0.001
## New.O.New -0.65496 0.75567 0.0263 0.001
## New.O.Old 0.65496 -0.75567 0.0263 0.001
##
## 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 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999# 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 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.
Q6 Since mid-March 2020, how many days have you gone fishing?
# Order the factor
rev2_wts$Q6_DaysFished <- as.factor(rev2_wts$Q6_DaysFished)
rev2_wts$Q6_DaysFished <- ordered(rev2_wts$Q6_DaysFished, levels = c("Once (1 day)", "2–3 days", "4–10 days","More than 10 days"))
svy_desn_Q6 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts,na.rm = TRUE) # New survey design
# Ordinal regression
Q6_mod1 <- svyolr(Q6_DaysFished ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age, 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)
# 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_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(model.frame(Q6_mod4), aes(Q6_DaysFished, fill = Age)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q6. Days Fished in Last 6 Months") + facet_wrap(~Gender, ncol = 2) 
# 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| 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+).
Q8 From mid-March 2020 until today, did you go fishing from a boat?
rev2_wts <- rev1_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 boatresponses_Q8_Y_N %>%
kbl(caption = "Q8A. Did you fish from a boat?",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") # 81% were already planning to go fishing| 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") # 81% were already planning to go fishing| Boat_Own | n | freq |
|---|---|---|
| Borrow | 674 | 0.64 |
| Own | 379 | 0.36 |
Anglers fished from a boat 56% of the time and 36% of those that fished from a boat owned the boat.
Q9. The last time you went fishing, was fishing the main purpose of the trip, or was it just one of the activities but not the main purpose of the trip?
rev2_wts$Q9_Purpose <- recode_factor(rev2_wts$Q9_Purpose,"Fishing was one activity, but not the main purpose of the trip." = "No","Fishing was the main purpose of the trip." = "Yes")
svy_desn_Q9 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts) # Survey design
Q9_mod1 <- svyglm(Q9_Purpose ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age, 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
Q4_AOR_65 <- 1/0.21 # 4.76 x
Q4_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_Purpose | n | freq |
|---|---|---|
| No | 468 | 0.25 |
| Yes | 1398 | 0.75 |
# Histogram of responses
ggplot(model.frame(Q9_mod2), aes(Q9_Purpose, fill = Age)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q9. Was Fishing the Main Purpose?") + facet_wrap(~Gender, ncol = 2)
# Summary table
Q9_summary_table %>%
kbl(caption = "Q9 Model Odds Ratios",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | Estimate | Std. Error | t value | Pr(>|t|) | Odds_Ratio | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| (Intercept) | 1.78 | 0.43 | 4.10 | 0.00 | 5.92 | 2.53 | 13.85 |
| 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.14 | 0.35 | 3.28 | 0.00 | 3.14 | 1.58 | 6.21 |
| IncomeBetween $40,000 and $59,999 | 0.49 | 0.27 | 1.81 | 0.07 | 1.64 | 0.96 | 2.80 |
| IncomeBetween $60,000 and $79,999 | 0.40 | 0.26 | 1.54 | 0.12 | 1.49 | 0.90 | 2.47 |
| IncomeBetween $80,000 and $99,999 | 0.01 | 0.24 | 0.06 | 0.95 | 1.01 | 0.63 | 1.63 |
| IncomeLess than $20,000 | -0.63 | 0.40 | -1.58 | 0.11 | 0.53 | 0.24 | 1.17 |
| IncomeOver $150,000 | -0.38 | 0.23 | -1.69 | 0.09 | 0.68 | 0.44 | 1.06 |
# 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, 3.14 times more likely for anglers with incomes between $20-$39K compared to anglers with higher incomes.
Q10. Who do you typically fish with?
# Evaluate frequencies by categorical responses
Q10_df1 <- rev2_wts %>%
dplyr::select(c("Q10A":"Q10J")) %>%
pivot_longer(cols = "Q10A":"Q10J",names_to = "Q10_answers",values_to = "Agreement") %>%
group_by(Q10_answers,Agreement) %>%
count() %>%
na.omit()
freq = prop.table(Q10_df1$n)
Q10_df1 <- data.frame(cbind(Q10_df1,freq=freq))
# Format data for analysis
Q10_df2 <- rev2_wts %>%
dplyr::select(Q10A,Q10F,Q10G,Q10H,Q10I,Q10J,Gender,Age,Race_Ethn,Income,SW_FW2,New.Old,Children,COVID,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()
# 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_tbl <- svytable(~Gender + Agreement, svy_desn_Q10)
Q10_Age_tbl <- svytable(~Age + Agreement, svy_desn_Q10)
Q10_Race_Ethn_tbl <- svytable(~Race_Ethn + Agreement, svy_desn_Q10)
Q10_Income_tbl <- svytable(~Income + Agreement, svy_desn_Q10)
Q10_SW_FW2_tbl <- svytable(~SW_FW2 + Agreement, svy_desn_Q10)
Q10_COVID_tbl <- svytable(~COVID + Agreement, svy_desn_Q10)
Q10_Children_tbl <- svytable(~Children + Agreement, svy_desn_Q10)# 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?") 
# Chi-square test results
Q10_Gender_tbl %>%
kbl(caption = "Q10 Who do you fish with?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | 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 | 159 | 239 | 281 | 99 | 98 |
| Male | 184 | 265 | 350 | 590 | 177 | 135 |
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-07corrplot(Q10_Gender_chisq$residuals, is.cor = FALSE)
Q10_Age_tbl %>%
kbl(caption = "Q10 Who do you fish with?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | 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 | 10 | 120 | 231 | 10 | 21 |
| 25-34 | 57 | 125 | 143 | 239 | 24 | 5 |
| 35-44 | 39 | 185 | 143 | 169 | 120 | 29 |
| 45-54 | 37 | 69 | 88 | 132 | 92 | 67 |
| 55-64 | 30 | 24 | 66 | 70 | 23 | 69 |
| 65 or older | 10 | 13 | 32 | 30 | 6 | 37 |
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-16corrplot(Q10_Age_chisq$residuals, is.cor = FALSE)
Q10_Race_Ethn_tbl %>%
kbl(caption = "Q10 Who do you fish with?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | 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 | 11 | 20 | 35 | 6 | 2 |
| Black or African American | 8 | 18 | 19 | 33 | 4 | 4 |
| Hispanic | 1 | 14 | 22 | 21 | 9 | 6 |
| Native American or Alaskan Native | 9 | 14 | 12 | 24 | 10 | 7 |
| Other Hispanic | 15 | 33 | 29 | 60 | 25 | 20 |
| White | 135 | 236 | 342 | 474 | 160 | 141 |
| White Hispanic | 32 | 60 | 104 | 171 | 40 | 36 |
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.00247corrplot(Q10_Race_Ethn_chisq$residuals, is.cor = FALSE)
Q10_Income_tbl %>%
kbl(caption = "Q10 Who do you fish with?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | 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 | |
|---|---|---|---|---|---|---|
| Between $100,000 and $149,999 | 28 | 69 | 91 | 137 | 59 | 43 |
| Between $20,000 and $39,999 | 27 | 42 | 59 | 94 | 14 | 14 |
| Between $40,000 and $59,999 | 35 | 42 | 62 | 88 | 22 | 17 |
| Between $60,000 and $79,999 | 31 | 63 | 68 | 117 | 29 | 27 |
| Between $80,000 and $99,999 | 25 | 52 | 59 | 85 | 36 | 29 |
| Less than $20,000 | 21 | 12 | 39 | 36 | 7 | 4 |
| Over $150,000 | 27 | 63 | 71 | 105 | 49 | 38 |
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.02125corrplot(Q10_Income_chisq$residuals, is.cor = FALSE)
Q10_SW_FW2_tbl %>%
kbl(caption = "Q10 Who do you fish with?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | 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 | 218 | 253 | 345 | 107 | 92 |
| SW | 27 | 84 | 137 | 209 | 68 | 68 |
| SW.FW | 77 | 172 | 251 | 392 | 127 | 101 |
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.0001434corrplot(Q10_SW_FW2_chisq$residuals, is.cor = FALSE)
Q10_COVID_tbl %>%
kbl(caption = "Q10 Who do you fish with?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | 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 | 387 | 532 | 765 | 253 | 230 |
| Yes | 54 | 88 | 116 | 186 | 49 | 33 |
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.04205corrplot(Q10_COVID_chisq$residuals, is.cor = FALSE)
Q10_Children_tbl %>%
kbl(caption = "Q10 Who do you fish with?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | 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 | 83 | 182 | 187 | 132 | 179 |
| No children | 126 | 13 | 207 | 425 | 18 | 30 |
| Under 17 | 24 | 118 | 85 | 119 | 28 | 14 |
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-16corrplot(Q10_Children_chisq$residuals, is.cor = FALSE)
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 <- 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)
Q11_rda_mod2 <- rda(Q11_dep2 ~ ., data = Q11_indep3)
envfit(Q11_rda_mod2,Q11_indep3) # select variables at alpha = 0.001 for next model
Q11_indep4 <- Q11_indep2 %>%
dplyr::select(Age.18.24, Age.25.34,Age.55.64, Age.65.or.older, Child.18.or.over, Child.No.children)
Q11_rda_mod3 <- rda(Q11_dep2 ~ ., data = Q11_indep4)
envfit(Q11_rda_mod3,Q11_indep4) # select variables at alpha = 0.001 for next model
Q11_rda_mod3$call
RsquareAdj(Q11_rda_mod3)
sqrt(vif.cca(Q11_rda_mod3))
# 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"),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)# Reorganize data to plot
Q11_a <- rev1_wts[c(which(colnames(rev1_wts)=="Q11A"):which(colnames(rev1_wts)=="Q11D"))]
colnames(Q11_a) <- c("Enjoyed_fishing","Too_expensive","Fishing_rules_confusing","Enjoyed_surroundings")
Q11_b <- Q11_a %>%
pivot_longer(cols = Enjoyed_fishing:Enjoyed_surroundings,names_to = "Q11_answers",values_to = "Agreement")
Q11_c <- Q11_b %>%
group_by(Q11_answers,Agreement) %>%
count() %>%
na.omit()
Q11_c$Agreement <- ordered(Q11_c$Agreement,levels = c("Strongly Agree", "Agree", "Disagree","Strongly Disagree") )
ggplot(Q11_c, aes(fill=Agreement, y=n, x=Q11_answers)) + geom_bar(position="fill", stat="identity") + theme_classic() + theme(axis.title = element_text(size=16), axis.text = element_text(size = 8),axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + labs(y = "Percent", x = "Q11. Agreement on Fishing Statements") 
# PERMANOVA
anova.cca(Q11_rda_mod3, step = 1000, by = "axis",permutations = 999)## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = Q11_dep2 ~ Age.18.24 + Age.25.34 + Age.55.64 + Age.65.or.older + Child.18.or.over + Child.No.children, data = Q11_indep4)
## Df Variance F Pr(>F)
## RDA1 1 1.91738 879.8543 0.001 ***
## RDA2 1 0.00205 0.9425 0.993
## RDA3 1 0.00115 0.5297 0.998
## RDA4 1 0.00044 0.2038 1.000
## Residual 954 2.07896
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1envfit(Q11_rda_mod3,Q11_indep4) # select variables at alpha = 0.001 for next model##
## ***VECTORS
##
## RDA1 RDA2 r2 Pr(>r)
## Age.18.24 0.999920 -0.012517 0.5785 0.001 ***
## Age.25.34 0.999850 -0.017438 0.0709 0.001 ***
## Age.55.64 -0.999890 0.014526 0.0948 0.001 ***
## Age.65.or.older -0.999900 0.014469 0.0491 0.001 ***
## Child.18.or.over -0.999940 0.011054 0.1972 0.001 ***
## Child.No.children 0.999880 -0.015624 0.1200 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999# 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 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.
Q12. How often did you catch fish during your trips?
# Order the factor
rev2_wts$Q12_CatchOften <- as.factor(rev2_wts$Q12_CatchOften)
rev2_wts$Q12_CatchOften <- ordered(rev2_wts$Q12_CatchOften, levels = c("Never","Rarely", "Sometimes","Most of the time","Always"))
svy_desn_Q12 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts,na.rm = TRUE) # New survey design
# Ordinal regression
Q12_mod1 <- svyolr(Q12_CatchOften ~ SW_FW2 + Age + Gender + Income + Children + New.Old + COVID + Race_Ethn, 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
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_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(model.frame(Q12_mod3), aes(Q12_CatchOften, fill = SW_FW2)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q12. How often do you catch fish?") + facet_wrap(~COVID, ncol = 2)
# Odds Ratio summary table
Q12_summary_table %>%
kbl(caption = "Q12 Model Odds Ratios",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | 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 already planning to fish pre-COVID compared to anglers that recruited during COVID, 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$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, 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# 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_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(model.frame(Q13_mod3), aes(Q13_KeepOften, fill = SW_FW2)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q13. How often do you keep fish?") + facet_wrap(~New.Old, ncol = 2)
# Odds Ratio summary table
Q13_summary_table %>%
kbl(caption = "Q13 Model Odds Ratios",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | Value | Std. Error | t value | p value | Odds_Ratio | |
|---|---|---|---|---|---|
| New.OldOld | 0.43 | 0.12 | 3.69 | 0.00 | 1.54 |
| SW_FW2SW | 1.21 | 0.15 | 8.30 | 0.00 | 3.36 |
| SW_FW2SW.FW | 0.78 | 0.14 | 5.79 | 0.00 | 2.18 |
| Age25-34 | 0.16 | 0.23 | 0.69 | 0.49 | 1.17 |
| Age35-44 | 0.32 | 0.22 | 1.45 | 0.15 | 1.38 |
| Age45-54 | 0.18 | 0.22 | 0.81 | 0.42 | 1.20 |
| Age55-64 | 0.57 | 0.23 | 2.51 | 0.01 | 1.77 |
| Age65 or older | 1.02 | 0.25 | 4.04 | 0.00 | 2.79 |
| Never|Rarely | -0.20 | 0.23 | -0.89 | 0.37 | 0.82 |
| Rarely|Sometimes | 0.96 | 0.23 | 4.12 | 0.00 | 2.61 |
| Sometimes|Most of the time | 2.54 | 0.25 | 10.15 | 0.00 | 12.63 |
| Most of the time|Always | 4.18 | 0.27 | 15.32 | 0.00 | 65.05 |
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.54 times more likely for for anglers that had licenses prior to the pandemic, 3.36 times more likely for saltwater anglers compared to freshwater anglers, 2.18 times more likely for anglers that fish in both freshwater and saltwater compared to freshwater anglers, 1.77-2.79 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
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")
# 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
# Contingency tables for significant factors
Q14_Gender_tbl <- svytable(~Gender + Agreement, svy_desn_Q14)
Q14_Age_tbl <- svytable(~Age + Agreement, svy_desn_Q14)
Q14_Race_Ethn_tbl <- svytable(~Race_Ethn + Agreement, svy_desn_Q14)
Q14_SW_FW2_tbl <- svytable(~SW_FW2 + Agreement, svy_desn_Q14)
Q14_New.Old_tbl <- svytable(~New.Old + Agreement, svy_desn_Q14)
Q14_COVID_tbl <- svytable(~COVID + Agreement, svy_desn_Q14)
Q14_Children_tbl <- svytable(~Children + Agreement, svy_desn_Q14)# 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?") 
# Chi-square test results
Q14_Gender_tbl %>%
kbl(caption = "Q14 Anything about fishing you didn't enjoy?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | Enjoyed everything | Finding close fishing spot | Needed equipment | Techniques for casting/reeling | Children had trouble | Didn’t catch anything | |
|---|---|---|---|---|---|---|
| Female | 285 | 125 | 27 | 102 | 32 | 128 |
| Male | 433 | 367 | 114 | 174 | 60 | 240 |
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-05corrplot(Q14_Gender_chisq$residuals, is.cor = FALSE)
Q14_Age_tbl %>%
kbl(caption = "Q14 Anything about fishing you didn't enjoy?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | Enjoyed everything | Finding close fishing spot | Needed equipment | Techniques for casting/reeling | Children had trouble | Didn’t catch anything | |
|---|---|---|---|---|---|---|
| 18-24 | 89 | 144 | 41 | 63 | 5 | 101 |
| 25-34 | 146 | 131 | 34 | 98 | 33 | 79 |
| 35-44 | 175 | 101 | 36 | 56 | 37 | 89 |
| 45-54 | 153 | 68 | 18 | 33 | 11 | 52 |
| 55-64 | 110 | 33 | 9 | 19 | 4 | 28 |
| 65 or older | 48 | 15 | 3 | 7 | 3 | 15 |
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-14corrplot(Q14_Age_chisq$residuals, is.cor = FALSE)
Q14_Race_Ethn_tbl %>%
kbl(caption = "Q14 Anything about fishing you didn't enjoy?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | Enjoyed everything | Finding close fishing spot | Needed equipment | Techniques for casting/reeling | Children had trouble | Didn’t catch anything | |
|---|---|---|---|---|---|---|
| Asian or Pacific Islander | 19 | 24 | 9 | 22 | 2 | 14 |
| Black or African American | 20 | 27 | 13 | 17 | 7 | 23 |
| Hispanic | 14 | 15 | 3 | 12 | 5 | 10 |
| Native American or Alaskan Native | 22 | 12 | 7 | 12 | 3 | 9 |
| Other Hispanic | 33 | 44 | 6 | 16 | 4 | 9 |
| White | 453 | 239 | 71 | 127 | 58 | 197 |
| White Hispanic | 112 | 91 | 31 | 67 | 8 | 73 |
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.001846corrplot(Q14_Race_Ethn_chisq$residuals, is.cor = FALSE)
Q14_SW_FW2_tbl %>%
kbl(caption = "Q14 Anything about fishing you didn't enjoy?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | Enjoyed everything | Finding close fishing spot | Needed equipment | Techniques for casting/reeling | Children had trouble | Didn’t catch anything | |
|---|---|---|---|---|---|---|
| FW | 309 | 230 | 56 | 137 | 44 | 186 |
| SW | 215 | 68 | 22 | 51 | 15 | 59 |
| SW.FW | 261 | 214 | 68 | 111 | 42 | 152 |
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-05corrplot(Q14_SW_FW2_chisq$residuals, is.cor = FALSE)
Q14_COVID_tbl %>%
kbl(caption = "Q14 Anything about fishing you didn't enjoy?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | 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 | 390 | 114 | 180 | 64 | 282 |
| Yes | 103 | 128 | 32 | 120 | 37 | 116 |
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-13corrplot(Q14_COVID_chisq$residuals, is.cor = FALSE)
Q14_Children_tbl %>%
kbl(caption = "Anything about fishing you didn't enjoy?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | Enjoyed everything | Finding close fishing spot | Needed equipment | Techniques for casting/reeling | Children had trouble | Didn’t catch anything | |
|---|---|---|---|---|---|---|
| 18 or over | 270 | 97 | 31 | 45 | 17 | 92 |
| No children | 221 | 233 | 74 | 152 | 5 | 151 |
| Under 17 | 81 | 81 | 19 | 32 | 31 | 57 |
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-16corrplot(Q14_Children_chisq$residuals, is.cor = FALSE)
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.
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$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
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_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(model.frame(Q16_mod2), aes(Q16_TravelTime, fill = Age)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q16. How long did it take you to get to your fishing spot?") + facet_wrap(~SW_FW2, ncol = 2)
ggplot(model.frame(Q16_mod2), aes(Q16_TravelTime, fill = Income)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q16. How long did it take you to get to your fishing spot?") 
ggplot(model.frame(Q16_mod2), aes(Q16_TravelTime, fill = Race_Ethn)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q16. How long did it take you to get to your fishing spot?") 
# Summary table
Q16_summary_table %>%
kbl(caption = "Q16 Model Odds Ratios",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | Value | Std. Error | t value | p value | Odds_Ratio | |
|---|---|---|---|---|---|
| IncomeBetween $20,000 and $39,999 | -0.60 | 0.28 | -2.18 | 0.03 | 0.55 |
| IncomeBetween $40,000 and $59,999 | -0.05 | 0.20 | -0.24 | 0.81 | 0.95 |
| IncomeBetween $60,000 and $79,999 | -0.26 | 0.20 | -1.30 | 0.20 | 0.77 |
| IncomeBetween $80,000 and $99,999 | -0.06 | 0.21 | -0.30 | 0.77 | 0.94 |
| IncomeLess than $20,000 | -0.75 | 0.35 | -2.13 | 0.03 | 0.47 |
| IncomeOver $150,000 | 0.21 | 0.19 | 1.08 | 0.28 | 1.23 |
| 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.50 | 0.39 | -1.29 | 0.20 | 0.61 |
| Between 30 and 60 minutes|Between 1 and 2 hours | 0.77 | 0.38 | 2.00 | 0.04 | 2.16 |
| Between 1 and 2 hours|More than 2 hours | 1.93 | 0.39 | 4.93 | 0.00 | 6.87 |
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.13 times more likely for anglers with an income between $100K-$150K 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$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, design=svy_desn_Q17)
summary(Q17_mod2) # examine model results
VIF(Q17_mod2) # collinearity
Q17_mod3 <- svyolr(Q17_MaxTravel ~ Income + COVID + SW_FW2, 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
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_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(model.frame(Q17_mod2), aes(Q17_MaxTravel, fill = Income)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q17. How far are you willing to travel to go fishing?") + facet_wrap(~SW_FW2, ncol = 2)
ggplot(model.frame(Q17_mod2), aes(Q17_MaxTravel, fill = COVID)) + geom_histogram(stat = "count") + theme_classic() + 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 = "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") | Value | Std. Error | t value | p value | Odds_ratio | |
|---|---|---|---|---|---|
| IncomeBetween $20,000 and $39,999 | -0.02 | 0.25 | -0.08 | 0.94 | 0.98 |
| IncomeBetween $40,000 and $59,999 | 0.17 | 0.20 | 0.85 | 0.40 | 1.19 |
| IncomeBetween $60,000 and $79,999 | 0.12 | 0.22 | 0.52 | 0.60 | 1.12 |
| IncomeBetween $80,000 and $99,999 | -0.17 | 0.16 | -1.04 | 0.30 | 0.84 |
| IncomeLess than $20,000 | -0.95 | 0.45 | -2.12 | 0.03 | 0.39 |
| IncomeOver $150,000 | 0.51 | 0.18 | 2.84 | 0.00 | 1.66 |
| COVIDYes | -0.41 | 0.15 | -2.70 | 0.01 | 0.66 |
| SW_FW2SW | 0.10 | 0.17 | 0.62 | 0.54 | 1.11 |
| SW_FW2SW.FW | 0.67 | 0.15 | 4.49 | 0.00 | 1.96 |
| Less than 30 minutes|Between 30 and 60 minutes | -3.32 | 0.23 | -14.53 | 0.00 | 0.04 |
| Between 30 and 60 minutes|Between 1 and 2 hours | -1.03 | 0.14 | -7.10 | 0.00 | 0.36 |
| Between 1 and 2 hours|Between 2 and 4 hours | 0.30 | 0.14 | 2.10 | 0.04 | 1.35 |
| Between 2 and 4 hours|More than 4 hours | 1.40 | 0.15 | 9.35 | 0.00 | 4.07 |
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.56 times more likely for anglers with an income between $100K-$150K compared to anglers with the lowest incomes (<$20K), 1.66 times more likely for anglers with an income over $150K compared to anglers with lower incomes, 1.52 times more likely for anglers that had licenses prior to the pandemic, and 1.96 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 <- 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)
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)
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=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"),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)# 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, data = Q18_indep4)
## Df Variance F Pr(>F)
## RDA1 1 0.01240 26.2410 0.001 ***
## RDA2 1 0.00486 10.2844 0.002 **
## RDA3 1 0.00260 5.5001 0.103
## RDA4 1 0.00251 5.3210 0.087 .
## RDA5 1 0.00099 2.1053 0.877
## RDA6 1 0.00037 0.7828 1.000
## RDA7 1 0.00015 0.3210 1.000
## RDA8 1 0.00011 0.2434 1.000
## RDA9 1 0.00008 0.1719 1.000
## Residual 949 0.44843
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1envfit(Q18_rda_mod3,Q18_indep4) # select variables at alpha = 0.001 for next model##
## ***VECTORS
##
## RDA1 RDA2 r2 Pr(>r)
## Incom.Over..150.000 0.93128 0.36431 0.0214 0.001 ***
## Age.25.34 -0.99940 -0.03474 0.0270 0.001 ***
## Age.65.or.older 0.99711 0.07597 0.0192 0.001 ***
## Child.No.children -0.50675 -0.86209 0.0210 0.001 ***
## Child.Under.17 -0.87368 0.48650 0.0501 0.001 ***
## New.O.New -0.79407 -0.60783 0.0267 0.001 ***
## New.O.Old 0.79407 0.60783 0.0267 0.001 ***
## Race_.Black.or.African.American -0.61014 -0.79229 0.0138 0.001 ***
## Race_.White 0.99582 0.09129 0.0476 0.001 ***
## Race_.White.Hispanic -0.98298 0.18369 0.0180 0.001 ***
## SW_FW.SW 0.28128 0.95962 0.0377 0.001 ***
## SW_FW.FW -0.24117 -0.97048 0.0175 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999# 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 that recruited during the pandemic and fish in freshwater were more likely to to agree that fish stockings, having peaceful areas, ADA accessibility, and shoreline access were important.
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, white, saltwater anglers, and 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.
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:
# 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)# 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.321
## RDA4 1 0.00136 2.3558 0.529
## RDA5 1 0.00112 1.9393 0.576
## RDA6 1 0.00040 0.6848 0.993
## 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 ' ' 1envfit(Q19_rda_mod2,Q19_indep3) # select variables at alpha = 0.001 for next model##
## ***VECTORS
##
## RDA1 RDA2 r2 Pr(>r)
## Age.25.34 -0.45223 0.89190 0.0232 0.001 ***
## Age.65.or.older 0.58642 -0.81001 0.0127 0.001 ***
## Child.18.or.over 0.42683 -0.90433 0.0306 0.001 ***
## Child.No.children 0.25079 0.96804 0.0602 0.001 ***
## Child.Under.17 -0.81407 -0.58077 0.0594 0.001 ***
## New.O.New -0.27017 0.96281 0.0399 0.001 ***
## New.O.Old 0.27017 -0.96281 0.0399 0.001 ***
## Race_.White 0.87387 0.48617 0.0297 0.001 ***
## Race_.White.Hispanic -0.95584 -0.29390 0.0157 0.001 ***
## SW_FW.SW 0.02461 -0.99970 0.0341 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999plot(Q19_rda_mod2,display = c("sp","bp"),const = -3,scaling = "symmetric")
# 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.
Q20. How likely are you to purchase another fishing license when your current license expires?
# Order the factor
rev2_wts$Q20_NewLicense <- recode_factor(rev2_wts$Q20_NewLicense,"Not likely at all" = "Not likely","Not very likely" = "Not likely")
rev2_wts$Q20_NewLicense <- ordered(rev2_wts$Q20_NewLicense, levels = c("Not likely", "Somewhat likely", "Very likely"))
svy_desn_Q20 <- svydesign(ids = ~0, weights = ~weights, data = rev2_wts,na.rm = TRUE) # New survey design
# Ordinal regression
Q20_mod1 <- svyolr(Q20_NewLicense ~ Gender + Income + Children + New.Old + COVID + Race_Ethn + SW_FW2 + Age, design=svy_desn_Q20)
summary(Q20_mod1) # examine model results
VIF(Q20_mod1) # collinearity
Q20_mod2 <- svyolr(Q20_NewLicense ~ New.Old + COVID + SW_FW2, design=svy_desn_Q20)
summary(Q20_mod2) # examine model results
VIF(Q20_mod2) # collinearity
# Survey response proportions
Q20_prop <- rev2_wts %>%
filter(!(is.na(Q20_NewLicense))) %>% # remove NA values from Gender & Age
count(Q20_NewLicense) %>%
mutate(freq = prop.table(n)) # get proportions
# Create summary table
Q20_coef_table <- coef(summary(Q20_mod2))
Q20_pval <- pnorm(abs(Q20_coef_table[, "t value"]),lower.tail = FALSE)* 2
Q20_AOR <- exp(coef(Q20_mod2))
Q20_summary_table <- cbind(Q20_coef_table, "p value" = round(Q20_pval,3),Odds_Ratio = Q20_AOR)# Table of responses
Q20_prop %>%
kbl(caption = "Q20 How likely are you to purchase another fishing license?",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | Q20_NewLicense | n | freq |
|---|---|---|
| Not likely | 18 | 0.01 |
| Somewhat likely | 139 | 0.07 |
| Very likely | 1798 | 0.92 |
# Histogram of responses
ggplot(model.frame(Q20_mod2), aes(Q20_NewLicense, fill = New.Old)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q20. How likely are you to buy another fishing license?") + facet_wrap(~SW_FW2, ncol = 2)
ggplot(model.frame(Q20_mod2), aes(Q20_NewLicense, fill = COVID)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q20. How likely are you to buy another fishing license?")
# Summary table
Q20_summary_table %>%
kbl(caption = "Q20 Model Odds Ratios",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | Value | Std. Error | t value | p value | Odds_Ratio | |
|---|---|---|---|---|---|
| New.OldOld | 1.19 | 0.31 | 3.86 | 0.00 | 3.27 |
| COVIDYes | -1.02 | 0.25 | -4.06 | 0.00 | 0.36 |
| SW_FW2SW | -0.43 | 0.29 | -1.50 | 0.13 | 0.65 |
| SW_FW2SW.FW | 0.63 | 0.29 | 2.19 | 0.03 | 1.88 |
| Not likely|Somewhat likely | -4.93 | 0.42 | -11.74 | 0.00 | 0.01 |
| Somewhat likely|Very likely | -2.58 | 0.23 | -11.44 | 0.00 | 0.08 |
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$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
# 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_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(model.frame(Q21_mod2), aes(Q21_AutoRenew, fill = Income)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q21. Likeliness to enroll in auto-renew program") 
ggplot(model.frame(Q21_mod2), aes(Q21_AutoRenew, fill = Race_Ethn)) + geom_histogram(stat = "count") + theme_classic() + 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 = "Q21. Likeliness to enroll in auto-renew program")
# Summary table
Q21_summary_table %>%
kbl(caption = "Q21 Model Odds Ratios",digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria") | Value | Std. Error | t value | p value | Odds_ratio | |
|---|---|---|---|---|---|
| IncomeBetween $20,000 and $39,999 | -0.75 | 0.26 | -2.85 | 0.00 | 0.47 |
| IncomeBetween $40,000 and $59,999 | -0.38 | 0.21 | -1.84 | 0.07 | 0.69 |
| IncomeBetween $60,000 and $79,999 | 0.00 | 0.18 | 0.02 | 0.98 | 1.00 |
| IncomeBetween $80,000 and $99,999 | -0.06 | 0.21 | -0.28 | 0.78 | 0.94 |
| IncomeLess than $20,000 | -0.69 | 0.39 | -1.77 | 0.08 | 0.50 |
| IncomeOver $150,000 | 0.47 | 0.18 | 2.60 | 0.01 | 1.60 |
| 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.72 | 0.48 | -11.79 | 0.00 | 0.00 |
| Unlikely|Not sure | -3.04 | 0.32 | -9.47 | 0.00 | 0.05 |
| Not sure|Likely | -1.17 | 0.30 | -3.96 | 0.00 | 0.31 |
| Likely|Extremely likely | 0.15 | 0.30 | 0.51 | 0.61 | 1.16 |
Most anglers (68%) reported that they were either likely (28%) or extremely likely (40%) to enroll in an auto-renew program, while only 7% were unlikely and 24% were unsure.
The odds of being likely to enroll in an autorenew program were: 1.02 times more likely for anglers with an income between $100-$150K compared to anglers with an income between $20K-$40K, 1.66 times more likely for anglers with an income over $150K compared to anglers with an income between $100-$150K, 1.96 times more likely for anglers that identified as an Other 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
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")
# 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_tbl <- svytable(~Age + Agreement, svy_desn_Q24)
Q24_Children_tbl <- svytable(~Children + Agreement, svy_desn_Q24)# 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. Who do you go fishing with?") 
# Chi-square test results
Q24_Age_tbl %>%
kbl(caption = "Q24 Who do you fish with?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | Not interested | Improve skills | What species | Baits and lures | Clean fish | Buy license | Available amenities | Understand Regulations | Where to fish | |
|---|---|---|---|---|---|---|---|---|---|
| 18-24 | 30 | 249 | 171 | 254 | 145 | 39 | 145 | 100 | 224 |
| 25-34 | 25 | 255 | 209 | 290 | 185 | 30 | 188 | 131 | 270 |
| 35-44 | 49 | 239 | 190 | 293 | 142 | 55 | 193 | 121 | 266 |
| 45-54 | 44 | 160 | 126 | 216 | 91 | 28 | 144 | 97 | 201 |
| 55-64 | 29 | 95 | 72 | 134 | 52 | 22 | 104 | 58 | 116 |
| 65 or older | 19 | 42 | 33 | 71 | 18 | 8 | 46 | 32 | 53 |
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-05corrplot(Q24_Age_chisq$residuals, is.cor = FALSE)
Q24_Children_tbl %>%
kbl(caption = "Q24 Who do you fish with?",digits = 0) %>%
kable_classic(full_width = F, html_font = "Cambria") | Not interested | Improve skills | What species | Baits and lures | Clean fish | Buy license | Available amenities | Understand Regulations | Where to fish | |
|---|---|---|---|---|---|---|---|---|---|
| 18 or over | 77 | 238 | 192 | 340 | 134 | 43 | 234 | 148 | 297 |
| No children | 68 | 455 | 338 | 490 | 282 | 69 | 297 | 205 | 459 |
| Under 17 | 19 | 160 | 112 | 182 | 116 | 34 | 130 | 90 | 173 |
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.0009995corrplot(Q24_Children_chisq$residuals, is.cor = FALSE)
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.