# Ashley Darst # MnDOT Pollinator Paper # December 18, 2023 # Load libraries library(tidyverse) library(data.table) library(vegan) library(expss) library(car) library(emmeans) library(glmmTMB) library(lme4) library(nlme) library(numDeriv) library(performance) library(stringi) library(anytime) #make a fn for replacing values in many cols f_dowle3natoblank = function(DT, x) { # or by number (slightly faster than by name) : for (j in x) set(DT,which(is.na(DT[[j]])),j,"") } #na to zero f_dowle3natozeros = function(DT, x) { # or by number (slightly faster than by name) : for (j in x) set(DT,which(is.na(DT[[j]])),j,"0") } # read in datasets -------------------------------------------------------- #read in the data we need to use (add pound signs in front of lines to "comment out" script so that it doesn't jam R): plants <- fread("plants_flat_file.csv") pollinator_occ <- fread("pollinator_occurrences.csv") cover <- fread("plants_early_season_cover.csv") # prep & rearrange data --------------------------------------------------- #prep some pollinator dataset metrics: #drop project 5 and 15 plants <- plants[!project %in% c(5,15)] pollinator_occ <- pollinator_occ[!Project %in% c(5,15)] # number of pollinators observed (this line counts how many rows there are that were for each pollinator species at each site-date survey) pollinator_taxa_counts <- pollinator_occ[ , .(.N) , .(Site,Date,Pollinator)] pollinator_taxa_counts_table <- pollinator_occ[ , .(.N) , .(Treatment,Pollinator)] #this line takes those "long" data and makes them "wide" pollinators_site_date <- dcast(pollinator_taxa_counts, Site + Date ~ Pollinator ) # where pollinator was blank, no bugs were seen, that got automatically named V1...change that names names(pollinators_site_date)[3] <- "no_pollinator_found" #keep goign with the comments to describe each line of code: plants$date <- anydate(plants$date) plants$date <- as.IDate(plants$date) pollinators_site_date$Date <- as.IDate(pollinators_site_date$Date) site_date_table <- plants[pollinators_site_date, on = .(site = Site, date = Date)] # survey time ------------------------------------------------------------- # survey time was same for both pollinator types on each site-date? Nope, each person stopped their #clock, meaing that we have different survey times for bees & butterflies summary(pollinator_occ[ , length(unique(`Total Survey Time (min)`)) , .(Site,Date)]) pollinator_occ[ , length(unique(`Total Survey Time (min)`)) , .(Site,Date, `Pollinator Type`)][V1>1] #luke has 2 survey times for this survey. We'll have to use the first() function to grab just the first from this time data surveytimes <- pollinator_occ[ , .(time_mins = first(`Total Survey Time (min)`)) , .(Site,Date, `Pollinator Type`)] surveytimes <- dcast(surveytimes, Site + Date ~ `Pollinator Type` ) names(surveytimes)[3:4] <- c("bee_surv_time_mins", "butterfly_survey_time_mins") surveytimes site_date_table <- site_date_table[surveytimes, on = .(site = Site, date = Date)] # Replace space with _ in columns colnames(site_date_table) <- gsub(" ", "_", colnames(site_date_table)) #Set NAs to zeros #f_dowle3natozeros(site_date_table, c(259:311)) f_dowle3natozeros(site_date_table, c("Anatrytone_logan", "Ancyloxypha_numitor", "Asterocampa_celtis", "Atalopedes_campestris", "Boloria_bellona", "Bombus_affinis", "Bombus_auricomus", "Bombus_auricomus/pensylvanicus", "Bombus_bimaculatus", "Bombus_borealis", "Bombus_citrinus", "Bombus_fervidus", "Bombus_griseocollis", "Bombus_impatiens", "Bombus_pensylvanicus", "Bombus_rufocinctus", "Bombus_sp.", "Bombus_vagans", "Celastrina_ladon", "Cercyonis_pegala", "Coenonympha_tullia", "Colias_eurytheme", "Colias_philodice", "Danaus_plexippus", "Enodia_anthedon", "Epargyreus_clarus", "Everes_comyntas", "Limenitis_archippus", "Lycaena_dione", "Lycaena_hyllus", "Lycaena_phlaeas", "Nymphalis_antiopa", "Papilio_glaucus", "Papilio_glaucus/canadensis", "Papilio_glaucus/polyxenes", "Phyciodes_selenis/tharos", "Pieris_rapae", "Polites_mystic", "Polites_peckius", "Polites_themistocles", "Pompeius_verna", "Satyrium_acadica", "Satyrium_caryaevorum", "Satyrodes_eurydice","Speyeria_cybele", "Thymelicus_lineola", "Vanessa_atalanta", "blue_sp.", "fritillary_sp.", "skipper_sp.", "sulphur_sp.")) # Add unique survey column site_date_table$survey <- paste(site_date_table$date, site_date_table$site, sep = "_") # Add round site_date_table[ , round := order(date) , by = .(site)] checks <- site_date_table[ , .(round, date)] #Pollinator Shannon Diversity #site_date_table[ , pol_shannon := diversity(site_date_table[,260:311],index = "shannon"), ] #site_date_table[ , pol_bee_shannon := diversity(site_date_table[,265:277],index = "shannon"), ] #site_date_table[ , pol_butt_shannon := diversity(site_date_table[,c(260:264, 278:311)],index = "shannon"), ] #site_date_table[ , pol_butt_shannon_nat := diversity(site_date_table[,c(260:264, 278:296, 298:305, 307:311)],index = "shannon"), ] site_date_table[ , pol_shannon := diversity(site_date_table[,c("Anatrytone_logan", "Ancyloxypha_numitor", "Asterocampa_celtis", "Atalopedes_campestris", "Boloria_bellona", "Bombus_affinis", "Bombus_auricomus", "Bombus_auricomus/pensylvanicus", "Bombus_bimaculatus", "Bombus_borealis", "Bombus_citrinus", "Bombus_fervidus", "Bombus_griseocollis", "Bombus_impatiens", "Bombus_pensylvanicus", "Bombus_rufocinctus", "Bombus_sp.", "Bombus_vagans", "Celastrina_ladon", "Cercyonis_pegala", "Coenonympha_tullia", "Colias_eurytheme", "Colias_philodice", "Danaus_plexippus", "Enodia_anthedon", "Epargyreus_clarus", "Everes_comyntas", "Limenitis_archippus", "Lycaena_dione", "Lycaena_hyllus", "Lycaena_phlaeas", "Nymphalis_antiopa", "Papilio_glaucus", "Papilio_glaucus/canadensis", "Papilio_glaucus/polyxenes", "Phyciodes_selenis/tharos", "Pieris_rapae", "Polites_mystic", "Polites_peckius", "Polites_themistocles", "Pompeius_verna", "Satyrium_acadica", "Satyrium_caryaevorum", "Satyrodes_eurydice","Speyeria_cybele", "Thymelicus_lineola", "Vanessa_atalanta", "blue_sp.", "fritillary_sp.", "skipper_sp.", "sulphur_sp.")],index = "shannon"), ] site_date_table[ , pol_bee_shannon := diversity(site_date_table[,c("Bombus_affinis", "Bombus_auricomus", "Bombus_auricomus/pensylvanicus", "Bombus_bimaculatus", "Bombus_borealis", "Bombus_citrinus", "Bombus_fervidus", "Bombus_griseocollis", "Bombus_impatiens", "Bombus_pensylvanicus", "Bombus_rufocinctus", "Bombus_sp.", "Bombus_vagans")],index = "shannon"), ] site_date_table[ , pol_butt_shannon := diversity(site_date_table[,c("Anatrytone_logan", "Ancyloxypha_numitor", "Asterocampa_celtis", "Atalopedes_campestris", "Boloria_bellona", "Celastrina_ladon", "Cercyonis_pegala", "Coenonympha_tullia", "Colias_eurytheme", "Colias_philodice", "Danaus_plexippus", "Enodia_anthedon", "Epargyreus_clarus", "Everes_comyntas", "Limenitis_archippus", "Lycaena_dione", "Lycaena_hyllus", "Lycaena_phlaeas", "Nymphalis_antiopa", "Papilio_glaucus", "Papilio_glaucus/canadensis", "Papilio_glaucus/polyxenes", "Phyciodes_selenis/tharos", "Pieris_rapae", "Polites_mystic", "Polites_peckius", "Polites_themistocles", "Pompeius_verna", "Satyrium_acadica", "Satyrium_caryaevorum", "Satyrodes_eurydice","Speyeria_cybele", "Thymelicus_lineola", "Vanessa_atalanta", "blue_sp.", "fritillary_sp.", "skipper_sp.", "sulphur_sp.")],index = "shannon"), ] site_date_table[ , pol_butt_shannon_nat := diversity(site_date_table[,c("Anatrytone_logan", "Ancyloxypha_numitor", "Asterocampa_celtis", "Atalopedes_campestris", "Boloria_bellona", "Celastrina_ladon", "Cercyonis_pegala", "Coenonympha_tullia", "Colias_eurytheme", "Colias_philodice", "Danaus_plexippus", "Enodia_anthedon", "Epargyreus_clarus", "Everes_comyntas", "Limenitis_archippus", "Lycaena_dione", "Lycaena_hyllus", "Lycaena_phlaeas", "Nymphalis_antiopa", "Papilio_glaucus", "Papilio_glaucus/canadensis", "Papilio_glaucus/polyxenes", "Phyciodes_selenis/tharos", "Polites_mystic", "Polites_peckius", "Polites_themistocles", "Pompeius_verna", "Satyrium_acadica", "Satyrium_caryaevorum", "Satyrodes_eurydice","Speyeria_cybele", "Vanessa_atalanta", "blue_sp.", "fritillary_sp.", "skipper_sp.", "sulphur_sp.")],index = "shannon"), ] #Percent quadrats with flowers over time site_date_table[ , um_percent_quad_floral := (um_n_quad_floral_res/N)] site_date_table[ , um_percent_quad_nat := (um_n_quad_floral_res_nat/N)] site_date_table[ , um_percent_quad_int := (um_n_quad_floral_res_int/N)] #Pollinator abundance #site_date_table[ , pol_abundance := rowSums(site_date_table[,260:311]) ] #site_date_table[ , pol_bee_abundance := rowSums(site_date_table[,265:277]) ] #site_date_table[ , pol_butt_abundance := rowSums(site_date_table[,c(260:264, 278:311)]) ] #site_date_table[ , pol_butt_abundance_nat := rowSums(site_date_table[,c(260:264, 278:296, 298:305, 307:311)]) ] site_date_table[ , pol_abundance := rowSums(site_date_table[,c("Anatrytone_logan", "Ancyloxypha_numitor", "Asterocampa_celtis", "Atalopedes_campestris", "Boloria_bellona", "Bombus_affinis", "Bombus_auricomus", "Bombus_auricomus/pensylvanicus", "Bombus_bimaculatus", "Bombus_borealis", "Bombus_citrinus", "Bombus_fervidus", "Bombus_griseocollis", "Bombus_impatiens", "Bombus_pensylvanicus", "Bombus_rufocinctus", "Bombus_sp.", "Bombus_vagans", "Celastrina_ladon", "Cercyonis_pegala", "Coenonympha_tullia", "Colias_eurytheme", "Colias_philodice", "Danaus_plexippus", "Enodia_anthedon", "Epargyreus_clarus", "Everes_comyntas", "Limenitis_archippus", "Lycaena_dione", "Lycaena_hyllus", "Lycaena_phlaeas", "Nymphalis_antiopa", "Papilio_glaucus", "Papilio_glaucus/canadensis", "Papilio_glaucus/polyxenes", "Phyciodes_selenis/tharos", "Pieris_rapae", "Polites_mystic", "Polites_peckius", "Polites_themistocles", "Pompeius_verna", "Satyrium_acadica", "Satyrium_caryaevorum", "Satyrodes_eurydice","Speyeria_cybele", "Thymelicus_lineola", "Vanessa_atalanta", "blue_sp.", "fritillary_sp.", "skipper_sp.", "sulphur_sp.")]) ] site_date_table[ , pol_bee_abundance := rowSums(site_date_table[,c("Bombus_affinis", "Bombus_auricomus", "Bombus_auricomus/pensylvanicus", "Bombus_bimaculatus", "Bombus_borealis", "Bombus_citrinus", "Bombus_fervidus", "Bombus_griseocollis", "Bombus_impatiens", "Bombus_pensylvanicus", "Bombus_rufocinctus", "Bombus_sp.", "Bombus_vagans")]) ] site_date_table[ , pol_butt_abundance := rowSums(site_date_table[,c("Anatrytone_logan", "Ancyloxypha_numitor", "Asterocampa_celtis", "Atalopedes_campestris", "Boloria_bellona", "Celastrina_ladon", "Cercyonis_pegala", "Coenonympha_tullia", "Colias_eurytheme", "Colias_philodice", "Danaus_plexippus", "Enodia_anthedon", "Epargyreus_clarus", "Everes_comyntas", "Limenitis_archippus", "Lycaena_dione", "Lycaena_hyllus", "Lycaena_phlaeas", "Nymphalis_antiopa", "Papilio_glaucus", "Papilio_glaucus/canadensis", "Papilio_glaucus/polyxenes", "Phyciodes_selenis/tharos", "Pieris_rapae", "Polites_mystic", "Polites_peckius", "Polites_themistocles", "Pompeius_verna", "Satyrium_acadica", "Satyrium_caryaevorum", "Satyrodes_eurydice","Speyeria_cybele", "Thymelicus_lineola", "Vanessa_atalanta", "blue_sp.", "fritillary_sp.", "skipper_sp.", "sulphur_sp.")]) ] site_date_table[ , pol_butt_abundance_nat := rowSums(site_date_table[,c("Anatrytone_logan", "Ancyloxypha_numitor", "Asterocampa_celtis", "Atalopedes_campestris", "Boloria_bellona", "Celastrina_ladon", "Cercyonis_pegala", "Coenonympha_tullia", "Colias_eurytheme", "Colias_philodice", "Danaus_plexippus", "Enodia_anthedon", "Epargyreus_clarus", "Everes_comyntas", "Limenitis_archippus", "Lycaena_dione", "Lycaena_hyllus", "Lycaena_phlaeas", "Nymphalis_antiopa", "Papilio_glaucus", "Papilio_glaucus/canadensis", "Papilio_glaucus/polyxenes", "Phyciodes_selenis/tharos", "Polites_mystic", "Polites_peckius", "Polites_themistocles", "Pompeius_verna", "Satyrium_acadica", "Satyrium_caryaevorum", "Satyrodes_eurydice","Speyeria_cybele", "Vanessa_atalanta", "blue_sp.", "fritillary_sp.", "skipper_sp.", "sulphur_sp.")]) ] # Richness site_date_table[ , pol_richness := count_row_if(gt(0),site_date_table[,c("Anatrytone_logan", "Ancyloxypha_numitor", "Asterocampa_celtis", "Atalopedes_campestris", "Boloria_bellona", "Bombus_affinis", "Bombus_auricomus", "Bombus_auricomus/pensylvanicus", "Bombus_bimaculatus", "Bombus_borealis", "Bombus_citrinus", "Bombus_fervidus", "Bombus_griseocollis", "Bombus_impatiens", "Bombus_pensylvanicus", "Bombus_rufocinctus", "Bombus_sp.", "Bombus_vagans", "Celastrina_ladon", "Cercyonis_pegala", "Coenonympha_tullia", "Colias_eurytheme", "Colias_philodice", "Danaus_plexippus", "Enodia_anthedon", "Epargyreus_clarus", "Everes_comyntas", "Limenitis_archippus", "Lycaena_dione", "Lycaena_hyllus", "Lycaena_phlaeas", "Nymphalis_antiopa", "Papilio_glaucus", "Papilio_glaucus/canadensis", "Papilio_glaucus/polyxenes", "Phyciodes_selenis/tharos", "Pieris_rapae", "Polites_mystic", "Polites_peckius", "Polites_themistocles", "Pompeius_verna", "Satyrium_acadica", "Satyrium_caryaevorum", "Satyrodes_eurydice","Speyeria_cybele", "Thymelicus_lineola", "Vanessa_atalanta", "blue_sp.", "fritillary_sp.", "skipper_sp.", "sulphur_sp.")]) ] site_date_table[ , pol_bee_richness := count_row_if(gt(0),site_date_table[,c("Bombus_affinis", "Bombus_auricomus", "Bombus_auricomus/pensylvanicus", "Bombus_bimaculatus", "Bombus_borealis", "Bombus_citrinus", "Bombus_fervidus", "Bombus_griseocollis", "Bombus_impatiens", "Bombus_pensylvanicus", "Bombus_rufocinctus", "Bombus_sp.", "Bombus_vagans")]) ] site_date_table[ , pol_butt_richness := count_row_if(gt(0),site_date_table[,c("Anatrytone_logan", "Ancyloxypha_numitor", "Asterocampa_celtis", "Atalopedes_campestris", "Boloria_bellona", "Celastrina_ladon", "Cercyonis_pegala", "Coenonympha_tullia", "Colias_eurytheme", "Colias_philodice", "Danaus_plexippus", "Enodia_anthedon", "Epargyreus_clarus", "Everes_comyntas", "Limenitis_archippus", "Lycaena_dione", "Lycaena_hyllus", "Lycaena_phlaeas", "Nymphalis_antiopa", "Papilio_glaucus", "Papilio_glaucus/canadensis", "Papilio_glaucus/polyxenes", "Phyciodes_selenis/tharos", "Pieris_rapae", "Polites_mystic", "Polites_peckius", "Polites_themistocles", "Pompeius_verna", "Satyrium_acadica", "Satyrium_caryaevorum", "Satyrodes_eurydice","Speyeria_cybele", "Thymelicus_lineola", "Vanessa_atalanta", "blue_sp.", "fritillary_sp.", "skipper_sp.", "sulphur_sp.")]) ] site_date_table[ , pol_butt_richness_nat := count_row_if(gt(0),site_date_table[,c("Anatrytone_logan", "Ancyloxypha_numitor", "Asterocampa_celtis", "Atalopedes_campestris", "Boloria_bellona", "Celastrina_ladon", "Cercyonis_pegala", "Coenonympha_tullia", "Colias_eurytheme", "Colias_philodice", "Danaus_plexippus", "Enodia_anthedon", "Epargyreus_clarus", "Everes_comyntas", "Limenitis_archippus", "Lycaena_dione", "Lycaena_hyllus", "Lycaena_phlaeas", "Nymphalis_antiopa", "Papilio_glaucus", "Papilio_glaucus/canadensis", "Papilio_glaucus/polyxenes", "Phyciodes_selenis/tharos", "Polites_mystic", "Polites_peckius", "Polites_themistocles", "Pompeius_verna", "Satyrium_acadica", "Satyrium_caryaevorum", "Satyrodes_eurydice","Speyeria_cybele", "Vanessa_atalanta", "blue_sp.", "fritillary_sp.", "skipper_sp.", "sulphur_sp.")])] #Percent of quadrats containing flowers that also have natives site_date_table[ , um_percent_quad_floral_nat := (um_n_quad_floral_res_nat/um_n_quad_floral_res)] #Mean Shannon diversity mean_um_shannon <- setDT(site_date_table)[, .(Mean = mean(um_shannons_floral), SD = sd(um_shannons_floral)), .(site,trt,age)] mean_pol_shannon <- setDT(site_date_table)[, .(Mean = mean(pol_shannon), SD = sd(pol_shannon)), .(site,trt,age)] mean_pol_butt_shannon <- setDT(site_date_table)[, .(Mean = mean(pol_butt_shannon), SD = sd(pol_butt_shannon)), .(site,trt,age)] mean_pol_bee_shannon <- setDT(site_date_table)[, .(Mean = mean(pol_bee_shannon), SD = sd(pol_bee_shannon)), .(site,trt,age)] mean_percent_floral <- setDT(site_date_table)[, .(Mean = mean(um_percent_quad_floral), SD = sd(um_percent_quad_floral)), .(site,trt,age)] mean_percent_floral_nat <- setDT(site_date_table)[, .(Mean = mean(um_percent_quad_floral_nat), SD = sd(um_percent_quad_floral_nat)), .(site,trt,age)] mean_pol_bee_abun <- setDT(site_date_table)[, .(Mean = mean(pol_bee_abundance), SD = sd(pol_bee_abundance)), .(site,trt,age)] mean_pol_butt_abun <- setDT(site_date_table)[, .(Mean = mean(pol_butt_abundance), SD = sd(pol_butt_abundance)), .(site,trt,age)] mean_pol_butt_abun_nat <- setDT(site_date_table)[, .(Mean = mean(pol_butt_abundance_nat), SD = sd(pol_butt_abundance_nat)), .(site,trt,age)] mean_pol_butt_shannon_nat <- setDT(site_date_table)[, .(Mean = mean(pol_butt_shannon_nat), SD = sd(pol_butt_shannon_nat)), .(site,trt,age)] mean_pol_butt_rich <- setDT(site_date_table)[, .(Mean = mean(pol_butt_richness)), .(site,trt,age)] mean_pol_bee_rich <- setDT(site_date_table)[, .(Mean = mean(pol_bee_richness)), .(site,trt,age)] mean_pol_butt_rich_nat <- setDT(site_date_table)[, .(Mean = mean(pol_butt_richness_nat)), .(site,trt,age)] mean_shannon <- merge(mean_um_shannon, mean_pol_shannon, by=c("site", "trt","age")) mean_shannon <- merge(mean_shannon, mean_pol_butt_shannon, by=c("site","trt","age")) mean_shannon <- merge(mean_shannon, mean_pol_bee_shannon, by=c("site","trt","age")) names(mean_shannon) <- c('site', 'trt', 'age', 'mean_shannons_floral', 'SD_shannons_floral', 'mean_pol_shannon', 'SD_pol_shannon', 'mean_pol_butt_shannon', 'SD_pol_butt_shannon', 'mean_pol_bee_shannon', 'SD_pol_bee_shannon') mean_shannon <- merge(mean_shannon, mean_percent_floral, by=c("site","trt","age")) mean_shannon <- merge(mean_shannon, mean_percent_floral_nat, by=c("site","trt","age")) names(mean_shannon) <- c('site', 'trt', 'age', 'mean_shannons_floral', 'SD_shannons_floral', 'mean_pol_shannon', 'SD_pol_shannon', 'mean_pol_butt_shannon', 'SD_pol_butt_shannon', 'mean_pol_bee_shannon', 'SD_pol_bee_shannon', 'mean_percent_floral', 'SD_percent_floral', 'mean_percent_floral_nat', 'SD_percent_floral_nat') mean_shannon <- merge(mean_shannon, mean_pol_bee_abun, by=c("site","trt","age")) mean_shannon <- merge(mean_shannon, mean_pol_butt_abun, by=c("site","trt","age")) names(mean_shannon) <- c('site', 'trt', 'age', 'mean_shannons_floral', 'SD_shannons_floral', 'mean_pol_shannon', 'SD_pol_shannon', 'mean_pol_butt_shannon', 'SD_pol_butt_shannon', 'mean_pol_bee_shannon', 'SD_pol_bee_shannon', 'mean_percent_floral', 'SD_percent_floral', 'mean_percent_floral_nat', 'SD_percent_floral_nat', 'mean_pol_bee_abun', 'SD_pol_bee_abun', 'mean_pol_butt_abun','SD_pol_butt_abun') mean_shannon <- merge(mean_shannon, mean_pol_butt_abun_nat, by=c("site","trt","age")) mean_shannon <- merge(mean_shannon, mean_pol_butt_shannon_nat, by=c("site","trt","age")) names(mean_shannon) <- c('site', 'trt', 'age', 'mean_shannons_floral', 'SD_shannons_floral', 'mean_pol_shannon', 'SD_pol_shannon', 'mean_pol_butt_shannon', 'SD_pol_butt_shannon', 'mean_pol_bee_shannon', 'SD_pol_bee_shannon', 'mean_percent_floral', 'SD_percent_floral', 'mean_percent_floral_nat', 'SD_percent_floral_nat', 'mean_pol_bee_abun', 'SD_pol_bee_abun', 'mean_pol_butt_abun','SD_pol_butt_abun', 'mean_pol_butt_abun_nat','SD_pol_butt_abun_nat', 'mean_pol_butt_shannon_nat','SD_pol_butt_shannon_nat') mean_shannon <- merge(mean_shannon, mean_pol_butt_rich, by=c("site","trt","age")) mean_shannon <- merge(mean_shannon, mean_pol_bee_rich, by=c("site","trt","age")) mean_shannon <- merge(mean_shannon, mean_pol_butt_rich_nat, by=c("site","trt","age")) names(mean_shannon) <- c('site', 'trt', 'age', 'mean_shannons_floral', 'SD_shannons_floral', 'mean_pol_shannon', 'SD_pol_shannon', 'mean_pol_butt_shannon', 'SD_pol_butt_shannon', 'mean_pol_bee_shannon', 'SD_pol_bee_shannon', 'mean_percent_floral', 'SD_percent_floral', 'mean_percent_floral_nat', 'SD_percent_floral_nat', 'mean_pol_bee_abun', 'SD_pol_bee_abun', 'mean_pol_butt_abun','SD_pol_butt_abun', 'mean_pol_butt_abun_nat','SD_pol_butt_abun_nat', 'mean_pol_butt_shannon_nat','SD_pol_butt_shannon_nat', 'mean_pol_butt_rich', 'mean_pol_bee_rich', 'mean_pol_butt_rich_nat') mean_shannons_floral_nat <- setDT(site_date_table)[, .(Mean = mean(um_shannons_floral_nat)), .(site,trt,age)] mean_shannons_floral_int <- setDT(site_date_table)[, .(Mean = mean(um_shannons_floral_int)), .(site,trt,age)] mean_shannon <- merge(mean_shannon, mean_shannons_floral_nat, by=c("site","trt","age")) mean_shannon <- merge(mean_shannon, mean_shannons_floral_int, by=c("site","trt","age")) names(mean_shannon) <- c('site', 'trt', 'age', 'mean_shannons_floral', 'SD_shannons_floral', 'mean_pol_shannon', 'SD_pol_shannon', 'mean_pol_butt_shannon', 'SD_pol_butt_shannon', 'mean_pol_bee_shannon', 'SD_pol_bee_shannon', 'mean_percent_floral', 'SD_percent_floral', 'mean_percent_floral_nat', 'SD_percent_floral_nat', 'mean_pol_bee_abun', 'SD_pol_bee_abun', 'mean_pol_butt_abun','SD_pol_butt_abun', 'mean_pol_butt_abun_nat','SD_pol_butt_abun_nat', 'mean_pol_butt_shannon_nat','SD_pol_butt_shannon_nat', 'mean_pol_butt_rich', 'mean_pol_bee_rich', 'mean_pol_butt_rich_nat', 'mean_shannons_floral_nat', 'mean_shannons_floral_int') mean_percent_int <- setDT(site_date_table)[, .(Mean = mean(um_percent_quad_int)), .(site,trt,age)] mean_shannon <- merge(mean_shannon, mean_percent_int, by=c("site","trt","age")) mean_percent_nat <- setDT(site_date_table)[, .(Mean = mean(um_percent_quad_nat)), .(site,trt,age)] mean_shannon <- merge(mean_shannon, mean_percent_nat, by=c("site","trt","age")) names(mean_shannon) <- c('site', 'trt', 'age', 'mean_shannons_floral', 'SD_shannons_floral', 'mean_pol_shannon', 'SD_pol_shannon', 'mean_pol_butt_shannon', 'SD_pol_butt_shannon', 'mean_pol_bee_shannon', 'SD_pol_bee_shannon', 'mean_percent_floral', 'SD_percent_floral', 'mean_percent_floral_nat', 'SD_percent_floral_nat', 'mean_pol_bee_abun', 'SD_pol_bee_abun', 'mean_pol_butt_abun','SD_pol_butt_abun', 'mean_pol_butt_abun_nat','SD_pol_butt_abun_nat', 'mean_pol_butt_shannon_nat','SD_pol_butt_shannon_nat', 'mean_pol_butt_rich', 'mean_pol_bee_rich', 'mean_pol_butt_rich_nat', 'mean_shannons_floral_nat', 'mean_shannons_floral_int', 'mean_percent_int', 'mean_percent_nat') # Presence/absence site_date_table$presence_absence_bee <- ifelse(site_date_table$pol_bee_abundance > 0, 1, 0) site_date_table$presence_absence_butt <- ifelse(site_date_table$pol_butt_abundance > 0, 1, 0) site_date_table$presence_absence_butt_nat <- ifelse(site_date_table$pol_butt_abundance_nat > 0, 1, 0) # Subset data for only presence site_date_table_bee <- subset(site_date_table, presence_absence_bee == 1) site_date_table_butt <- subset(site_date_table, presence_absence_butt == 1) site_date_table_butt_nat <- subset(site_date_table, presence_absence_butt_nat == 1) #Questions--------------------- # All questions in one model and check collineaerity ---- butt.m.all <- glmer(presence_absence_butt ~ trt + um_shannons_floral + um_percent_quad_floral + butterfly_survey_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) butt_nat.m.all <- glmer(presence_absence_butt_nat ~ trt + um_shannons_floral + um_percent_quad_floral + um_percent_quad_floral_nat + butterfly_survey_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) bee.m.all <- glmer(presence_absence_bee ~ trt + um_shannons_floral + um_percent_quad_floral + um_percent_quad_floral_nat + bee_surv_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) vif(butt.m.all) vif(butt_nat.m.all) vif(bee.m.all) butt.m.all2 <- lme(pol_butt_shannon ~ trt + um_shannons_floral + um_percent_quad_floral + butterfly_survey_time_mins, random = ~1|round/project, data = site_date_table_butt, na.action = na.omit) butt_nat.m.all2 <- lme(pol_butt_shannon_nat ~ trt + um_shannons_floral + um_percent_quad_floral + um_percent_quad_floral_nat + butterfly_survey_time_mins, random = ~1|round/project, data = site_date_table_butt_nat, na.action = na.omit) bee.m.all2 <- lme(pol_bee_shannon ~ trt + um_shannons_floral + um_percent_quad_floral + um_percent_quad_floral_nat + bee_surv_time_mins, random = ~1|round/project, data = site_date_table_bee, na.action = na.omit) vif(butt.m.all2) vif(butt_nat.m.all2) vif(bee.m.all2) summary(butt.m.all2) summary(butt_nat.m.all2) summary(bee.m.all2) butt_abun.m.all <- glmmTMB(pol_butt_abundance ~ trt + um_shannons_floral + um_percent_quad_floral + butterfly_survey_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) butt_abun_nat.m.all <- glmmTMB(pol_butt_abundance_nat ~ trt + um_shannons_floral + um_percent_quad_floral + um_percent_quad_floral_nat + butterfly_survey_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) bee_abun.m.all <- glmmTMB(pol_bee_abundance ~ trt + um_shannons_floral + um_percent_quad_floral + um_percent_quad_floral_nat + bee_surv_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) check_collinearity(butt_abun.m.all) check_collinearity(butt_abun_nat.m.all) check_collinearity(bee_abun.m.all) # Plant summary ---- # Average plot size area_summary <- site_date_table %>% group_by(site, trt) %>% filter(!is.na(area_m2)) %>% summarize(area = mean(area_m2)) area_summary %>% ungroup() %>% summarize(mean_area = mean(area), median_area = median(area), sd_area = sd(area)) %>% as.data.frame(.) %>% mutate_if(is.numeric, round, 1) area_summary %>% group_by(trt) %>% summarize(mean_area = mean(area), median_area = median(area), sd_area = sd(area)) %>% as.data.frame(.) %>% mutate_if(is.numeric, round, 1) # Number of quadrats site_date_table %>% filter(!is.na(N)) %>% summarize(mean_quad = mean(N), median_quad = median(N), max_quad = max(N), min_quad = min(N)) # Early cover cover <- cover %>% mutate(year_category = case_when(year >= 2013 ~ "0-8", year <= 2012 ~ "9-20")) cover$year_category[cover$trt == "typ"] <- NA cover_long <- cover %>% dplyr::select(!(notes:totcover)) %>% rename(grass = grass_p, ground = ground_p, forb = forb_p) %>% pivot_longer(cols = grass:forb, names_to = "cover_class", values_to = "cover_p") # pdf("cover.pdf", width=6, height=5) cover_long %>% ggplot(aes(x = 2021 - year, y = cover_p, col = cover_class)) + geom_point(alpha = 0.1) + geom_smooth(method = "lm") + labs(y = "Proportion cover per quadrat", x = "Plot age", color = "Cover type") + theme_bw() + scale_color_manual(values = c("lightcoral", "darkgreen", "burlywood4")) # dev.off() # pdf("cover_boxplot.pdf", width=6, height=5) cover_long %>% ggplot(aes(x = trt, y = cover_p, fill = cover_class)) + geom_boxplot() + labs(x = "Seeding treatment", y = "Proportion cover", fill = "Cover type") + scale_x_discrete(labels = c("Native\nseed mix", "Non-native\nseed mix", "Typical\nroadside")) + stat_summary(fun="mean", shape = 23, position = position_dodge(width = .75)) + theme_bw() + scale_fill_manual(values = c("lightcoral", "darkgreen", "burlywood4")) # dev.off() cover_summary <- cover %>% group_by(trt, year_category) %>% summarize(n = n(), mean_grass = mean(grass_p), sd_grass = sd(grass_p), mean_forb = mean(forb_p), sd_forb = sd(forb_p), mean_ground = mean(ground_p), sd_ground = sd(ground_p)) # Plant community summary table # forb richness for different native and non-native seed mixes seed_mix_rich_young_nat <- c(20, 18, 18, 18, 20, 11, 18, 11, 11, 18) seed_mix_rich_young_non <- c(0, 2, 2, 2, 2, 2, 0, 2, 2, 2) seed_mix_rich_old_nat <- c(20, 6, 20, 20, 18, 20, 18, 9, 11) seed_mix_rich_old_non <- c(2, 2, 2, 2, 2, 2, 4, 2, 2) mean(seed_mix_rich_young_nat) sd(seed_mix_rich_young_nat) mean(seed_mix_rich_young_non) sd(seed_mix_rich_young_non) mean(seed_mix_rich_old_nat) sd(seed_mix_rich_old_nat) mean(seed_mix_rich_old_non) sd(seed_mix_rich_old_non) site_date_table <- site_date_table %>% mutate(year_category = case_when(Year >= 2013 ~ "0-8", Year <= 2012 ~ "9-20")) site_date_table$year_category[site_date_table$trt == "typ"] <- NA # Summarize floral counts by site (cumulative richness at each site) cum_richness <- site_date_table %>% select(site, trt, year_category, c(Achillea_millefolium:Alliaria_petiolata), c(Anemone_canadensis:Berteroa_incana), c(Carduus_acanthoides:Echinocystis_lobata), c(Epilobium_ciliatum:Oxalis_dillenii), c(Pastinaca_sativa:Persicaria_amphibia), Plantago_lanceolata, c(Polygonum_spp.:Saponaria_officinalis), Securigera_varia, c(Silene_csereii:Sonchus_arvensis), c(Spergularia_rubra:Zizia_aurea)) %>% group_by(site, trt, year_category) %>% summarize_if(is.numeric, sum, na.rm = TRUE) cum_richness <- cum_richness %>% ungroup() %>% mutate(richness = rowSums(.[4:123]!=0)) cum_richness %>% group_by(trt, year_category) %>% summarize(mean_richness = mean(richness), sd_richness = sd(richness)) cum_richness_nat <- cum_richness %>% select(site, trt, year_category, Achillea_millefolium, Agastache_foeniculum, Anemone_canadensis, Apocynum_cannabinum, c(Asclepias_syriaca:Astragalus_canadensis), Cirsium_discolor, Comandra_umbellata, Conyza_canadensis, Coreopsis_leavenworthii, Dalea_candida, Dalea_purpurea, c(Desmodium_canadense:Euphorbia_corollata), Geranium_maculatum, Geum_aleppicum, c(Grindelia_hirsutula:Heliopsis_helianthoides), Impatiens_capensis, Lobelia_siphilitica, Lycopus_americanus, Mirabilis_nyctaginea, Monarda_fistulosa, c(Oenothera_biennis:Oxalis_dillenii), Persicaria_amphibia, Polygonum_spp., Potentilla_norvegica, Pseudognaphalium_obtusifolium, c(Ratibida_columnifera:Rudbeckia_hirta), Sambucus_canadensis, Silphium_perfoliatum, Solidago_spp., Stachys_palustris, Teucrium_canadense, Thalictrum_dasycarpum, c(Verbena_hastata:Vernonia_fasciculata), Zizia_aurea, Symphyotrichum_novae_angliae) cum_richness_nat <- cum_richness_nat %>% ungroup() %>% mutate(richness = rowSums(.[4:62]!=0)) cum_richness_nat %>% group_by(trt, year_category) %>% summarize(mean_richness_nat = mean(richness), sd_richness_nat = sd(richness)) cum_richness_non <- cum_richness %>% select(-c(richness, Achillea_millefolium, Agastache_foeniculum, Anemone_canadensis, Apocynum_cannabinum, c(Asclepias_syriaca:Astragalus_canadensis), Cirsium_discolor, Comandra_umbellata, Conyza_canadensis, Coreopsis_leavenworthii, Dalea_candida, Dalea_purpurea, c(Desmodium_canadense:Euphorbia_corollata), Geranium_maculatum, Geum_aleppicum, c(Grindelia_hirsutula:Heliopsis_helianthoides), Impatiens_capensis, Lobelia_siphilitica, Lycopus_americanus, Mirabilis_nyctaginea, Monarda_fistulosa, c(Oenothera_biennis:Oxalis_dillenii), Persicaria_amphibia, Polygonum_spp., Potentilla_norvegica, Pseudognaphalium_obtusifolium, c(Ratibida_columnifera:Rudbeckia_hirta), Sambucus_canadensis, Silphium_perfoliatum, Solidago_spp., Stachys_palustris, Teucrium_canadense, Thalictrum_dasycarpum, c(Verbena_hastata:Vernonia_fasciculata), Zizia_aurea, Symphyotrichum_novae_angliae)) cum_richness_non <- cum_richness_non %>% ungroup() %>% mutate(richness = rowSums(.[4:64]!=0)) cum_richness_non %>% group_by(trt, year_category) %>% summarize(mean_richness_non = mean(richness), sd_richness_non = sd(richness)) plant_summary <- site_date_table %>% ungroup() %>% filter(!is.na(N)) %>% group_by(trt, year_category) %>% summarize(n = n(), mean_floral_richness = mean(um_forb_richness), sd_floral_richness = sd(um_forb_richness), mean_native_floral_richness = mean(um_richness_floral_nat), sd_native_floral_richness = sd(um_richness_floral_nat), mean_non_floral_richness = mean(um_richness_floral_int), sd_non_floral_richness = sd(um_richness_floral_int), mean_floral_shannon = mean(um_shannons_floral), sd_floral_shannon = sd(um_shannons_floral), mean_native_floral_shannon = mean(um_shannons_floral_nat), sd_native_floral_shannon = sd(um_shannons_floral_nat), mean_non_floral_shannon = mean(um_shannons_floral_int), sd_non_floral_shannon = sd(um_shannons_floral_int), mean_proportion_flowers = mean(um_percent_quad_floral), sd_proportion_flowers = sd(um_percent_quad_floral)) plant_summary2 <- site_date_table %>% filter(!is.na(N)) %>% group_by(site, trt) %>% summarize(mean_floral_richness = mean(um_forb_richness), mean_proportion_flowers = mean(um_percent_quad_floral)) # Question 4 ----- # Formerly Q1 # Explore data (both diversity and abundance zero-inflated) hist(site_date_table$pol_butt_shannon) hist(site_date_table$pol_butt_shannon_nat) hist(site_date_table$pol_bee_shannon) hist(site_date_table$pol_butt_abundance) hist(site_date_table$pol_butt_abundance_nat) hist(site_date_table$pol_bee_abundance) # Q1 diversity models # Did NOT drop zeros (these three models are not used in manuscript) Q1_butt.m <- lme(pol_butt_shannon ~ trt + butterfly_survey_time_mins, random = ~1|project/site, data = site_date_table, na.action = na.omit) Q1_butt_nat.m <- lme(pol_butt_shannon_nat ~ trt + butterfly_survey_time_mins, random = ~1|project/site, data = site_date_table, na.action = na.omit) Q1_bee.m <- lme(pol_bee_shannon ~ trt + bee_surv_time_mins, random = ~1|project/site, data = site_date_table, na.action = na.omit) # Check models (zero-inflated, bad) plot(Q1_butt.m) hist(resid(Q1_butt.m)) plot(Q1_butt_nat.m) hist(resid(Q1_butt_nat.m)) plot(Q1_bee.m) hist(resid(Q1_bee.m)) # Q1 drop zeros # Presence/absence Q1_butt.m1 <- glmer(presence_absence_butt ~ trt + butterfly_survey_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) # Check singularity (it's fine; 0.001) tt <- getME(Q1_butt.m1,"theta") ll <- getME(Q1_butt.m1,"lower") min(tt[ll==0]) # Check gradient calculations (0.00106 is larger than the tolerance of 0.001) derivs1 <- Q1_butt.m1@optinfo$derivs sc_grad1 <- with(derivs1,solve(Hessian,gradient)) max(abs(sc_grad1)) max(pmin(abs(sc_grad1),abs(derivs1$gradient))) # Redo gradient calculations to double check (not that different) dd <- update(Q1_butt.m1,devFunOnly=TRUE) pars <- unlist(getME(Q1_butt.m1,c("theta","fixef"))) grad2 <- grad(dd,pars) hess2 <- hessian(dd,pars) sc_grad2 <- solve(hess2,grad2) max(pmin(abs(sc_grad2),abs(grad2))) # Restart model with max interations ss <- getME(Q1_butt.m1,c("theta","fixef")) # Model works, but singular fit Q1_butt.m1.1 <- update(Q1_butt.m1,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4))) Q1_butt_nat.m1 <- glmer(presence_absence_butt_nat ~ trt + butterfly_survey_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) Q1_bee.m1 <- glmer(presence_absence_bee ~ trt + bee_surv_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) # Anovas presence/absence (type II for unbalanced ancovas) Anova(Q1_butt.m1.1, type = 2) Anova(Q1_butt_nat.m1, type = 2) Anova(Q1_bee.m1, type = 2) # Post-hoc for Q1 bee presence/abundance emm_options(opt.digits = FALSE) emmeans(Q1_bee.m1, pairwise ~ trt, adjust = "tukey") # Diversity with no zeros Q1_butt.m2 <- lme(pol_butt_shannon ~ trt + butterfly_survey_time_mins, random = ~1|round/project, data = site_date_table_butt, na.action = na.omit) Q1_butt_nat.m2 <- lme(pol_butt_shannon_nat ~ trt + butterfly_survey_time_mins, random = ~1|round/project, data = site_date_table_butt_nat, na.action = na.omit) Q1_bee.m2 <- lme(pol_bee_shannon ~ trt + bee_surv_time_mins, random = ~1|round/project, data = site_date_table_bee, na.action = na.omit) # Anovas diversity (type II for unbalanced ancovas) Anova(Q1_butt.m2, type = 2) Anova(Q1_butt_nat.m2, type = 2) Anova(Q1_bee.m2, type = 2) # Post-hoc for Q1 butterfly diversity emmeans(Q1_butt.m2, pairwise ~ trt, adjust = "tukey") # Q1 abundance models (zero-inflated) Q1_butt_abun.m <- glmmTMB(pol_butt_abundance ~ trt + butterfly_survey_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) Q1_butt_abun_nat.m <- glmmTMB(pol_butt_abundance_nat ~ trt + butterfly_survey_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) Q1_bee_abun.m <- glmmTMB(pol_bee_abundance ~ trt + bee_surv_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) # Anovas abundance (type II for unbalanced ancovas) Anova(Q1_butt_abun.m, type = 2) Anova(Q1_butt_abun_nat.m, type = 2) Anova(Q1_bee_abun.m, type = 2) # Post-hoc abundance emmeans(Q1_butt_abun_nat.m, pairwise ~ trt, adjust = "tukey") emmeans(Q1_bee_abun.m, pairwise ~ trt, adjust = "tukey") # Q1 boxplots boxplot_data <- mean_shannon %>% gather(pol_type, pol_diversity, mean_pol_butt_shannon_nat, mean_pol_bee_shannon) boxplot_data2 <- mean_shannon %>% gather(pol_type, pol_abundance, mean_pol_butt_abun_nat, mean_pol_bee_abun) # pdf("boxplot_div.pdf", width=6, height=5) ggplot(data = subset(boxplot_data, !is.na(trt)), aes(x = trt, y = pol_diversity, fill = factor(pol_type, labels = c("Bumble bee", "Butterfly"))))+ geom_boxplot()+ theme_classic()+ labs(x = "Seeding treatment", y= "Native Pollinator Shannon Diversity Index", fill = "Pollinator Type")+ scale_x_discrete(labels = c("Native\nseed mix", "Non-native\nseed mix", "Typical\nroadside"))+ stat_summary(fun="mean", shape = 23, position = position_dodge(width = .75)) + scale_fill_grey(start=0.8, end = 0.4) # dev.off() # pdf("boxplot_abun.pdf", width=6, height=5) ggplot(data = subset(boxplot_data2, !is.na(trt)), aes(x = trt, y = pol_abundance, fill = factor(pol_type, labels = c("Bumble bee", "Butterfly"))))+ geom_boxplot()+ theme_classic()+ labs(x = "Seeding treatment", y= "Native Pollinator Abundance", fill = "Pollinator Type")+ scale_x_discrete(labels = c("Native\nseed mix", "Non-native\nseed mix", "Typical\nroadside"))+ stat_summary(fun="mean", shape = 23, position = position_dodge(width = .75)) + scale_fill_grey(start=0.8, end = 0.4) # dev.off() # Question 1 ----- # Formerly Q2 # Assign colors for plots colors <- c("Bumble bee" = "red", "Butterfly" = "blue") # Explore data # Native floral diversity vs pollinator diversity ggplot(data = mean_shannon, aes(x = mean_shannons_floral_nat))+ geom_point(aes(y = mean_pol_bee_shannon, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_shannon, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_shannon, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_shannon, color="Butterfly"), linetype = "dashed")+ labs(x="Native Floral Shannon Diversity Index", y="Pollinator Shannon Diversity Index", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # Introduced floral diversity vs pollinator diversity ggplot(data = mean_shannon, aes(x = mean_shannons_floral_int))+ geom_point(aes(y = mean_pol_bee_shannon, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_shannon, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_shannon, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_shannon, color="Butterfly"), linetype = "dashed")+ labs(x="Introduced Floral Shannon Diversity Index", y="Pollinator Shannon Diversity Index", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # Q2 presence/absence Q2_butt.m <- glmer(presence_absence_butt ~ um_shannons_floral + butterfly_survey_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) # Check singularity (it's fine; 0.002) tt <- getME(Q2_butt.m,"theta") ll <- getME(Q2_butt.m,"lower") min(tt[ll==0]) # Check gradient calculations (0.002 is larger than the tolerance of 0.001) derivs1 <- Q2_butt.m@optinfo$derivs sc_grad1 <- with(derivs1,solve(Hessian,gradient)) max(abs(sc_grad1)) max(pmin(abs(sc_grad1),abs(derivs1$gradient))) # Redo gradient calculations to double check (not that different) dd <- update(Q2_butt.m,devFunOnly=TRUE) pars <- unlist(getME(Q2_butt.m,c("theta","fixef"))) grad2 <- grad(dd,pars) hess2 <- hessian(dd,pars) sc_grad2 <- solve(hess2,grad2) max(pmin(abs(sc_grad2),abs(grad2))) # Restart model with max interations ss <- getME(Q2_butt.m,c("theta","fixef")) # Model works, but singular fit Q2_butt.m1 <- update(Q2_butt.m,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4))) Q2_butt_nat.m <- glmer(presence_absence_butt_nat ~ um_shannons_floral + butterfly_survey_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) Q2_bee.m <- glmer(presence_absence_bee ~ um_shannons_floral + bee_surv_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) # Q2 diversity models Q2_butt.m2 <- lme(pol_butt_shannon ~ um_shannons_floral + butterfly_survey_time_mins, random = ~1|round/project, data = site_date_table_butt, na.action = na.omit) Q2_butt_nat.m2 <- lme(pol_butt_shannon_nat ~ um_shannons_floral + butterfly_survey_time_mins, random = ~1|round/project, data = site_date_table_butt_nat, na.action = na.omit) Q2_bee.m2 <- lme(pol_bee_shannon ~ um_shannons_floral + bee_surv_time_mins, random = ~1|round/project, data = site_date_table_bee, na.action = na.omit) # Q2 abundance models Q2_butt_abun.m <- glmmTMB(pol_butt_abundance ~ um_shannons_floral + butterfly_survey_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) Q2_butt_abun_nat.m <- glmmTMB(pol_butt_abundance_nat ~ um_shannons_floral + butterfly_survey_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) Q2_bee_abun.m <- glmmTMB(pol_bee_abundance ~ um_shannons_floral + bee_surv_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) # Q2 plots # jpeg("diversity.jpeg", units="in", width=6, height=5, res=300) ggplot(data = mean_shannon, aes(x = mean_shannons_floral))+ geom_point(aes(y = mean_pol_bee_shannon, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_shannon, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_shannon, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_shannon, color="Butterfly"), linetype = "dashed")+ labs(x="Floral Shannon Diversity Index", y="Pollinator Shannon Diversity Index", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # dev.off() # pdf("diversity_nat.pdf", width=6, height=5) ggplot(data = mean_shannon, aes(x = mean_shannons_floral))+ geom_point(aes(y = mean_pol_bee_shannon, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_shannon_nat, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_shannon, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_shannon_nat, color="Butterfly"), linetype = "dashed")+ labs(x="Floral Shannon Diversity Index", y="Native Pollinator Shannon Diversity Index", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # dev.off() # pdf("diversity_abun.pdf", width=6, height=5) ggplot(data = mean_shannon, aes(x = mean_shannons_floral))+ geom_point(aes(y = mean_pol_bee_abun, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_abun, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_abun, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_abun, color="Butterfly"), linetype = "dashed")+ labs(x="Floral Shannon Diversity Index", y="Pollinator Abundance", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # dev.off() # pdf("diversity_abun_nat.pdf", width=6, height=5) ggplot(data = mean_shannon, aes(x = mean_shannons_floral))+ geom_point(aes(y = mean_pol_bee_abun, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_abun_nat, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_abun, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_abun_nat, color="Butterfly"), linetype = "dashed")+ labs(x="Floral Shannon Diversity Index", y="Native Pollinator Abundance", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # dev.off() # Question 2 ----- # Formerly Q3 # Explore data # Native flowers ggplot(data = mean_shannon, aes(x = mean_percent_nat))+ geom_point(aes(y = mean_pol_bee_shannon, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_shannon, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_shannon, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_shannon, color="Butterfly"), linetype = "dashed")+ labs(x="Proportion of quadrats containing native flowers", y="Pollinator Shannon Diversity Index", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # Introduced flowers ggplot(data = mean_shannon, aes(x = mean_percent_int))+ geom_point(aes(y = mean_pol_bee_shannon, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_shannon, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_shannon, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_shannon, color="Butterfly"), linetype = "dashed")+ labs(x="Proportion of quadrats containing introduced flowers", y="Pollinator Shannon Diversity Index", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # Q3 Presence/absence models Q3_butt.m <- glmer(presence_absence_butt ~ um_percent_quad_floral + butterfly_survey_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) Q3_butt_nat.m <- glmer(presence_absence_butt_nat ~ um_percent_quad_floral + butterfly_survey_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) Q3_bee.m <- glmer(presence_absence_bee ~ um_percent_quad_floral + bee_surv_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) # Q3 diversity models (no zeros) Q3_butt.m2 <- lme(pol_butt_shannon ~ um_percent_quad_floral + butterfly_survey_time_mins, random = ~1|round/project, data = site_date_table_butt, na.action = na.omit) Q3_butt_nat.m2 <- lme(pol_butt_shannon_nat ~ um_percent_quad_floral + butterfly_survey_time_mins, random = ~1|round/project, data = site_date_table_butt_nat, na.action = na.omit) Q3_bee.m2 <- lme(pol_bee_shannon ~ um_percent_quad_floral + bee_surv_time_mins, random = ~1|round/project, data = site_date_table_bee, na.action = na.omit) # Q3 abundance models Q3_butt_abun.m <- glmmTMB(pol_butt_abundance ~ um_percent_quad_floral + butterfly_survey_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) Q3_butt_abun_nat.m <- glmmTMB(pol_butt_abundance_nat ~ um_percent_quad_floral + butterfly_survey_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) Q3_bee_abun.m <- glmmTMB(pol_bee_abundance ~ um_percent_quad_floral + bee_surv_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) # Q3 plots # pdf("perquad_diversity.pdf", width=6, height=5) ggplot(data = mean_shannon, aes(x = mean_percent_floral))+ geom_point(aes(y = mean_pol_bee_shannon, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_shannon, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_shannon, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_shannon, color="Butterfly"), linetype = "dashed")+ labs(x="Proportion of quadrats containing all flowers", y="Pollinator Shannon Diversity Index", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # dev.off() # pdf("perquad_diversity_nat.pdf", width=6, height=5) ggplot(data = mean_shannon, aes(x = mean_percent_floral))+ geom_point(aes(y = mean_pol_bee_shannon, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_shannon_nat, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_shannon, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_shannon_nat, color="Butterfly"), linetype = "dashed")+ labs(x="Proportion of quadrats containing flowers", y="Native Pollinator Shannon Diversity Index", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # dev.off() # jpeg("perquad_abundance.jpeg", units="in", width=6, height=5, res=300) ggplot(data = mean_shannon, aes(x = mean_percent_floral))+ geom_point(aes(y = mean_pol_bee_abun, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_abun, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_abun, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_abun, color="Butterfly"), linetype = "dashed")+ labs(x="Proportion of quadrats containing flowers", y="Pollinator Abundance", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # dev.off() # pdf("perquad_abundance_nat.pdf", width=6, height=5) ggplot(data = mean_shannon, aes(x = mean_percent_floral))+ geom_point(aes(y = mean_pol_bee_abun, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_abun_nat, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_abun, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_abun_nat, color="Butterfly"), linetype = "dashed")+ labs(x="Proportion of quadrats containing flowers", y="Native Pollinator Abundance", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # dev.off() #Question 3 ----- # Formerly Q4 # Q4 Presence/absence models Q4_butt_nat.m <- glmer(presence_absence_butt_nat ~ um_percent_quad_floral_nat + butterfly_survey_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) Q4_bee.m <- glmer(presence_absence_bee ~ um_percent_quad_floral_nat + bee_surv_time_mins + (1|round/project), family = "binomial", data = site_date_table, na.action = na.omit) # Q4 diversity models Q4_butt_nat.m2 <- lme(pol_butt_shannon_nat ~ um_percent_quad_floral_nat + butterfly_survey_time_mins, random = ~1|round/project, data = site_date_table_butt_nat, na.action = na.omit) Q4_bee.m2 <- lme(pol_bee_shannon ~ um_percent_quad_floral_nat + bee_surv_time_mins, random = ~1|round/project, data = site_date_table_bee, na.action = na.omit) # Q4 abundance models Q4_butt_abun_nat.m <- glmmTMB(pol_butt_abundance_nat ~ um_percent_quad_floral_nat + butterfly_survey_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) Q4_bee_abun.m <- glmmTMB(pol_bee_abundance ~ um_percent_quad_floral_nat + bee_surv_time_mins + (1|round/project), family = "nbinom2", data = site_date_table, na.action = na.omit) # Q4 plots # pdf("nat_diversity.pdf", width=6, height=5) ggplot(data = mean_shannon, aes(x = mean_percent_floral_nat))+ geom_point(aes(y = mean_pol_bee_shannon, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_shannon_nat, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_shannon, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_shannon_nat, color="Butterfly"), linetype = "dashed")+ labs(x="Proportion of flowers that are native", y="Native Pollinator Shannon Diversity Index", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # dev.off() # pdf("nat_abundance.pdf", width=6, height=5) ggplot(data = mean_shannon, aes(x = mean_percent_floral_nat))+ geom_point(aes(y = mean_pol_bee_abun, color = "Bumble bee")) + geom_point(aes(y = mean_pol_butt_abun_nat, color="Butterfly"))+ geom_smooth(method = "lm", aes(y = mean_pol_bee_abun, color = "Bumble bee"))+ geom_smooth(method = "lm", aes(y = mean_pol_butt_abun_nat, color="Butterfly"), linetype = "dashed")+ labs(x="Proportion of flowers that are native", y="Native Pollinator Abundance", color = "Pollinator Type")+ theme_classic()+ scale_color_manual(values = colors) # dev.off() # Site age ----- # Exploratory plots # jpeg("bee_age.jpeg", units="in", width=6, height=5, res=300) ggplot(data = subset(mean_shannon, !is.na(age)), aes(x = age, y = mean_pol_bee_shannon, group = trt, color = trt))+ geom_point()+ geom_smooth(method = "lm")+ labs(x = "Age (years)", y = "Bumble bee Shannon Diversity Index", color = "Seeding Treatment")+ theme_classic()+ scale_color_discrete(labels = c("Native seed mix", "Non-native seed mix")) # dev.off() # jpeg("butt_age.jpeg", units="in", width=6, height=5, res=300) ggplot(data = subset(mean_shannon, !is.na(age)), aes(x = age, y = mean_pol_butt_shannon, group = trt, color = trt))+ geom_point()+ geom_smooth(method = "lm")+ labs(x = "Age (years)", y = "Butterfly Shannon Diversity Index", color = "Seeding Treatment")+ theme_classic()+ scale_color_discrete(labels = c("Native seed mix", "Non-native seed mix")) # dev.off() # jpeg("butt_age_nat.jpeg", units="in", width=6, height=5, res=300) ggplot(data = subset(mean_shannon, !is.na(age)), aes(x = age, y = mean_pol_butt_shannon_nat, group = trt, color = trt))+ geom_point()+ geom_smooth(method = "lm")+ labs(x = "Age (years)", y = "Native Butterfly Shannon Diversity Index", color = "Seeding Treatment")+ theme_classic()+ scale_color_discrete(labels = c("Native seed mix", "Non-native seed mix")) # dev.off() # jpeg("bee_age_abun.jpeg", units="in", width=6, height=5, res=300) ggplot(data = subset(mean_shannon, !is.na(age)), aes(x = age, y = mean_pol_bee_abun, group = trt, color = trt))+ geom_point()+ geom_smooth(method = "lm")+ labs(x = "Age (years)", y = "Bumble bee Abundance", color = "Seeding Treatment")+ theme_classic()+ scale_color_discrete(labels = c("Native seed mix", "Non-native seed mix")) # dev.off() # jpeg("butt_age_abun.jpeg", units="in", width=6, height=5, res=300) ggplot(data = subset(mean_shannon, !is.na(age)), aes(x = age, y = mean_pol_butt_abun, group = trt, color = trt))+ geom_point()+ geom_smooth(method = "lm")+ labs(x = "Age (years)", y = "Butterfly Abundance", color = "Seeding Treatment")+ theme_classic()+ scale_color_discrete(labels = c("Native seed mix", "Non-native seed mix")) # dev.off() # jpeg("butt_age_abun_nat.jpeg", units="in", width=6, height=5, res=300) ggplot(data = subset(mean_shannon, !is.na(age)), aes(x = age, y = mean_pol_butt_abun_nat, group = trt, color = trt))+ geom_point()+ geom_smooth(method = "lm")+ labs(x = "Age (years)", y = "Native Butterfly Abundance", color = "Seeding Treatment")+ theme_classic()+ scale_color_discrete(labels = c("Native seed mix", "Non-native seed mix")) # dev.off() # Seasonal pollinatory diversity ----- # Exploratory plots # jpeg("bee_seasonal.jpeg", units="in", width=6, height=5, res=300) ggplot(data = subset(site_date_table, !is.na(trt)), aes(x = date, y = pol_bee_shannon, group = trt, color = trt))+ geom_point()+ geom_smooth()+ labs(x = "Date", y = "Bumble bee Shannon Diversity Index", color = "Seeding Treatment")+ theme_classic() # dev.off() # jpeg("butt_seasonal.jpeg", units="in", width=6, height=5, res=300) ggplot(data = subset(site_date_table, !is.na(trt)), aes(x = date, y = pol_butt_shannon, group = trt, color = trt))+ geom_point()+ geom_smooth()+ labs(x = "Date", y = "Butterfly Shannon Diversity Index", color = "Seeding Treatment")+ theme_classic()+ scale_color_discrete(labels = c("Native seed mix", "Non-native seed mix", "Typical roadside")) # dev.off() # Occurrences tables ----- report_table <- pollinator_taxa_counts_table %>% group_by(Pollinator) %>% summarise(Frequency = sum(N))