# Load packages library(weights) library(questionr) library(survey) library(Rmisc) library(tidyr) library(dplyr) library(data.table) library(ggplot2) library(sf) library(tigris) library(noncensus) # Read in data swp = read.csv("./swp_data.csv", header=T, stringsAsFactors=F) data(zip_codes) data(counties) coCounties = counties("CO", cb=T, resolution="5m", class="sf") # Process data ## Bin community size swp$commGroup = swp$commSize swp$commGroup[swp$commGroup == 1 | swp$commGroup == 2 | swp$commGroup == 3 | swp$commGroup == 4] = "City\n(25k-250k people)" swp$commGroup[swp$commGroup == 5 | swp$commGroup == 6 | swp$commGroup == 7] = "Town\n(5k-25k people)" swp$commGroup[swp$commGroup == 8] = "Farm or Rural Area\n(< 5k people)" ## Binarize children in household swp$children = ifelse(!is.na(swp$nChildren), ifelse(swp$nChildren < 1, "No Children", "Children"), NA) ## Extract dog or cat ownership swp$petOwn = rep(NA,nrow(swp)) for (i in 1:nrow(swp)) { if (grepl("1", swp$pets[i]) | grepl("2", swp$pets[i])) { swp$petOwn[i] = "Own\nDogs or Cats" } else { swp$petOwn[i] = "Do Not Own\nDogs or Cats" } } ## Binarize identification swp$idWildlife = swp$identifyWildAdv swp$idWildlife[swp$idWildlife == 3 | swp$idWildlife == 4] = 2 swp$idAnimalRights = swp$identifyAnimRightsAdv swp$idAnimalRights[swp$idAnimalRights == 3 | swp$idAnimalRights == 4] = 2 swp$idGun = swp$identifyGunRightsAdv swp$idGun[swp$idGun == 3 | swp$idGun == 4] = 2 swp$idProperty = swp$identifyPropRightsAdv swp$idProperty[swp$idProperty == 3 | swp$idProperty == 4] = 2 swp$idHunter = swp$identifyHunter swp$idHunter[swp$idHunter == 3 | swp$idHunter == 4] = 2 swp$idRancher = swp$identifyRancher swp$idRancher[swp$idRancher == 3 | swp$idRancher == 4] = 2 swp$idConservationist = swp$identifyConserv swp$idConservationist[swp$idConservationist == 3 | swp$idConservationist == 4] = 2 ## Bin perceived impact swp$impact = rep(NA,nrow(swp)) for (i in 1:nrow(swp)) { if (swp$impactExtent[i] < 4) { swp$impact[i] = "Negative" } else if (swp$impactExtent[i] == 4) { swp$impact[i] = "None" } else if (swp$impactExtent[i] > 4) { swp$impact[i] = "Positive" } } # Define survey design swpSvy = svydesign(id=~1, strata=~region, weight=~weight, data=swp) # Define map theme ## Original source: ## https://timogrossenbacher.ch/2016/12/beautiful-thematic-maps-with-ggplot2-only/ theme_map <- function(...) { theme_minimal() + theme( text=element_text(color="#22211d"), axis.line=element_blank(), axis.text.x=element_blank(), axis.text.y=element_blank(), axis.ticks=element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), panel.grid=element_line(color="transparent"), panel.background=element_rect(fill="transparent", color=NA), panel.border=element_blank(), panel.spacing=unit(c(-.1,0.2,.2,0.2), "cm"), plot.background=element_rect(fill="transparent", color=NA), plot.title=element_text(hjust=0.5, color="#4e4d47", face="bold"), plot.subtitle=element_text(hjust=0.5, color="#4e4d47", margin=margin(b=-0.1, t=-0.1, l=2, unit="cm"), debug=F), plot.margin=unit(c(.5,.5,.2,.5), "cm"), plot.caption=element_text(size=6, hjust=0.92, margin=margin(t=0.2, b=0, unit = "cm"), color = "#939184"), ... ) } # Add region for counties eastPlainsCounties = c("Baca","Bent","Cheyenne","Crowley","Elbert","Huerfano", "Kiowa","Kit Carson","Las Animas","Lincoln","Logan", "Morgan","Otero","Phillips","Prowers","Sedgwick", "Washington","Yuma") westSlopeCounties = c("Alamosa","Archuleta","Chaffee","Clear Creek","Conejos", "Costilla","Custer","Delta","Dolores","Eagle","Fremont", "Garfield","Gilpin","Grand","Gunnison","Hinsdale", "Jackson","La Plata","Lake","Mesa","Mineral","Moffat", "Montezuma","Montrose","Ouray","Park","Pitkin", "Rio Blanco","Rio Grande","Routt","Saguache","San Juan", "San Miguel","Summit","Teller") frontRangeCounties = c("Adams","Arapahoe","Boulder","Denver","Douglas", "El Paso","Jefferson","Larimer","Pueblo","Weld", "Broomfield") coCounties$region = NA for (i in 1:nrow(coCounties)) { if (coCounties$NAME[i] %in% eastPlainsCounties) { coCounties$region[i] = "Eastern Plains" } else { if (coCounties$NAME[i] %in% westSlopeCounties) { coCounties$region[i] = "Western Slope" } else { if (coCounties$NAME[i] %in% frontRangeCounties) { coCounties$region[i] = "Front Range" } } } } # Aggregate region shapes coRegions = coCounties %>% group_by(region) %>% dplyr::summarize() # Summarize sample by region swpRegions = swp %>% group_by(region) %>% dplyr::summarize(regionN = n(), regionSupport = sum(votingIntention == 1) / length(votingIntention)) %>% as.data.frame() coCounties = left_join(coCounties, swpRegions, by=c("region")) coRegions = left_join(coRegions, swpRegions, by=c("region")) # Set region labels coRegions$label = paste0(coRegions$region, "\n(N = ", coRegions$regionN, ")") frontRangeLabel = coCounties[coCounties$NAME == "Larimer" | coCounties$NAME == "Weld" | coCounties$NAME == "Boulder" | coCounties$NAME == "Broomfield",] %>% group_by(region) %>% dplyr::summarize() frontRangeLabel$regionN = coRegions$regionN[coRegions$region == "Front Range"] frontRangeLabel$label = coRegions$label[coRegions$region == "Front Range"] easternPlainsLabel = coCounties[coCounties$NAME == "Lincoln" | coCounties$NAME == "Cheyenne" | coCounties$NAME == "Crowley" | coCounties$NAME == "Kiowa" | coCounties$NAME == "Otero",] %>% group_by(region) %>% dplyr::summarize() easternPlainsLabel$regionN = coRegions$regionN[coRegions$region == "Eastern Plains"] easternPlainsLabel$label = coRegions$label[coRegions$region == "Eastern Plains"] # Reorder region labels coCounties$region = factor(coCounties$region, levels=c("Western Slope","Front Range","Eastern Plains")) coRegions$region = factor(coRegions$region, levels=c("Western Slope","Front Range","Eastern Plains")) coRegions$label = factor(coRegions$label, levels=c("Western Slope\n(N = 277)", "Front Range\n(N = 365)", "Eastern Plains\n(N = 92)")) # FIGURE 1: State map coCountiesPlot = ggplot(coCounties) + geom_sf(aes(geometry=geometry), color=NA) + geom_sf(aes(geometry=geometry, fill=region), coRegions, color=NA) + geom_sf(aes(geometry=geometry), fill=NA, color="#555555", size=0.1) + geom_sf(aes(geometry=geometry), coRegions, fill=NA, color="#000000", size=0.5) + geom_sf_text(aes(label=NAME), size=2.5) + geom_sf_text(aes(label=label), coRegions[coRegions$region == "Western Slope",], size=4) + geom_sf_text(aes(label=label), frontRangeLabel, size=4) + geom_sf_text(aes(label=label), easternPlainsLabel, size=4) + scale_fill_manual(name="Region", values=c("#BBCCEE","#EEEEBB","#FFCCCC")) + theme_map() + theme(legend.position="none", plot.margin=unit(rep(0,4), "in")) ggsave("./paper_Fig1.pdf", coCountiesPlot, "pdf", width=8.5, height=6.5, units="in") ggsave("./paper_Fig1.jpg", coCountiesPlot, "jpg", width=8.5, height=6.5, units="in") ggsave("./paper_Fig1.png", coCountiesPlot, "png", width=8.5, height=6.5, units="in", dpi=600) # Summarize overall voting intention ### Unweighted # voteOverall = summarySE(swp, "votingIntention") # voteOverall = rbind(voteOverall, voteOverall) # voteOverall[2,3] = 1 - voteOverall[1,3] # voteOverall[,1] = c("Yes","No") # voteOverall$.id = factor(voteOverall$.id, levels=voteOverall$.id) # voteOverall$type = "Overall" ### Weighted voteOverall = as.data.frame(svymean(~factor(votingIntention), swpSvy)) voteOverall = cbind(.id=rownames(voteOverall), data.frame(voteOverall, row.names=NULL)) voteOverall = cbind(voteOverall, data.frame(confint(svymean(~factor(votingIntention), swpSvy)), row.names=NULL)) voteOverall = voteOverall[c(2,1),] colnames(voteOverall) = c(".id","votingIntention","se","cilower","ciupper") voteOverall[,1] = c("Yes","No") voteOverall$.id = factor(voteOverall$.id, levels=voteOverall$.id) voteOverall$type = "Overall" # FIGURE 2: Overall voting intention voteOverallPlot = ggplot(aes(x=.id, y=votingIntention), data=voteOverall) + geom_bar(stat="identity") + geom_errorbar(aes(ymin=cilower, ymax=ciupper), width=0.4) + geom_text(aes(y=0.975, label=paste0(formatC(round(voteOverall$votingIntention,3)*100, format="f", digits=1),"%"))) + scale_x_discrete(name="Voting Intention") + scale_y_continuous(name="Weighted Proportion of Respondents", limits=c(0,1), breaks=c(seq(0,1,0.25)), expand=c(0,0)) + theme_minimal() + theme(panel.grid.major.x=element_blank(), axis.text.x=element_text(margin=margin(b=10))) ggsave("./paper_Fig2.pdf", voteOverallPlot, "pdf", width=4, height=4, units="in") ggsave("./paper_Fig2.jpg", voteOverallPlot, "jpg", width=4, height=4, units="in") ggsave("./paper_Fig2.png", voteOverallPlot, "png", width=4, height=4, units="in", dpi=600) # Summarize voting intention by demographics ## By region ### Unweighted (Weighted N/A) sweep(table(swp$votingIntention, swp$region, useNA="no"),2, colSums(table(swp$votingIntention, swp$region, useNA="no")),'/') voteRegion = summarySE(swp, "votingIntention", "region") voteRegion = voteRegion[,-c(2,4)] voteRegion$ciupper = voteRegion$votingIntention + voteRegion$ci voteRegion$ci = voteRegion$votingIntention - voteRegion$ci colnames(voteRegion) = c(".id","votingIntention","se","cilower","ciupper") voteRegion[,1] = c("Eastern Plains","Front Range","Western Slope") voteRegion = voteRegion[c(3,2,1),] voteRegion$.id = factor(voteRegion$.id, levels=voteRegion$.id) voteRegion$type = "Region" ## By community size (binned) ### Unweighted # voteComm = summarySE(swp, "votingIntention", "commGroup") # colnames(voteComm)[1] = ".id" # voteComm = voteComm[c(2,3,1),] # voteComm$.id = factor(voteComm$.id, levels=voteComm$.id) # voteComm$type = "Community Size" ### Weighted voteComm = as.data.frame(svyby(~votingIntention, ~factor(commGroup), swpSvy, svymean, vartype=c("se","ci"))) colnames(voteComm) = c(".id","votingIntention","se","cilower","ciupper") voteComm = voteComm[c(2,3,1),] voteComm$.id = factor(voteComm$.id, levels=voteComm$.id) voteComm$type = "Community Size" ## By children in household ### Unweighted # voteChildren = summarySE(swp, "votingIntention", "children") # voteChildren = voteChildren[-3,] # colnames(voteChildren)[1] = ".id" # voteChildren$.id = factor(voteChildren$.id, levels=voteChildren$.id) # voteChildren$type = "Children in Household" ### Weighted voteChildren = as.data.frame(svyby(~votingIntention, ~factor(children), swpSvy, svymean, vartype=c("se","ci"))) colnames(voteChildren) = c(".id","votingIntention","se","cilower","ciupper") voteChildren$.id = factor(voteChildren$.id, levels=voteChildren$.id) voteChildren$type = "Children in Household" ## By dog and cat ownership ### Unweighted # votePetOwn = summarySE(swp, "votingIntention", "petOwn") # colnames(votePetOwn)[1] = ".id" # votePetOwn = votePetOwn[c(2,1),] # votePetOwn$.id = factor(votePetOwn$.id, levels=votePetOwn$.id) # votePetOwn$type = "Pet Ownership" ### Weighted votePetOwn = as.data.frame(svyby(~votingIntention, ~factor(petOwn), swpSvy, svymean, vartype=c("se","ci"))) colnames(votePetOwn) = c(".id","votingIntention","se","cilower","ciupper") votePetOwn = votePetOwn[c(2,1),] votePetOwn$.id = factor(votePetOwn$.id, levels=votePetOwn$.id) votePetOwn$type = "Pet Ownership" ## By gender ### Unweighted # voteGender = summarySE(swp, "votingIntention", "gender") # voteGender = voteGender[-4,] # colnames(voteGender)[1] = ".id" # voteGender[,1] = c("Male","Female","Non-Binary") # voteGender$.id = factor(voteGender$.id, levels=voteGender$.id) # voteGender$type = "Gender" ### Weighted voteGender = as.data.frame(svyby(~votingIntention, ~factor(gender), swpSvy, svymean, vartype=c("se","ci"))) colnames(voteGender) = c(".id","votingIntention","se","cilower","ciupper") voteGender[,1] = c("Male","Female","Non-Binary") voteGender$.id = factor(voteGender$.id, levels=voteGender$.id) voteGender$type = "Gender" ## By age group ### Unweighted # voteAge = summarySE(swp, "votingIntention", "ageGroup") # colnames(voteAge)[1] = ".id" # voteAge$.id = factor(voteAge$.id, levels=voteAge$.id) # voteAge$type = "Age Group" ### Weighted voteAge = as.data.frame(svyby(~votingIntention, ~factor(ageGroup), swpSvy, svymean, vartype=c("se","ci"))) colnames(voteAge) = c(".id","votingIntention","se","cilower","ciupper") voteAge$.id = factor(voteAge$.id, levels=voteAge$.id) voteAge$type = "Age Group" ## By income ### Unweighted # voteIncome = summarySE(swp, "votingIntention", "income") # voteIncome = voteIncome[-7,] # colnames(voteIncome)[1] = ".id" # voteIncome[,1] = c("< 10k","10k-25k","25k-50k","50k-100k","100k-250k","> 250k") # voteIncome$.id = factor(voteIncome$.id, levels=voteIncome$.id) # voteIncome$type = "Income" ### Weighted voteIncome = as.data.frame(svyby(~votingIntention, ~factor(income), swpSvy, svymean, vartype=c("se","ci"))) colnames(voteIncome) = c(".id","votingIntention","se","cilower","ciupper") voteIncome[,1] = c("< 10k","10k-25k","25k-50k","50k-100k","100k-250k","> 250k") voteIncome$.id = factor(voteIncome$.id, levels=voteIncome$.id) voteIncome$type = "Income" ## By education ### Unweighted # voteEducation = summarySE(swp, "votingIntention", "education") # colnames(voteEducation)[1] = ".id" # voteEducation[,1] = c("Less than\nHigh School","High School","Associate's", # "Bachelor's","Graduate/\nProfessional") # voteEducation$.id = factor(voteEducation$.id, levels=voteEducation$.id) # voteEducation$type = "Education" ### Weighted voteEducation = as.data.frame(svyby(~votingIntention, ~factor(education), swpSvy, svymean, vartype=c("se","ci"))) colnames(voteEducation) = c(".id","votingIntention","se","cilower","ciupper") voteEducation[,1] = c("Less than\nHigh School","High School","Associate's", "Bachelor's","Graduate/\nProfessional") voteEducation$.id = factor(voteEducation$.id, levels=voteEducation$.id) voteEducation$type = "Education" # Combine voting intention by demographics voteCombined = rbindlist(list(voteRegion, voteComm, voteChildren, votePetOwn, voteGender[voteGender$.id != "Non-Binary",], voteAge, voteIncome, voteEducation), use.names=F) voteCombined$type = factor(voteCombined$type, levels=unique(voteCombined$type)) # FIGURE 3: Voting intention by demographics voteCombinedPlot = ggplot(aes(x=.id, y=votingIntention), data=voteCombined) + geom_bar(stat="identity") + geom_errorbar(aes(ymin=cilower, ymax=ifelse(ciupper <= 1, ciupper, 1), width=0.4)) + geom_text(aes(y=0.975, label=paste0(formatC(round(voteCombined$votingIntention,3)*100, format="f", digits=1),"%")), size=2.5) + facet_wrap(facets=vars(type), nrow=2, drop=T, scales="free") + scale_x_discrete(name="Demographic Group") + scale_y_continuous(name="Proportion that Would Vote in Favor of Wolf Reintroduction", limits=c(0,1), breaks=c(seq(0,1,0.25)), expand=c(0,0)) + theme_minimal() + theme(panel.grid.major.x=element_blank(), axis.text.x=element_text(angle=45, hjust=1, vjust=0.9, margin=margin(b=10)), strip.background=element_rect(fill="grey85", color=NA)) ggsave("./paper_Fig3.pdf", voteCombinedPlot, "pdf", width=10, height=8, units="in") ggsave("./paper_Fig3.jpg", voteCombinedPlot, "jpg", width=10, height=8, units="in") ggsave("./paper_Fig3.png", voteCombinedPlot, "png", width=10, height=8, units="in", dpi=600) # Summarize voting intention by identity ## Wildlife advocate ### Unweighted # idWildlife = summarySE(swp, "votingIntention", "idWildlife") # colnames(idWildlife)[1] = ".id" # idWildlife[,1] = c("Do Not\nIdentify As","Identify As") # idWildlife = idWildlife[c(2,1),] # idWildlife$.id = factor(idWildlife$.id, levels=idWildlife$.id) # idWildlife$type = "Wildlife\nAdvocate" ### Weighted idWildlife = as.data.frame(svyby(~votingIntention, ~factor(idWildlife), swpSvy, svymean, vartype=c("se","ci"))) colnames(idWildlife) = c(".id","votingIntention","se","cilower","ciupper") idWildlife[,1] = c("Do Not\nIdentify As","Identify As") idWildlife = idWildlife[c(2,1),] idWildlife$.id = factor(idWildlife$.id, levels=idWildlife$.id) idWildlife$type = "Wildlife\nAdvocate" ## Animal rights advocate ### Unweighted # idAnimalRights = summarySE(swp, "votingIntention", "idAnimalRights") # colnames(idAnimalRights)[1] = ".id" # idAnimalRights[,1] = c("Do Not\nIdentify As","Identify As") # idAnimalRights = idAnimalRights[c(2,1),] # idAnimalRights$.id = factor(idAnimalRights$.id, levels=idAnimalRights$.id) # idAnimalRights$type = "Animal Rights\nAdvocate" ### Weighted idAnimalRights = as.data.frame(svyby(~votingIntention, ~factor(idAnimalRights), swpSvy, svymean, vartype=c("se","ci"))) colnames(idAnimalRights) = c(".id","votingIntention","se","cilower","ciupper") idAnimalRights[,1] = c("Do Not\nIdentify As","Identify As") idAnimalRights = idAnimalRights[c(2,1),] idAnimalRights$.id = factor(idAnimalRights$.id, levels=idAnimalRights$.id) idAnimalRights$type = "Animal Rights\nAdvocate" ## Gun rights advocate ### Unweighted # idGun = summarySE(swp, "votingIntention", "idGun") # colnames(idGun)[1] = ".id" # idGun[,1] = c("Do Not\nIdentify As","Identify As") # idGun = idGun[c(2,1),] # idGun$.id = factor(idGun$.id, levels=idGun$.id) # idGun$type = "Gun Rights\nAdvocate" ### Weighted idGun = as.data.frame(svyby(~votingIntention, ~factor(idGun), swpSvy, svymean, vartype=c("se","ci"))) colnames(idGun) = c(".id","votingIntention","se","cilower","ciupper") idGun[,1] = c("Do Not\nIdentify As","Identify As") idGun = idGun[c(2,1),] idGun$.id = factor(idGun$.id, levels=idGun$.id) idGun$type = "Gun Rights\nAdvocate" ## Property rights advocate ### Unweighted # idProperty = summarySE(swp, "votingIntention", "idProperty") # colnames(idProperty)[1] = ".id" # idProperty[,1] = c("Do Not\nIdentify As","Identify As") # idProperty = idProperty[c(2,1),] # idProperty$.id = factor(idProperty$.id, levels=idProperty$.id) # idProperty$type = "Property Rights\nAdvocate" ### Weighted idProperty = as.data.frame(svyby(~votingIntention, ~factor(idProperty), swpSvy, svymean, vartype=c("se","ci"))) colnames(idProperty) = c(".id","votingIntention","se","cilower","ciupper") idProperty[,1] = c("Do Not\nIdentify As","Identify As") idProperty = idProperty[c(2,1),] idProperty$.id = factor(idProperty$.id, levels=idProperty$.id) idProperty$type = "Property Rights\nAdvocate" ## Hunter ### Unweighted # idHunter = summarySE(swp, "votingIntention", "idHunter") # colnames(idHunter)[1] = ".id" # idHunter[,1] = c("Do Not\nIdentify As","Identify As") # idHunter = idHunter[c(2,1),] # idHunter$.id = factor(idHunter$.id, levels=idHunter$.id) # idHunter$type = "Hunter" ### Weighted idHunter = as.data.frame(svyby(~votingIntention, ~factor(idHunter), swpSvy, svymean, vartype=c("se","ci"))) colnames(idHunter) = c(".id","votingIntention","se","cilower","ciupper") idHunter[,1] = c("Do Not\nIdentify As","Identify As") idHunter = idHunter[c(2,1),] idHunter$.id = factor(idHunter$.id, levels=idHunter$.id) idHunter$type = "Hunter" ## Rancher ### Unweighted # idRancher = summarySE(swp, "votingIntention", "idRancher") # colnames(idRancher)[1] = ".id" # idRancher[,1] = c("Do Not\nIdentify As","Identify As") # idRancher = idRancher[c(2,1),] # idRancher$.id = factor(idRancher$.id, levels=idRancher$.id) # idRancher$type = "Rancher" ### Weighted idRancher = as.data.frame(svyby(~votingIntention, ~factor(idRancher), swpSvy, svymean, vartype=c("se","ci"))) colnames(idRancher) = c(".id","votingIntention","se","cilower","ciupper") idRancher[,1] = c("Do Not\nIdentify As","Identify As") idRancher = idRancher[c(2,1),] idRancher$.id = factor(idRancher$.id, levels=idRancher$.id) idRancher$type = "Rancher" ## Rancher ### Unweighted # idConservationist = summarySE(swp, "votingIntention", "idConservationist") # colnames(idConservationist)[1] = ".id" # idConservationist[,1] = c("Do Not\nIdentify As","Identify As") # idConservationist = idConservationist[c(2,1),] # idConservationist$.id = factor(idConservationist$.id, levels=idConservationist$.id) # idConservationist$type = "Conservationist" ### Weighted idConservationist = as.data.frame(svyby(~votingIntention, ~factor(idConservationist), swpSvy, svymean, vartype=c("se","ci"))) colnames(idConservationist) = c(".id","votingIntention","se","cilower","ciupper") idConservationist[,1] = c("Do Not\nIdentify As","Identify As") idConservationist = idConservationist[c(2,1),] idConservationist$.id = factor(idConservationist$.id, levels=idConservationist$.id) idConservationist$type = "Conservationist" # Combine voting intention by identity idCombined = rbindlist(list(idWildlife, idAnimalRights, idGun, idProperty, idHunter, idRancher, idConservationist), use.names=F) idCombined$type = factor(idCombined$type, levels=unique(idCombined$type)) # FIGURE 4: Voting intention by identity idCombinedPlot = ggplot(aes(x=.id, y=votingIntention), data=idCombined) + geom_bar(stat="identity") + geom_errorbar(aes(ymin=cilower, ymax=ciupper), width=0.4) + geom_text(aes(y=0.975, label=paste0(formatC(round(idCombined$votingIntention,3)*100, format="f", digits=1),"%")), size=3.25) + facet_grid(cols=vars(type), drop=T, scales="free", space="free") + scale_x_discrete(name="Self Identification") + scale_y_continuous(name="Proportion that Would Vote in Favor of Wolf Reintroduction", limits=c(0,1), breaks=c(seq(0,1,0.25)), expand=c(0,0)) + theme_minimal() + theme(panel.grid.major.x=element_blank(), axis.text.x=element_text(angle=45, hjust=1, vjust=0.9, margin=margin(b=10)), strip.background=element_rect(fill="grey85", color=NA)) ggsave("./paper_Fig4.pdf", idCombinedPlot, "pdf", width=8, height=6, units="in") ggsave("./paper_Fig4.jpg", idCombinedPlot, "jpg", width=8, height=6, units="in") ggsave("./paper_Fig4.png", idCombinedPlot, "png", width=8, height=6, units="in", dpi=600) # Summarize support for management options ## Limit for deer pop decline ### Unweighted # table(factor(ifelse(swp$mgmtLimitDeerElk > 4, 1, 0))) # prop.table(table(factor(ifelse(swp$mgmtLimitDeerElk > 4, 1, 0)))) ### Weighted mgmtDeerPop = as.data.frame(svymean(~factor(ifelse(swp$mgmtLimitDeerElk > 4, 1, 0)), swpSvy)) mgmtDeerPop = cbind(.id=rownames(mgmtDeerPop), data.frame(mgmtDeerPop, row.names=NULL)) mgmtDeerPop = cbind(mgmtDeerPop, data.frame(confint(svymean(~factor(ifelse(mgmtLimitDeerElk > 4, 1, 0)), swpSvy)), row.names=NULL)) colnames(mgmtDeerPop) = c(".id","prop","se","cilower","ciupper") mgmtDeerPop = mgmtDeerPop[c(2,1),] mgmtDeerPop[,1] = c("Yes","No") mgmtDeerPop$.id = factor(mgmtDeerPop$.id, levels=mgmtDeerPop$.id) mgmtDeerPop$type = "Limit for Deer and Elk\nPopulation Decline" ## Lethal removal for livestock loss ### Unweighted # table(factor(ifelse(swp$mgmtLethalRemoval > 4, 1, 0))) # prop.table(table(factor(ifelse(swp$mgmtLethalRemoval > 4, 1, 0)))) ### Weighted mgmtLethalLvstck = as.data.frame(svymean(~factor(ifelse(swp$mgmtLethalRemoval > 4, 1, 0)), swpSvy)) mgmtLethalLvstck = cbind(.id=rownames(mgmtLethalLvstck), data.frame(mgmtLethalLvstck, row.names=NULL)) mgmtLethalLvstck = cbind(mgmtLethalLvstck, data.frame(confint(svymean(~factor(ifelse(swp$mgmtLethalRemoval > 4, 1, 0)), swpSvy)), row.names=NULL)) colnames(mgmtLethalLvstck) = c(".id","prop","se","cilower","ciupper") mgmtLethalLvstck = mgmtLethalLvstck[c(2,1),] mgmtLethalLvstck[,1] = c("Yes","No") mgmtLethalLvstck$.id = factor(mgmtLethalLvstck$.id, levels=mgmtLethalLvstck$.id) mgmtLethalLvstck$type = "Lethal Removal for\nLivestock Loss" ## Compensate for livestock loss ### Unweighted # table(factor(ifelse(swp$mgmtCompLivestock > 4, 1, 0))) # prop.table(table(factor(ifelse(swp$mgmtCompLivestock > 4, 1, 0)))) ### Weighted mgmtCompLvstck = as.data.frame(svymean(~factor(ifelse(swp$mgmtCompLivestock > 4, 1, 0)), swpSvy)) mgmtCompLvstck = cbind(.id=rownames(mgmtCompLvstck), data.frame(mgmtCompLvstck, row.names=NULL)) mgmtCompLvstck = cbind(mgmtCompLvstck, data.frame(confint(svymean(~factor(ifelse(swp$mgmtCompLivestock > 4, 1, 0)), swpSvy)), row.names=NULL)) colnames(mgmtCompLvstck) = c(".id","prop","se","cilower","ciupper") mgmtCompLvstck = mgmtCompLvstck[c(2,1),] mgmtCompLvstck[,1] = c("Yes","No") mgmtCompLvstck$.id = factor(mgmtCompLvstck$.id, levels=mgmtCompLvstck$.id) mgmtCompLvstck$type = "Compensate for\nLivestock Loss" ## Compensate for livestock loss with license rev ### Unweighted # table(factor(ifelse(swp$mgmtCompLivestockLicense > 4, 1, 0))) # prop.table(table(factor(ifelse(swp$mgmtCompLivestockLicense > 4, 1, 0)))) ### Weighted mgmtCompLvstckLicense = as.data.frame(svymean(~factor(ifelse(swp$mgmtCompLivestockLicense > 4, 1, 0)), swpSvy)) mgmtCompLvstckLicense = cbind(.id=rownames(mgmtCompLvstckLicense), data.frame(mgmtCompLvstckLicense, row.names=NULL)) mgmtCompLvstckLicense = cbind(mgmtCompLvstckLicense, data.frame(confint(svymean(~factor(ifelse(swp$mgmtCompLivestockLicense > 4, 1, 0)), swpSvy)), row.names=NULL)) colnames(mgmtCompLvstckLicense) = c(".id","prop","se","cilower","ciupper") mgmtCompLvstckLicense = mgmtCompLvstckLicense[c(2,1),] mgmtCompLvstckLicense[,1] = c("Yes","No") mgmtCompLvstckLicense$.id = factor(mgmtCompLvstckLicense$.id, levels=mgmtCompLvstckLicense$.id) mgmtCompLvstckLicense$type = "Compensate for\nLivestock Loss with\nLicensing Revenue" ## Compensate for livestock loss with tax rev ### Unweighted # table(factor(ifelse(swp$mgmtCompLivestockTax > 4, 1, 0))) # prop.table(table(factor(ifelse(swp$mgmtCompLivestockTax > 4, 1, 0)))) ### Weighted mgmtCompLvstckTax = as.data.frame(svymean(~factor(ifelse(swp$mgmtCompLivestockTax > 4, 1, 0)), swpSvy)) mgmtCompLvstckTax = cbind(.id=rownames(mgmtCompLvstckTax), data.frame(mgmtCompLvstckTax, row.names=NULL)) mgmtCompLvstckTax = cbind(mgmtCompLvstckTax, data.frame(confint(svymean(~factor(ifelse(swp$mgmtCompLivestockTax > 4, 1, 0)), swpSvy)), row.names=NULL)) colnames(mgmtCompLvstckTax) = c(".id","prop","se","cilower","ciupper") mgmtCompLvstckTax = mgmtCompLvstckTax[c(2,1),] mgmtCompLvstckTax[,1] = c("Yes","No") mgmtCompLvstckTax$.id = factor(mgmtCompLvstckTax$.id, levels=mgmtCompLvstckTax$.id) mgmtCompLvstckTax$type = "Compensate for\nLivestock Loss with\nTax Revenue" ## Rec hunting ### Unweighted # table(factor(ifelse(swp$mgmtRecHunt > 4, 1, 0))) # prop.table(table(factor(ifelse(swp$mgmtRecHunt > 4, 1, 0)))) ### Weighted mgmtRecHunt = as.data.frame(svymean(~factor(ifelse(swp$mgmtRecHunt > 4, 1, 0)), swpSvy)) mgmtRecHunt = cbind(.id=rownames(mgmtRecHunt), data.frame(mgmtRecHunt, row.names=NULL)) mgmtRecHunt = cbind(mgmtRecHunt, data.frame(confint(svymean(~factor(ifelse(swp$mgmtRecHunt > 4, 1, 0)), swpSvy)), row.names=NULL)) colnames(mgmtRecHunt) = c(".id","prop","se","cilower","ciupper") mgmtRecHunt = mgmtRecHunt[c(2,1),] mgmtRecHunt[,1] = c("Yes","No") mgmtRecHunt$.id = factor(mgmtRecHunt$.id, levels=mgmtRecHunt$.id) mgmtRecHunt$type = "Recreational\nHunting" # Combine support for management options mgmtCombined = rbindlist(list(mgmtDeerPop, mgmtLethalLvstck, mgmtCompLvstck, mgmtCompLvstckLicense, mgmtCompLvstckTax, mgmtRecHunt), use.names=F) mgmtCombined = mgmtCombined[mgmtCombined$.id != "No",-1] mgmtCombined$type = factor(mgmtCombined$type, levels=unique(mgmtCombined$type)) # FIGURE 5: Support for management options mgmtCombinedPlot = ggplot(aes(x=reorder(type, -prop), y=prop), data=mgmtCombined) + geom_bar(stat="identity") + geom_errorbar(aes(ymin=cilower, ymax=ciupper), width=0.4) + geom_text(aes(y=0.975, label=paste0(formatC(round(mgmtCombined$prop,3)*100, format="f", digits=1),"%"))) + scale_x_discrete(name="Wolf Management Options") + scale_y_continuous(name="Weighted Proportion of Respondents in Support", limits=c(0,1), breaks=c(seq(0,1,0.25)), expand=c(0,0)) + theme_minimal() + theme(panel.grid.major.x=element_blank(), axis.text.x=element_text(angle=45, hjust=1, vjust=0.9, margin=margin(b=10)), strip.background=element_rect(fill="grey85", color=NA)) ggsave("./paper_Fig5.pdf", mgmtCombinedPlot, "pdf", width=8, height=6, units="in") ggsave("./paper_Fig5.jpg", mgmtCombinedPlot, "jpg", width=8, height=6, units="in") ggsave("./paper_Fig5.png", mgmtCombinedPlot, "png", width=8, height=6, units="in", dpi=600) # Summarize perceived impacts by region ## Unweighted (Weighted N/A) sweep(table(swp$impact, swp$region, useNA="no"),2, colSums(table(swp$impact, swp$region, useNA="no")),'/') impactRegion = as.data.frame(svyby(~factor(impact), ~factor(region), swpSvy, svymean, vartype=c("se","ci"))) impactRegion = data.frame(region=rep(c("Eastern Plains","Front Range","Western Slope"), each=3), impact=rep(c("Negative\nImpact","No\nImpact","Positive\nImpact"),3), mean=as.numeric(c(impactRegion[1,2:4], impactRegion[2,2:4], impactRegion[3,2:4])), se=as.numeric(c(impactRegion[1,5:7], impactRegion[2,5:7], impactRegion[3,5:7])), cilower=as.numeric(c(impactRegion[1,8:10], impactRegion[2,8:10], impactRegion[3,8:10])), ciupper=as.numeric(c(impactRegion[1,11:13], impactRegion[2,11:13], impactRegion[3,11:13]))) impactRegion = impactRegion[c(7:9,4:6,1:3),] impactRegion$region = factor(impactRegion$region, levels=unique(impactRegion$region)) impactRegion$impact = factor(impactRegion$impact, levels=unique(impactRegion$impact)) impactRegion$type = "Impact" # FIGURE 6: Perceived impacts by region impactRegionPlot = ggplot(aes(x=impact, y=mean), data=impactRegion) + geom_bar(stat="identity") + geom_errorbar(aes(ymin=cilower, ymax=ciupper), width=0.4) + geom_text(aes(y=0.725, label=paste0(formatC(round(impactRegion$mean,3)*100, format="f", digits=1),"%"))) + facet_grid(cols=vars(region), drop=T, scales="free", space="free") + scale_x_discrete(name="Perceived Personal Impact of Wolf Reintroduction") + scale_y_continuous(name="Proportion of Respondents", limits=c(0,0.75), breaks=c(seq(0,0.75,0.25)), expand=c(0,0)) + theme_minimal() + theme(panel.grid.major.x=element_blank(), axis.text.x=element_text(angle=45, hjust=1, vjust=0.9, margin=margin(b=10)), strip.background=element_rect(fill="grey85", color=NA)) ggsave("./paper_Fig6.pdf", impactRegionPlot, "pdf", width=8, height=6, units="in") ggsave("./paper_Fig6.jpg", impactRegionPlot, "jpg", width=8, height=6, units="in") ggsave("./paper_Fig6.png", impactRegionPlot, "png", width=8, height=6, units="in", dpi=600) # Summarize voting intention among hunters and ranchers ## Hunter ### Unweighted # identifyHunter = summarySE(swp, "votingIntention", "identifyHunter") # colnames(identifyHunter)[1] = ".id" # identifyHunter[,1] = c("Do Not Identify\nAt All","Identify\nSlight Amount", # "Identify\nModerate Amount","Identify\nGreat Deal") # identifyHunter$.id = factor(identifyHunter$.id, levels=identifyHunter$.id) # identifyHunter$type = "Hunter" ### Weighted identifyHunter = as.data.frame(svyby(~votingIntention, ~factor(identifyHunter), swpSvy, svymean, vartype=c("se","ci"))) colnames(identifyHunter) = c(".id","votingIntention","se","cilower","ciupper") identifyHunter[,1] = c("Do Not Identify\nAt All","Identify\nSlight Amount", "Identify\nModerate Amount","Identify\nGreat Deal") identifyHunter$.id = factor(identifyHunter$.id, levels=identifyHunter$.id) identifyHunter$type = "Hunter" ## Rancher ### Unweighted # identifyRancher = summarySE(swp, "votingIntention", "identifyRancher") # colnames(identifyRancher)[1] = ".id" # identifyRancher[,1] = c("Do Not Identify\nAt All","Identify\nSlight Amount", # "Identify\nModerate Amount","Identify\nGreat Deal") # identifyRancher$.id = factor(identifyRancher$.id, levels=identifyRancher$.id) # identifyRancher$type = "Rancher" ### Weighted identifyRancher = as.data.frame(svyby(~votingIntention, ~factor(identifyRancher), swpSvy, svymean, vartype=c("se","ci"))) colnames(identifyRancher) = c(".id","votingIntention","se","cilower","ciupper") identifyRancher[,1] = c("Do Not Identify\nAt All","Identify\nSlight Amount", "Identify\nModerate Amount","Identify\nGreat Deal") identifyRancher$.id = factor(identifyRancher$.id, levels=identifyRancher$.id) identifyRancher$type = "Rancher" # Combine voting intention by identity identifyHR = rbindlist(list(identifyHunter, identifyRancher), use.names=F) identifyHR$type = factor(identifyHR$type, levels=unique(identifyHR$type)) # FIGURE S-1: Voting intention by identity among hunters and ranchers identifyHRPlot = ggplot(aes(x=.id, y=votingIntention), data=identifyHR) + geom_bar(stat="identity") + geom_errorbar(aes(ymin=cilower, ymax=ciupper), width=0.4) + geom_text(aes(y=0.975, label=paste0(formatC(round(identifyHR$votingIntention,3)*100, format="f", digits=1),"%")), size=3.25) + facet_grid(cols=vars(type), drop=T, scales="free", space="free") + scale_x_discrete(name="Self Identification") + scale_y_continuous(name="Proportion that Would Vote in Favor of Wolf Reintroduction", limits=c(0,1), breaks=c(seq(0,1,0.25)), expand=c(0,0)) + theme_minimal() + theme(panel.grid.major.x=element_blank(), axis.text.x=element_text(angle=45, hjust=1, vjust=0.9, margin=margin(b=10)), strip.background=element_rect(fill="grey85", color=NA)) ggsave("./paper_FigS-1.pdf", identifyHRPlot, "pdf", width=8, height=6, units="in") ggsave("./paper_FigS-1.jpg", identifyHRPlot, "jpg", width=8, height=6, units="in") ggsave("./paper_FigS-1.png", identifyHRPlot, "png", width=8, height=6, units="in", dpi=600)