--- title: "tiede2022code" author: "Ben Iuliano" date: "12/2/2021" output: html_document --- #Setup ##Install packages ```{r} library(tidyverse) library(ggeffects) library(fitdistrplus) library(car) library(effects) library(nlme) library(Hmisc) library(MuMIn) library(sf) library(GGally) library(wec) library(effectsize) library(glmmTMB) library(DHARMa) library(ggpubr) library(RColorBrewer) library(latticeExtra) ``` ##Color settings ```{r} myPalette <- c("#009E73", "#D55E00", "#E69F00", "#CC79A7", "#0072B2") myPalette2 <- c("#000000", "#009E73", "#D55E00", "#E69F00", "#CC79A7", "#0072B2") myPalette2.rev <- rev(myPalette2) myPalette.noC7 <- c("#009E73", "#D55E00", "#E69F00", "#0072B2") ``` ##Load data ```{r} beetle.data <- read.csv("beetle.data.csv", stringsAsFactors = T) landscape.data <- read.csv("landscape.data.csv", stringsAsFactors = T) data <- merge(beetle.data, landscape.data) ``` #Main Analysis ##Create standardized dataframes for each body condition metric ###Body size ```{r} size <- droplevels(subset(data, size>0)) #remove missing values size.z <- data.frame(crop = standardize(size$crop), sidi = standardize(size$sidi), ed = standardize(size$ed), species = size$species, species.wec.1 = size$species, #for weighted effects coding species.wec.2 = size$species, #for weighted effects coding cd = size$cd, siteID = size$siteID, sex = size$sex, size.z = standardize(size$size)) #standardize #weighted effects coding contrasts(size.z$species.wec.1) <- contr.wec(size.z$species.wec.1, "C.ma") #selecte C. maculata as the omitted category for weighted effects coding contrasts(size.z$species.wec.2) <- contr.wec(size.z$species.wec.2, "H.a") #selecte H. axyridis as the omitted category for weighted effects coding #create interaction term with intensive cropland size.z$crop.species.inter1 <- wec.interact(size.z$species.wec.1, size.z$crop) size.z$crop.species.inter2 <- wec.interact(size.z$species.wec.2, size.z$crop) #create interaction term with landscape diversity size.z$sidi.species.inter1 <- wec.interact(size.z$species.wec.1, size.z$sidi) size.z$sidi.species.inter2 <- wec.interact(size.z$species.wec.2, size.z$sidi) #create interaction term with edge density size.z$ed.species.inter1 <- wec.interact(size.z$species.wec.1, size.z$ed) size.z$ed.species.inter2 <- wec.interact(size.z$species.wec.2, size.z$ed) ``` ###Density ```{r} density <- droplevels(subset(data, density>0 & tenerous=="no")) #remove missing values and tenerous beetles density.z <- data.frame(crop = standardize(density$crop), sidi = standardize(density$sidi), ed=standardize(density$ed), species=density$species, species.wec.1=density$species, #for weighted effects coding species.wec.2=density$species, #for weighted effects coding cd=density$cd, siteID=density$siteID, sex=density$sex, density.z=standardize(log(density$density))) #ln-transform and standardize #weighted effects coding contrasts(density.z$species.wec.1) <- contr.wec(density.z$species.wec.1, "C.ma") #selecte C. maculata as the omitted category for weighted effects coding contrasts(density.z$species.wec.2) <- contr.wec(density.z$species.wec.2, "H.a") #selecte H. axyridis as the omitted category for weighted effects coding density.z$crop.species.inter1 <- wec.interact(density.z$species.wec.1, density.z$crop) density.z$crop.species.inter2 <- wec.interact(density.z$species.wec.2, density.z$crop) density.z$sidi.species.inter1 <- wec.interact(density.z$species.wec.1, density.z$sidi) density.z$sidi.species.inter2 <- wec.interact(density.z$species.wec.2, density.z$sidi) density.z$ed.species.inter1 <- wec.interact(density.z$species.wec.1, density.z$ed) density.z$ed.species.inter2 <- wec.interact(density.z$species.wec.2, density.z$ed) ``` ###Lipid content ```{r} lipid <- droplevels(subset(data, lipid>0)) #remove missing values lipid.z <- data.frame(crop=standardize(lipid$crop), sidi=standardize(lipid$sidi), ed=standardize(lipid$ed), species=lipid$species, species.wec.1=lipid$species, #for weighted effects coding species.wec.2=lipid$species, #for weighted effects coding cd=lipid$cd, siteID=lipid$siteID, sex=lipid$sex, lipid.z=standardize(log(100*lipid$lipid))) #convert to percent, ln-transform, standardize contrasts(lipid.z$species.wec.1) <- contr.wec(lipid.z$species.wec.1, "C.ma") #selecte C. maculata as the omitted category for weighted effects coding contrasts(lipid.z$species.wec.2) <- contr.wec(lipid.z$species.wec.2, "H.a") #selecte H. axyridis as the omitted category for weighted effects coding lipid.z$crop.species.inter1 <- wec.interact(lipid.z$species.wec.1, lipid.z$crop) lipid.z$crop.species.inter2 <- wec.interact(lipid.z$species.wec.2, lipid.z$crop) lipid.z$sidi.species.inter1 <- wec.interact(lipid.z$species.wec.1, lipid.z$sidi) lipid.z$sidi.species.inter2 <- wec.interact(lipid.z$species.wec.2, lipid.z$sidi) lipid.z$ed.species.inter1 <- wec.interact(lipid.z$species.wec.1, lipid.z$ed) lipid.z$ed.species.inter2 <- wec.interact(lipid.z$species.wec.2, lipid.z$ed) ``` ###Eggs ```{r} eggs <- droplevels(subset(data, eggs>=0 & cd>130 & tenerous=="no" & species!="H.v" & density > 0)) #remove missing values, newly emerged & tenerous beetles, H. variegata) eggs$species.wec.1 <- eggs$species #for weighted effects coding eggs$species.wec.2 <- eggs$species #for weighted effects coding #for lipid content eggs.noC7 <- droplevels(filter(eggs, species != "C.s")) #for lipid analysis (no C. 7 measured for lipids) contrasts(eggs.noC7$species.wec.1) <- contr.wec(eggs.noC7$species.wec.1, "C.ma") #for weighted effects coding contrasts(eggs.noC7$species.wec.2) <- contr.wec(eggs.noC7$species.wec.2, "H.a") #for weighted effects coding eggs.noC7$lipid.species.inter1 <- wec.interact(eggs.noC7$species.wec.1, eggs.noC7$lipid) eggs.noC7$lipid.species.inter2 <- wec.interact(eggs.noC7$species.wec.2, eggs.noC7$lipid) #for body size and density contrasts(eggs$species.wec.1) <- contr.wec(eggs$species.wec.1, "C.ma") #select C. maculata as the omitted category for weighted effects coding contrasts(eggs$species.wec.2) <- contr.wec(eggs$species.wec.2, "H.a") #select H. axyridis as the omitted category for weighted effects coding eggs$size.species.inter1 <- wec.interact(eggs$species.wec.1, eggs$size) eggs$size.species.inter2 <- wec.interact(eggs$species.wec.2, eggs$size) eggs$density.species.inter1 <- wec.interact(eggs$species.wec.1, eggs$density) eggs$density.species.inter2 <- wec.interact(eggs$species.wec.2, eggs$density) ``` ##Create models and output dataframes ###Body size ```{r} size_mods <- function(df) { mod.crop.1 <- lme(size.z ~ crop + species.wec.1 + crop.species.inter1 + cd*species.wec.1, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.crop.2 <- lme(size.z ~ crop + species.wec.2 + crop.species.inter2 + cd*species.wec.2, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.sidi.1 <- lme(size.z ~ sidi + species.wec.1 + sidi.species.inter1 + cd*species.wec.1, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.sidi.2 <- lme(size.z ~ sidi + species.wec.2 + sidi.species.inter2 + cd*species.wec.2, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.ed.1 <- lme(size.z ~ ed + species.wec.1 + ed.species.inter1 + cd*species.wec.1, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.ed.2 <- lme(size.z ~ ed + species.wec.2 + ed.species.inter2 + cd*species.wec.2, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) predictor = as.factor(c(rep.int("% Intensive cropland", 1+length(unique(df$species))), rep.int("Landscape diversity", 1+length(unique(df$species))), rep.int("Edge density", 1+length(unique(df$species))) ) ) response = as.factor(rep.int("Body size", 3*(1+length(unique(df$species)))) ) species = as.factor(rep.int(c("All species", "C. maculata", "H. axyridis", "C. munda", "C. 7-punctata", "H. variegata"), 3) ) estimate <- c( #intensive cropland intervals(mod.crop.1, which="fixed")[1]$fixed[2,2], #all species intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[7,2], #C.ma intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[7,2], #H.a intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[8,2], #C.mu intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[9,2], #C.s intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[10,2], #H.v #landscape diversity intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2], #all species intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[7,2], #C.ma intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[7,2], #H.a intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[8,2], #C.mu intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[9,2], #C.s intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[10,2], #H.v #edge density intervals(mod.ed.1, which="fixed")[1]$fixed[2,2], #all species intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[7,2], #C.ma intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[7,2], #H.a intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[8,2], #C.mu intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[9,2], #C.s intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[10,2] #H.v ) lower <- c( #intensive cropland intervals(mod.crop.1, which="fixed")[1]$fixed[2,1], #all species intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[7,1], #C.ma intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[7,1], #H.a intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[8,1], #C.mu intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[9,1], #C.s intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[10,1], #H.v #landscape diversity intervals(mod.sidi.1, which="fixed")[1]$fixed[2,1], #all species intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[7,1], #C.ma intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[7,1], #H.a intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[8,1], #C.mu intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[9,1], #C.s intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[10,1], #H.v #edge density intervals(mod.ed.1, which="fixed")[1]$fixed[2,1], #all species intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[7,1], #C.ma intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[7,1], #H.a intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[8,1], #C.mu intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[9,1], #C.s intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[10,1] #H.v ) upper <- c( #intensive cropland intervals(mod.crop.1, which="fixed")[1]$fixed[2,3], #all species intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[7,3], #C.ma intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[7,3], #H.a intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[8,3], #C.mu intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[9,3], #C.s intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[10,3], #H.v #landscape diversity intervals(mod.sidi.1, which="fixed")[1]$fixed[2,3], #all species intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[7,3], #C.ma intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[7,3], #H.a intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[8,3], #C.mu intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[9,3], #C.s intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[10,3], #H.v #edge density intervals(mod.ed.1, which="fixed")[1]$fixed[2,3], #all species intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[7,3], #C.ma intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[7,3], #H.a intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[8,3], #C.mu intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[9,3], #C.s intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[10,3] #H.v ) #summary table summary.table <-bind_rows(as.data.frame(summary(mod.crop.2)[20]$tTable[c(2,7),]), as.data.frame(summary(mod.crop.1)[20]$tTable[c(7:10),]), as.data.frame(summary(mod.sidi.2)[20]$tTable[c(2,7),]), as.data.frame(summary(mod.sidi.1)[20]$tTable[c(7:10),]), as.data.frame(summary(mod.ed.2)[20]$tTable[c(2,7),]), as.data.frame(summary(mod.ed.1)[20]$tTable[c(7:10),])) summary.table_size <<- as_tibble(summary.table) %>% mutate(response=as.factor(rep.int("Body size", nrow(summary.table))), variable=row.names(summary.table), .before=Value ) #ANOVA table anova.table <- bind_rows(Anova(mod.crop.1), Anova(mod.sidi.1), Anova(mod.ed.1)) anova.table_size <<- as_tibble(anova.table) %>% mutate(response= as.factor(rep.int("Body size", nrow(Anova(mod.crop.1))+nrow(Anova(mod.sidi.1))+nrow(Anova(mod.ed.1)))), predictor=as.factor(c(rep.int("% Intensive cropland", nrow(Anova(mod.crop.1))), rep.int("Landscape diversity", nrow(Anova(mod.sidi.1))), rep.int("Edge density", nrow(Anova(mod.ed.1)))) ), variable= row.names(anova.table), .before = Chisq ) figure.df <- data.frame(predictor, response, species, estimate, lower, upper) return(figure.df) } ``` ###Density ```{r} density_mods <- function(df) { mod.crop.1 <- lme(density.z ~ crop + species.wec.1 + crop.species.inter1 + poly(cd,2)*species.wec.1, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.crop.2 <- lme(density.z ~ crop + species.wec.2 + crop.species.inter2 + poly(cd,2)*species.wec.2, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.sidi.1 <- lme(density.z ~ sidi + species.wec.1 + sidi.species.inter1 + poly(cd,2)*species.wec.1, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.sidi.2 <- lme(density.z ~ sidi + species.wec.2 + sidi.species.inter2 + poly(cd,2)*species.wec.2, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.ed.1 <- lme(density.z ~ ed + species.wec.1 + ed.species.inter1 + poly(cd,2)*species.wec.1, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.ed.2 <- lme(density.z ~ ed + species.wec.2 + ed.species.inter2 + poly(cd,2)*species.wec.2, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) predictor = as.factor(c(rep.int("% Intensive cropland", 1+length(unique(df$species))), rep.int("Landscape diversity", 1+length(unique(df$species))), rep.int("Edge density", 1+length(unique(df$species))) ) ) response = as.factor(rep.int("Density", 3*(1+length(unique(df$species)))) ) species = as.factor(rep.int(c("All species", "C. maculata", "H. axyridis", "C. munda", "C. 7-punctata", "H. variegata"), 3) ) estimate <- c( #intensive cropland intervals(mod.crop.1, which="fixed")[1]$fixed[2,2], #all species intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[7,2], #C.ma intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[7,2], #H.a intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[8,2], #C.mu intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[9,2], #C.s intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[10,2], #H.v #landscape diversity intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2], #all species intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[7,2], #C.ma intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[7,2], #H.a intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[8,2], #C.mu intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[9,2], #C.s intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[10,2], #H.v #edge density intervals(mod.ed.1, which="fixed")[1]$fixed[2,2], #all species intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[7,2], #C.ma intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[7,2], #H.a intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[8,2], #C.mu intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[9,2], #C.s intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[10,2] #H.v ) lower <- c( #intensive cropland intervals(mod.crop.1, which="fixed")[1]$fixed[2,1], #all species intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[7,1], #C.ma intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[7,1], #H.a intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[8,1], #C.mu intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[9,1], #C.s intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[10,1], #H.v #landscape diversity intervals(mod.sidi.1, which="fixed")[1]$fixed[2,1], #all species intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[7,1], #C.ma intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[7,1], #H.a intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[8,1], #C.mu intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[9,1], #C.s intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[10,1], #H.v #edge density intervals(mod.ed.1, which="fixed")[1]$fixed[2,1], #all species intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[7,1], #C.ma intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[7,1], #H.a intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[8,1], #C.mu intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[9,1], #C.s intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[10,1] #H.v ) upper <- c( #intensive cropland intervals(mod.crop.1, which="fixed")[1]$fixed[2,3], #all species intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[7,3], #C.ma intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[7,3], #H.a intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[8,3], #C.mu intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[9,3], #C.s intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[10,3], #H.v #landscape diversity intervals(mod.sidi.1, which="fixed")[1]$fixed[2,3], #all species intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[7,3], #C.ma intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[7,3], #H.a intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[8,3], #C.mu intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[9,3], #C.s intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[10,3], #H.v #edge density intervals(mod.ed.1, which="fixed")[1]$fixed[2,3], #all species intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[7,3], #C.ma intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[7,3], #H.a intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[8,3], #C.mu intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[9,3], #C.s intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[10,3] #H.v ) #summary statistics summary.table <-bind_rows(as.data.frame(summary(mod.crop.2)[20]$tTable[c(2,7),]), as.data.frame(summary(mod.crop.1)[20]$tTable[c(7:10),]), as.data.frame(summary(mod.sidi.2)[20]$tTable[c(2,7),]), as.data.frame(summary(mod.sidi.1)[20]$tTable[c(7:10),]), as.data.frame(summary(mod.ed.2)[20]$tTable[c(2,7),]), as.data.frame(summary(mod.ed.1)[20]$tTable[c(7:10),])) summary.table_density <<- as_tibble(summary.table) %>% mutate(response=as.factor(rep.int("Body density", nrow(summary.table))), variable=row.names(summary.table), .before=Value ) anova.table <- bind_rows(Anova(mod.crop.1), Anova(mod.sidi.1), Anova(mod.ed.1)) anova.table_density <<- as_tibble(anova.table) %>% mutate(response= as.factor(rep.int("Body density", nrow(Anova(mod.crop.1))+nrow(Anova(mod.sidi.1))+nrow(Anova(mod.ed.1)))), predictor=as.factor(c(rep.int("% Intensive cropland", nrow(Anova(mod.crop.1))), rep.int("Landscape diversity", nrow(Anova(mod.sidi.1))), rep.int("Edge density", nrow(Anova(mod.ed.1)))) ), variable= row.names(anova.table), .before = Chisq ) figure.df <- data.frame(predictor, response, species, estimate, lower, upper) return(figure.df) } ``` ###Lipid content ```{r} lipid_mods <- function(df) { mod.crop.1 <- lme(lipid.z ~ crop + species.wec.1 + crop.species.inter1 + poly(cd,2)*species.wec.1, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.crop.2 <- lme(lipid.z ~ crop + species.wec.2 + crop.species.inter2 + poly(cd,2)*species.wec.2, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.sidi.1 <- lme(lipid.z ~ sidi + species.wec.1 + sidi.species.inter1 + poly(cd,2)*species.wec.1, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.sidi.2 <- lme(lipid.z ~ sidi + species.wec.2 + sidi.species.inter2 + poly(cd,2)*species.wec.2, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.ed.1 <- lme(lipid.z ~ ed + species.wec.1 + ed.species.inter1 + poly(cd,2)*species.wec.1, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) mod.ed.2 <- lme(lipid.z ~ ed + species.wec.2 + ed.species.inter2 + poly(cd,2)*species.wec.2, random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=df) predictor = as.factor(c(rep.int("% Intensive cropland", 1+length(unique(df$species))), rep.int("Landscape diversity", 1+length(unique(df$species))), rep.int("Edge density", 1+length(unique(df$species))) ) ) response = as.factor(rep.int("Lipid content", 3*(1+length(unique(df$species)))) ) species = as.factor(rep.int(c("All species", "C. maculata", "H. axyridis", "C. munda", "H. variegata"), 3) ) estimate <- c( #intensive cropland intervals(mod.crop.1, which="fixed")[1]$fixed[2,2], #all species intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[6,2], #C.ma intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[6,2], #H.a intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[7,2], #C.mu intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[8,2], #H.v #landscape diversity intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2], #all species intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[6,2], #C.ma intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[6,2], #H.a intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[7,2], #C.mu intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[8,2], #H.v #edge density intervals(mod.ed.1, which="fixed")[1]$fixed[2,2], #all species intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[6,2], #C.ma intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[6,2], #H.a intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[7,2], #C.mu intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[8,2] #H.v ) lower <- c( #intensive cropland intervals(mod.crop.1, which="fixed")[1]$fixed[2,1], #all species intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[6,1], #C.ma intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[6,1], #H.a intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[7,1], #C.mu intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[8,1], #H.v #landscape diversity intervals(mod.sidi.1, which="fixed")[1]$fixed[2,1], #all species intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[6,1], #C.ma intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[6,1], #H.a intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[7,1], #C.mu intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[8,1], #H.v #edge density intervals(mod.ed.1, which="fixed")[1]$fixed[2,1], #all species intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[6,1], #C.ma intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[6,1], #H.a intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[7,1], #C.mu intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[8,1] #H.v ) upper <- c( #intensive cropland intervals(mod.crop.1, which="fixed")[1]$fixed[2,3], #all species intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[6,3], #C.ma intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[6,3], #H.a intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.1, which="fixed")[1]$fixed[7,3], #C.mu intervals(mod.crop.1, which="fixed")[1]$fixed[2,2] + intervals(mod.crop.2, which="fixed")[1]$fixed[8,3], #H.v #landscape diversity intervals(mod.sidi.1, which="fixed")[1]$fixed[2,3], #all species intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[6,3], #C.ma intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[6,3], #H.a intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.1, which="fixed")[1]$fixed[7,3], #C.mu intervals(mod.sidi.1, which="fixed")[1]$fixed[2,2] + intervals(mod.sidi.2, which="fixed")[1]$fixed[8,3], #H.v #edge density intervals(mod.ed.1, which="fixed")[1]$fixed[2,3], #all species intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[6,3], #C.ma intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[6,3], #H.a intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.1, which="fixed")[1]$fixed[7,3], #C.mu intervals(mod.ed.1, which="fixed")[1]$fixed[2,2] + intervals(mod.ed.2, which="fixed")[1]$fixed[8,3] #H.v ) #summary statistics summary.table <-bind_rows(as.data.frame(summary(mod.crop.2)[20]$tTable[c(2,6),]), as.data.frame(summary(mod.crop.1)[20]$tTable[c(6:8),]), as.data.frame(summary(mod.sidi.2)[20]$tTable[c(2,6),]), as.data.frame(summary(mod.sidi.1)[20]$tTable[c(6:8),]), as.data.frame(summary(mod.ed.2)[20]$tTable[c(2,6),]), as.data.frame(summary(mod.ed.1)[20]$tTable[c(6:8),])) summary.table_lipid <<- as_tibble(summary.table) %>% mutate(response=as.factor(rep.int("Lipid content", nrow(summary.table))), variable=row.names(summary.table), .before=Value ) anova.table <- bind_rows(Anova(mod.crop.1), Anova(mod.sidi.1), Anova(mod.ed.1)) anova.table_lipid <<- as_tibble(anova.table) %>% mutate(response= as.factor(rep.int("Lipid content", nrow(Anova(mod.crop.1))+nrow(Anova(mod.sidi.1))+nrow(Anova(mod.ed.1)))), predictor=as.factor(c(rep.int("% Intensive cropland", nrow(Anova(mod.crop.1))), rep.int("Landscape diversity", nrow(Anova(mod.sidi.1))), rep.int("Edge density", nrow(Anova(mod.ed.1)))) ), variable= row.names(anova.table), .before = Chisq ) figure.df <- data.frame(predictor, response, species, estimate, lower, upper) return(figure.df) } ``` ###Eggs ```{r} #Body size model eggs.size=glmmTMB(eggs ~ size*species + cd*species + (1 |siteID/species), family=binomial, data = eggs) #weighted effects coding eggs.size.1=glmmTMB(eggs ~ size + species.wec.1 + size.species.inter1 + cd*species.wec.1 + (1 |siteID/species), family=binomial, data = eggs) eggs.size.2=glmmTMB(eggs ~ size + species.wec.2 + size.species.inter2 + cd*species.wec.2 + (1 |siteID/species), family=binomial, data = eggs) #Body density model eggs.density=glmmTMB(eggs ~ density*species + cd*species + (1 |siteID/species), family=binomial, data = eggs) #weighted effects coding eggs.density.1=glmmTMB(eggs ~ density + species.wec.1 + density.species.inter1 + cd*species.wec.1 + (1 |siteID/species), family=binomial, data = eggs) eggs.density.2=glmmTMB(eggs ~ density + species.wec.2 + density.species.inter2 + cd*species.wec.2 + (1 |siteID/species), family=binomial, data = eggs) #Lipid content model eggs.lipid=glmmTMB(eggs ~ 100*lipid*species + cd*species + (1 |siteID/species), family=binomial, data = eggs.noC7) #weighted effects coding eggs.lipid.1=glmmTMB(eggs ~ 100*lipid + species.wec.1 + lipid.species.inter1 + cd*species.wec.1 + (1 |siteID/species), family=binomial, data = eggs.noC7) eggs.lipid.2=glmmTMB(eggs ~ 100*lipid + species.wec.2 + lipid.species.inter2 + cd*species.wec.2 + (1 |siteID/species), family=binomial, data = eggs.noC7) ``` ##Results ```{r} size.results <- size_mods(size.z) density.results <- density_mods(density.z) lipid.results <- lipid_mods(lipid.z) anova.table <- bind_rows(anova.table_size, anova.table_density, anova.table_lipid) summary.table <- bind_rows(summary.table_size, summary.table_density, summary.table_lipid) #create dataframe for Figure 2 all.results.df <- bind_rows(size.results, density.results, lipid.results) all.results.df$significant <- as.factor(ifelse(all.results.df$lower < 0 & all.results.df$upper > 0, "n", "y")) #significant if 95% CI doesn't cross 0 all.results.df$point <- as.factor(ifelse(all.results.df$species == "All species", "big", "small")) #distinguish points for overall effects vs. species-specific responses all.results.df$point <- factor(all.results.df$point, levels=c("small", "big")) #reorder levels #order factor levels all.results.df$predictor <- factor(all.results.df$predictor, levels=c("% Intensive cropland", "Landscape diversity", "Edge density")) all.results.df$response <- factor(all.results.df$response, levels=c("Lipid content", "Density", "Body size")) all.results.df$species <- factor(all.results.df$species, levels=c("All species", "C. maculata", "H. axyridis", "C. munda", "C. 7-punctata", "H. variegata")) ``` ##Figures ###Figure 2 ```{r} Fig2 <- ggplot(data=all.results.df, aes(x=response, y=estimate, ymin=-1, ymax=1)) + facet_wrap(~predictor, strip.position = "left", nrow=3, scales="free_y") + geom_hline(aes(fill=response), yintercept = 0, linetype=2) + xlab("") + ylab("Standardized effect size") + coord_flip() + theme_minimal() + theme(strip.placement = "outside", text = element_text(size=18)) + labs(col="Species") + geom_point(aes(col = fct_rev(species), alpha = significant, size=point, shape=point), position=position_dodge(width=.7)) + geom_errorbar(aes(ymin=lower, ymax=upper, col=fct_rev(species), alpha = significant), width=0, cex=1.5, position=position_dodge(width=.7)) + guides(color=guide_legend(reverse = T)) + scale_color_manual(values = myPalette2.rev, labels = c("*H. variegata*", "*C. 7-punctata*","*C. munda*", "*H. axyridis*", "*C. maculata*", "All species")) + scale_alpha_discrete(range=c(.2,1), guide="none") + scale_size_discrete(range=c(3,7), guide="none") + scale_shape_manual(values=c(16,18), guide="none") + theme(legend.text = element_markdown(), axis.title = element_text(size=12), axis.text = element_text(size=12), strip.text = element_text(size=12)) ``` ###Figure 3 ####A. Cropland vs. body density ```{r} density.crop.mod <- lme(log(density)~crop*species + poly(cd,2)*species,random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=density) density.crop.fit <- ggeffect(density.crop.mod, terms=c("crop [all]", "species")) density.crop.df <- data.frame(cropland=density$crop, density=(resid(density.crop.mod)+predict(density.crop.mod)), species=density$species) density.crop.plot <- ggplot(density.crop.df, aes(x=cropland, y=exp(density), color=species)) + geom_point(alpha=0.15) + geom_line(data=density.crop.fit, aes(x=x, y=exp(predicted), color=group), size=2) + geom_ribbon(data=density.crop.fit, aes(x=x, y=exp(predicted), ymin=exp(conf.low), ymax=exp(conf.high), fill=group), alpha=0.1, color=NA, show.legend = F) + xlab("% Intensive cropland") + ylab(expression(bold(Body~density~(mg/mm^"2")))) + labs(color="Species") + theme_classic() + scale_color_manual(values=myPalette, labels=c("C. maculata", "H. axyridis", "C. munda", "C. 7-punctata", "H. variegata")) + scale_fill_manual(values=myPalette) + theme(axis.title = element_text(size=12, face='bold'), axis.text = element_text(size=12, face = 'bold'), legend.text = element_text(size=16, face = 'bold'), legend.title = element_text(size=20, face = 'bold')) ``` ####B. Edge density vs. body density ```{r} density.edge.mod <- lme(log(density)~ed*species + poly(cd,2)*species,random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=density) density.edge.fit <- ggeffect(density.edge.mod, terms=c("ed [all]", "species")) density.edge.df <- data.frame(edge=density$ed, density=(resid(density.edge.mod)+predict(density.edge.mod)), species=density$species) density.edge.plot <- ggplot(density.edge.df, aes(x=edge, y=exp(density), color=species)) + geom_point(alpha=0.15) + geom_line(data=density.edge.fit, aes(x=x, y=exp(predicted), color=group), size=2) + geom_ribbon(data=density.edge.fit, aes(x=x, y=exp(predicted), ymin=exp(conf.low), ymax=exp(conf.high), fill=group), alpha=0.1, color=NA, show.legend = F) + xlab("Edge density (m/hectare)") + ylab(expression(bold(Body~density~(mg/mm^"2")))) + labs(color="Species") + theme_classic() + scale_color_manual(values=myPalette, labels=c("C. maculata", "H. axyridis", "C. munda", "C. 7-punctata", "H. variegata")) + scale_fill_manual(values=myPalette) + theme(axis.title = element_text(size=12, face='bold'), axis.text = element_text(size=12, face = 'bold'), legend.position = "none") ``` ####C. Cropland vs. lipid content ```{r} lipid.crop.mod <- lme(log(100*lipid)~crop*species + poly(cd,2)*species,random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=lipid) lipid.crop.fit <- ggeffect(lipid.crop.mod, terms=c("crop [all]", "species")) lipid.crop.df <- data.frame(crop=lipid$crop, lipid=(resid(lipid.crop.mod)+predict(lipid.crop.mod)), species=lipid$species) lipid.crop.plot <- ggplot(lipid.crop.df, aes(x=crop, y=exp(lipid), color=species)) + geom_point(alpha=0.15) + geom_line(data=lipid.crop.fit, aes(x=x, y=exp(predicted), color=group), size=2) + geom_ribbon(data=lipid.crop.fit, aes(x=x, y=exp(predicted), ymin=exp(conf.low), ymax=exp(conf.high), fill=group), alpha=0.1, color=NA, show.legend = F) + xlab("% Intensive cropland") + ylab("Lipid content (%)") + labs(color="Species") + theme_classic() + scale_color_manual(values=myPalette.noC7) + scale_fill_manual(values=myPalette, labels=c("C. maculata", "H. axyridis", "C. munda", "C. 7-punctata", "H. variegata")) + theme(axis.title = element_text(size=12, face='bold'), axis.text = element_text(size=12, face = 'bold'), legend.position = "none") ``` ####D. Landscape diversity vs. lipid content ```{r} lipid.sidi.mod <- lme(log(100*lipid)~sidi*species + poly(cd,2)*species,random=~1|siteID/species/sex, weights = varIdent(form = ~1 |species), data=lipid) lipid.sidi.fit <- ggeffect(lipid.sidi.mod, terms=c("sidi [all]", "species")) lipid.sidi.df <- data.frame(sidi=lipid$sidi, lipid=(resid(lipid.sidi.mod)+predict(lipid.sidi.mod)), species=lipid$species) lipid.sidi.plot <- ggplot(lipid.sidi.df, aes(x=sidi, y=exp(lipid), color=species)) + geom_point(alpha=0.15) + geom_line(data=lipid.sidi.fit, aes(x=x, y=exp(predicted), color=group), size=2) + geom_ribbon(data=lipid.sidi.fit, aes(x=x, y=exp(predicted), ymin=exp(conf.low), ymax=exp(conf.high), fill=group), alpha=0.1, color=NA, show.legend = F) + xlab("Simpson's diversity index") + ylab("Lipid content (%)") + labs(color="Species") + theme_classic() + scale_color_manual(values=myPalette.noC7, labels=c("C. maculata", "H. axyridis", "C. munda", "H. variegata")) + scale_fill_manual(values=myPalette) + theme(axis.title = element_text(size=12, face='bold'), axis.text = element_text(size=12, face = 'bold'), legend.position = "none") ``` ####Combined figure ```{r} Fig3 <- ggarrange(density.crop.plot, density.edge.plot, lipid.crop.plot, lipid.sidi.plot, labels = c("A", "B", "C"), ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom") ``` #Suplementary Material ##Figures ###Figure S1 ```{r} body.vars <- data.frame(size=data$size, density=data$density, lipid=100*data$lipid, species=factor(data$species, levels = c("C.ma", "H.a", "C.mu", "C.s", "H.v"))) #create dataframe body.vars[complete.cases(body.vars[ , 1:2]),] #remove missing values names(body.vars) <- c("Body size", "Body density", "Lipid content", "Species") levels(body.vars$Species) <- c("C. maculata", "H. axyridis", "C. munda", "C. 7-punctata", "H. variegata") body.plot <- ggpairs(body.vars, lower = list(continuous = wrap("points", alpha=0.2)), upper = list(continuous = wrap("cor", size = 5)), diag = list(continuous = wrap("densityDiag", alpha=0.5)), ggplot2::aes(color=Species)) + theme(text = element_text(size=12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), strip.text = element_text(size = 14), strip.background = element_blank(), axis.text.x = element_text(angle = 60, hjust = 1)) #change color palette for(i in 1:body.plot$nrow) { for(j in 1:body.plot$ncol){ body.plot[i,j] <- body.plot[i,j] + scale_fill_manual(values=myPalette) + scale_color_manual(values=myPalette) } } #adjust color palette for panels missing data for C. 7-punctata body.plot[3,3] <- body.plot[3,3] + scale_color_manual(values=myPalette.noC7) + scale_fill_manual(values=myPalette.noC7) body.plot[4,3] <- body.plot[4,3] + scale_color_manual(values=myPalette.noC7) + scale_fill_manual(values=myPalette.noC7) body.plot[3,4] <- body.plot[3,4] + scale_color_manual(values=myPalette.noC7) + scale_fill_manual(values=myPalette.noC7) body.plot[1,3] <- body.plot[1,3] + scale_color_manual(values=myPalette.noC7) + scale_fill_manual(values=myPalette.noC7) body.plot[2,3] <- body.plot[2,3] + scale_color_manual(values=myPalette.noC7) + scale_fill_manual(values=myPalette.noC7) #adjust y axis on histograms body.plot[4,1] <- body.plot[4,1] + scale_y_continuous(breaks = c(0, 50, 100)) body.plot[4,2] <- body.plot[4,2] + scale_y_continuous(breaks = c(0, 50, 100)) body.plot[4,3] <- body.plot[4,3] + scale_y_continuous(breaks = c(0, 50, 100)) #add counts to bar chart body.plot[4,4] <- body.plot[4,4] + geom_text(stat='count', aes(label=..count..), vjust=-0.1) ``` ###Figure S2 ```{r} landscape.vars <- landscape.data[,c(3:5)] names(landscape.vars) <- c("% Intensive crop", "Landscape diversity", "Edge density") landscape.plot <- ggpairs(landscape.vars, upper = list(continuous = wrap("cor", size = 8))) + theme(text = element_text(size=15), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), strip.text = element_text(size = 16), strip.background = element_blank(), axis.text.x = element_text(angle = 60, hjust = 1)) #remove axis from top right landscape.plot[1,1] <- landscape.plot[1,1] + theme(axis.text.y = element_blank(), axis.ticks = element_blank()) ``` ###Figure S3 ```{r} #create density model eggs.density.fit <- ggeffect(eggs.density, terms=c("density [all]", "species")) levels(eggs.density.fit$group) <- c("C. maculata", "H. axyridis", "C. munda", "C. 7-punctata") eggs.density.df <- data.frame(density=eggs$density, eggs=eggs$eggs, group=eggs$species) levels(eggs.density.df$group) <- c("C. maculata", "H. axyridis", "C. munda", "C. 7-punctata") FigS3 <- ggplot(eggs.density.df, aes(x=density, y=eggs, color=group)) + geom_jitter(alpha=0.25, height = 0.01) + geom_smooth(method = "glm", method.args = list(family = "binomial"), se=FALSE, data=eggs.density.fit, aes(x=x, y=predicted, color=group), size=2) + geom_ribbon(data=eggs.density.fit, aes(x=x, y=predicted, ymin=conf.low, ymax=conf.high, fill=group), alpha=0.1, color=NA, show.legend = F) + facet_wrap(~group, ncol = 2) + xlab(expression(bold(Body~density~(mg/mm^"2")))) + ylab("Probability of egg presence") + labs(color="Species") + theme_classic() + scale_color_manual(values=myPalette) + scale_fill_manual(values=myPalette) + theme(axis.title = element_text(size=12, face='bold'), axis.text = element_text(size=12, face = 'bold'), legend.position = "none", strip.background = element_blank(), strip.text = element_text(face = "bold", size = 12)) ```