knitr::opts_chunk$set(echo = TRUE, out.width = '100%')
setwd("~/PhD")

Description of research

This study seeks to address how strongly laboratory leaf consumption assays with a generalist herbivore predict naturally occurring arthropod counts, mass and densities on native plants. The metadata used for our analyses was collected each year in April-May from 2021-2023 in a common garden at Concordia University Irvine planted with woody shrubs native to Southern California. The common garden is made up of 20 separate plots measuring 6m x 5m, each with 1 individual plant per species. Plants were grown from seed in a greenhouse and transplanted after 1 year. Data collected, abbreviations used and variable type assigned for the metadata can be found in the table below.

Data Description Column Name Variable Type
Plant # Individual plant Plant Factor
Plant Species Plant Species Species Factor
Watering Treatment Water added to plots ranging from Ambient (1) to High (6) Watering_Treatment Factor
Plot Individual Plot Plot Factor
Year Year data collected Year Factor
Spodoptera percent survival % Spodoptera surviving at end of trial Spod_perc_survival Numeric
Average Spodoptera mass (mg) Average mass of surviving Spodoptera Avg_spod_mass_mg Numeric
Individual Spodoptera mass (mg) Individual mass of surviving Spodoptera Spod_mass_mg Numeric
Individual Spodopgera survial Spodoptera alive/dead at end of trial Spod_survival Factor
Total arthropod count Total number of arthropods Abundance_Total Numeric
Herbivore count Total number of herbivores Abundance_Herb Numeric
Predator count Total number of Predators Abundance_Pred Numeric
Omnivore count Total number of Omnivores Abundance_Omni Numeric
Detritivore count Total number of Detritivores Abundance_Detr Numeric
Palynivore count Total number of Palynivores Abundance_Paly Numeric
Total arthropod mass (mg) Combined mass of all arthropods Arth_Mass_Total_mg Numeric
Herbivore mass (mg) Combined mass of all herbivbores Mass_Herb_mg Numeric
Predator mass (mg) Combined mass of all predators Mass_Pred_mg Numeric
Omnivore mass (mg) Combined mass of all omnivores Mass_Omni_mg Numeric
Detritivore mass (mg) Combined mass of all detritivores Mass_Detr_mg Numeric
Palynivore mass (mg) Combined mass of all palynivores Mass_Paly_mg Numeric
Estimated dry plant biomass (g) Estimated dry biomass of aboveground portion of plant Dry_Plant_Biomass_g Numeric
Estimated plant biomass sampled (g) Estimated dry biomass of plant sampled during arthropod collection Dry_Biomass_Sampled_g Numeric
Cartesian coordinates of plot X,Y coordinates of plot Plot_X/Y Factor
Change in plant mass Percent change from year to year in estimted dry plant biomass Plant_mass_perc_change Numeric
Plant survival Plant alive/dead at annual sampling Plant_surv Factor

Methods

Spodoptera data was collected from a 10-day feeding assay where Spodoptera exigua, a generalist herbivore, were hatched from eggs and immediately fed leaves from each individual plant. At the end of 10 days, the number of surviving larvae and final mass of each larvae were recorded. See You, An, and Li (2020) for details. Plant and arthropod data were collected from each surviving plant from 2021-2023. Plant estimated biomass was collected using a branch-biomass estimation method as previously described in Nell and Mooney (2019). Arthropods were collected from plants using a vaccum method as previously described in Pratt et al. (2017), sorted from plant material and identified to feeding guild. For 2023 arthropod samples, lengths of arthropods were measured and used to estimate arthropod dry biomass. See Rogers, Hinds, and Buschbom (1976) for details.

1- Load packages and data tables

1.1 Packages

library(dplyr)
library(lme4)
library(RColorBrewer)
library(devtools)
library(sjPlot)
library(tidyr)
library(broom.mixed)
library(knitr)
library(broom)
library(jtools)
library(kableExtra)
library(huxtable)
library(dplyr)
library(gtsummary)
library(car)
library(lmerTest)
library(ggplot2)
library(sjmisc)
library(sjstats)
library(sjlabelled)
library(glmmTMB)
library(DHARMa)
library(viridisLite)
library(data.table)
library(effects)
library(emmeans)
library(kableExtra)
library(MuMIn)
library(performance)
library(gridExtra)

1.2 Data tables

Data frames are from arthropod, plant biomass and Spodoptera bioassays from 2021-2023. Several merged data frames were created:

  1. arthorder_spodavg - Arthropod counts and biomass classified by taxonomic ORDER with AVG spododoptera mass and percent survival
  2. arthorder_spodindiv - Arthropod counts and biomass classified by taxonomic ORDER with INDIVIDUAL spodoptera mass and survival
  3. arthguild_spodavg- Arthropod counts and biomass by feeding GUILD with AVG spododoptera mass and percent survival
  4. arthguild_spodindiv - Arthropod counts and biomass by feeding guild with INDIVIDUAL spodoptera mass and survival

2- Visualize raw data

2.1 Histogram

Plot histograms of Herbivore abundance, Herbivore mass, Spodoptera survival and Avg. Spodoptera mass to visualize distribution.

TAKEAWAY- Strongly skewed distributions. Consider the need for transformation.

2.2 Linear regressions

Herbivore abundance ~ Spodoptera mass

TAKEAWAY- Seems like a neutral to positive relationship between average spodoptera mass and herbivore abundance.

Herbivore abundance ~ Spodoptera Percent Survival

TAKEAWAY Mixed bag with some showing slightly positive, some showing slightly negative, some neutral patterns.

3- Perform Linear Mixed Effects Regression analysis

I will be using the general formula

glmmTMB(arthropod ~ bioassay * species * year + (1|Watering_treatment) + (1|Plot) + (1|Plant), data = total_plant)

for all analyses with arthropod being count or mass of numbers within each feeding guild or total arthropods on a plant and bioassay being either the average spodoptera mass or number/percent of spodoptera surviving at the end of the trial. Plot, Watering Treatment and Plant will be added as random factors to the model.

These results will help to determine the strength by which feeding assay data is predictive of naturally occurring arthropod communities.

The first set of analyses will be on the general equation:

Herbivore count ~ Avg spodoptera mass

I Will use the following transformations on the predictor and response variables:

The second set of analyses will be on the general equation:

Herbivore count ~ Percent spoptera survival

I will use the following transformations on the predictor and response variables:

3.1 log Herbivore count ~ log Avg spod biomass

ANOVA results of log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq)
log(Avg_spod_mass_mg + 1) 1.874393 1 0.1709728
Species 230.350532 13 0.0000000
Year 148.281133 2 0.0000000
log(Avg_spod_mass_mg + 1):Species 33.448593 13 0.0014587
log(Avg_spod_mass_mg + 1):Year 1.550634 2 0.4605578
Species:Year 200.506749 23 0.0000000
log(Avg_spod_mass_mg + 1):Species:Year 16.796649 18 0.5371254

NOTE: The two residual plots show the following: Panel A- uniform QQ plot with a 1-1 line as predicted if distribution is normal. Panel B- Visual homogeneity of the residuals in both vertical and horizontal directions with quantile test line overlay. Assumption is that there will be no pattern as shown by quantile lines.

TAKEAWAY- Residulas look pretty good. Significant effects are:

  • Species
  • Year
  • Spod mass x species
  • Species x year

3.2 log Herbivore count ~ Avg spod biomass

ANOVA results of log Herbivore count ~ Avg Spod mass
Chisq Df Pr(>Chisq)
Avg_spod_mass_mg 2.676319 1 0.1018509
Species 232.236056 13 0.0000000
Year 156.732061 2 0.0000000
Avg_spod_mass_mg:Species 24.607695 13 0.0259789
Avg_spod_mass_mg:Year 4.783648 2 0.0914627
Species:Year 179.830028 23 0.0000000
Avg_spod_mass_mg:Species:Year 20.128666 18 0.3256226

TAKEAWAY Similar results to log transformed Avg spodoptera mass.

3.3 Herbivore count (poisson) ~ Avg spod biomass

TAKEAWAY Residuals look ok on the QQ plot, however for rank-transformed model predictions, residuals deviate signifanctly from predicted. Furthermore, ANOVA results look real funky with probabilities >Chisq being 0 for most values. This does not seem right. I don’t think this transformation will be appropriate to use.

3.4 Herbivore count (poisson) ~ log Avg spod biomass

ANOVA results of Herbivore count (poisson) ~ log Avg Spod mass
Chisq Df Pr(>Chisq)
log(Avg_spod_mass_mg + 1) 11.47389 1 0.0007058
Species 80.08979 13 0.0000000
Year 381.67255 2 0.0000000
log(Avg_spod_mass_mg + 1):Species 144.70192 13 0.0000000
log(Avg_spod_mass_mg + 1):Year 10.51086 2 0.0052191
Species:Year 846.28406 23 0.0000000
log(Avg_spod_mass_mg + 1):Species:Year 66.77889 18 0.0000002

TAKEAWAY Similar to other poisson transformation on Herbivore count- QQ plot looks good but residuals show deviation from

3.5 Herbivore count (negative binomial) ~ Avg spod biomass

ANOVA results of Herbivore count (negative binomial) ~ Avg Spod mass
Chisq Df Pr(>Chisq)
Avg_spod_mass_mg 2.222755 1 0.1359902
Species 170.239507 13 0.0000000
Year 107.434231 2 0.0000000
Avg_spod_mass_mg:Species 17.249404 13 0.1881522
Avg_spod_mass_mg:Year 2.942055 2 0.2296893
Species:Year 199.104932 23 0.0000000
Avg_spod_mass_mg:Species:Year 21.417490 18 0.2588756

TAKEAWAY Results are similar to 3.1- log herbivore ~ log avg spodoptera biomass with significance in the same categories and residuals showing no deviation from normal distribution both plots.

3.6 Herbivore count (negative binomial) ~ log Avg spod biomass

ANOVA results of Herbivore count (negative binomial) ~ log Avg Spod mass
Chisq Df Pr(>Chisq)
log(Avg_spod_mass_mg + 1) 0.4818349 1 0.4875924
Species 162.2162114 13 0.0000000
Year 91.0724496 2 0.0000000
log(Avg_spod_mass_mg + 1):Species 25.7316789 13 0.0184668
log(Avg_spod_mass_mg + 1):Year 0.9284642 2 0.6286176
Species:Year 227.3617113 23 0.0000000
log(Avg_spod_mass_mg + 1):Species:Year 18.3417522 18 0.4333617

TAKEAWAY Same results as non-log transformed avg spod mass. It seems that we can use either log transformation of both count and spodoptera mass or herbivore count (negative binomial) ~ log and non-log transformed avg spodoptera mass.

3.7 log Herbivore count ~ Percent Spodoptera survival

ANOVA results of log Herbivore count ~ Percent spodoptera survival
Chisq Df Pr(>Chisq)
Spod_perc_survival 2.666787 1 0.1024627
Species 419.955824 13 0.0000000
Year 270.294526 2 0.0000000
Spod_perc_survival:Species 17.648823 13 0.1713143
Spod_perc_survival:Year 3.625561 2 0.1631997
Species:Year 297.120638 26 0.0000000
Spod_perc_survival:Species:Year 51.545609 22 0.0003611

TAKEAWAY Residuals look ok, some deviation from expected values. Significant effects for the following:

  • Spod percent survival
  • Species
  • Year
  • Species x year
  • Survial x species x year

Does seem to be a slight effect of survival alone on herbivore abundance, stronger species and year only effects and species x year effects. Moderate species x survival x year effect. NO survival x species effect or Survival x year.

3.8 Herbivore count (poisson) ~ Percent Spodoptera survival

ANOVA results of Herbivore count (poisson) ~ Percent spodoptera survival
Chisq Df Pr(>Chisq)
Spod_perc_survival 4.17443 1 0.0410384
Species 261.03216 13 0.0000000
Year 3562.32246 2 0.0000000
Spod_perc_survival:Species 287.34118 13 0.0000000
Spod_perc_survival:Year 30.73284 2 0.0000002
Species:Year 5615.42196 26 0.0000000
Spod_perc_survival:Species:Year 923.84375 22 0.0000000

TAKEAWAY Data are funky with P-values all being 0 which seems strange. Furthermore the QQ plots don’t show completely normalized distributions and deviate from expected with outliers.

3.9 Herbivore count (negative binomial) ~ Percent spodoptera survial

ANOVA results of Herbivore count (negative binomial) ~ Percent spodoptera survival
Chisq Df Pr(>Chisq)
Spod_perc_survival 3.344583 1 0.0674265
Species 224.006551 13 0.0000000
Year 130.915138 2 0.0000000
Spod_perc_survival:Species 45.242953 13 0.0000191

TAKEAWAY Residulas look very good, however results are slightly different from 2.3.1 with log-transformed herbivore count in that Percent survival alone is no longer significant. Other significant effects remain the same.

3.10 Conclusions

All Herbivore count ~ Spodoptera mass give qualitatively the same significant results: Species, Year, Mass * Species, Species * Year except for Poisson transformations which doesn’t make too much sense.

log and negative binomial Herbivore count ~ Spodoptera survival give qualitatively the same significant results: Species, Year, Species x Year, Survival x Species x Year. Again, Poisson gives strange results. Since there is no 2-way Species x Survival interaction but a 3-way Species x Survival x Year interaction, this suggests that the 2-way varies among years.

4- Fixed-effect level analysis

Since there are both 2-way and 3-way effects of “Year”, “Species” “Spodoptera mass/survival” on herbivbore counts, I want to dig deeper by analyses which species and/or which years when combined with Spodoptera have an effect on Herbivores.

To do this, I will investigate the following:

4.1 Herbivore count ~ Avg Spodoptera mass (separate species)

Since there is a significant interactive effect of both Spodoptera mass x Species and Species x Year, I will keep species constant and look at solo effects of mass and year as well as interactive effects of mass x year.

Acmispon glaber

ANOVA results of ACMGLA log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 9.340308 1 0.0022417 0.89935746919726
Year 109.757190 1 0.0000000
log(Avg_spod_mass_mg + 1):Year 1.382564 1 0.2396649

TAKEAWAY

  • Significance in Average Spodoptera Mass and Year but not in interactive effect
  • R2 = .9

Artemesia californica

ANOVA results of ARTCAL log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 0.0012954 1 0.9712885 0.35447470605641
Year 4.1044616 2 0.1284480
log(Avg_spod_mass_mg + 1):Year 0.1278513 2 0.9380747

TAKEAWAY

  • No significant effect of mass or year (or interactive effects) on herbivore count
  • R2 = .3

Encelia californica

ANOVA results of ENCCAL log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 1.113671 1 0.2912854 0.921800118612757
Year 91.841180 2 0.0000000

TAKEAWAY

  • Significant effect of year on herbivore count
  • No significant effect of mass or interactive effects on herbivore count
  • Log transformation does not result in normally distributed data
  • R2 = .74

Eriogonum fasciculatum

ANOVA results of ERIFAS log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 4.956263 1 0.0259963 0.251820703961022
Year 6.219642 2 0.0446089
log(Avg_spod_mass_mg + 1):Year 2.487654 2 0.2882789

TAKEAWAY

  • Significant effects of Spodoptera mass, year and interactive effects on Herbivore count
  • R2 = .5

Grindelia camporum

ANOVA results of GRICAM log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 1.764172 1 0.1841050 0.614973257287861
Year 15.029106 2 0.0005451

TAKEAWAY

  • Interactive effect not computable- will ignore.
  • Spodoptera mass insignificant, year significant.
  • Non-normal residual distribution.
  • R2 = .57

Isocoma menziesii

ANOVA results of ISOMEN log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 0.3680424 1 0.5440733 0.895423341283106
Year 32.0362231 2 0.0000001

TAKEAWAY

  • Interactive effect not computable- will ignore.
  • Year is significant
  • R2 = .88

Malacothamnus_fasciculatus

ANOVA results of MALFAS log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 181.649 1 0 0.999582420289744
Year 11474.989 2 0

TAKEAWAY

  • Interactive effect not computable- will ignore.
  • Mass and Year are strongly significant
  • Variance not normally distributed.
  • R2 = .92

Malacothrix_saxatilis

ANOVA results of MALSAX log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 6.441778 1 0.0111467 0.348050897277051
Year 8.537983 2 0.0139959

TAKEAWAY

  • Interactive effect not computable- will ignore.
  • Mass and year significant
  • R2 = .91

Diplacus aurantiacus

ANOVA results of MIMAUR log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 0.3207125 1 0.5711798 0.0549121966468426
Year 1.0291497 2 0.5977547

TAKEAWAY

  • Interactive effect not computable- will ignore.
  • Neither mass or year are significant
  • Data are normally distributed
  • R2 = .30

Mirabilis laevis

ANOVA results of MIRLAE log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 0.0027771 1 0.9579726 0.873378670562878
Year 46.8643081 2 0.0000000

TAKEAWAY

  • Interactive effect not computable- will ignore.
  • Spodoptera mass insignicant but year is significant.
  • Data normally distributed.
  • R2 = .81

Salvia apiana

ANOVA results of SALAPI log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 0.1287276 1 0.7197549 0.803740361576943
Year 68.4867543 2 0.0000000
log(Avg_spod_mass_mg + 1):Year 4.1860685 2 0.1233124

TAKEAWAY

  • Year significant but mass and interactive effects not significant.
  • Data normally distributed.
  • R2 = .76

Sisyrinchium_bellum

ANOVA results of SISBEL log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 0.6929966 1 0.4051470 0.814830030135829
Year 32.5974782 2 0.0000001

TAKEAWAY

  • Interactive effect not computable- will ignore.
  • Year significant but mass is not
  • Data normally distributed
  • R2 = .60

Stachys ajugoides var rigida (NOT COMPUTABLE- WILL IGNORE)

ANOVA Significance Table

ANOVA results for Herbivore Count ~ Spodoptera Mass * Year for each species
Species Mass_spod Year Mass_spod.year
Acmispon_glaber ** *** ns
Artemesia_californica ns ns ns
Encelia_californica ns *** ns
Eriogonum_fasciculatum ** *** **
Grindelia_camporum ns ** n/a
Isocoma_menziesii ns *** n/a
Malacothamnus_fasciculatus ** *** n/a
Malacothrix_saxatilis *** *** n/a
Diplacus_aurantiacus ns ns n/a
Mirabilis_laevis ns *** n/a
Salvia_apiana ns *** ns
Sisyrinchium_bellum ns *** n/a
Stachys_ajugoides n/a n/a n/a

4.2 Herbivore count ~ Avg Spodoptera mass (separate years)

2021 log Herbivore count ~ log Avg Spodoptera mass

ANOVA results of log Herbivore Count ~ log Average Spodoptera mass- 2021
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 0.0321914 1 0.8576081 0.923142278229597
Species 301.4838302 12 0.0000000
log(Avg_spod_mass_mg + 1):Species 22.2665189 12 0.0346391

TAKEAWAY

  • Significant effects for species and mass*species
  • Data are normally distributed
  • R2 = .93
  • 1 species dropped from analysis

2022 log Herbivore count ~ log Avg Spodoptera mass

ANOVA results of log Herbivore Count ~ log Average Spodoptera mass- 2022
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 0.442487 1 0.5059245 0.686117369837719
Species 140.582650 13 0.0000000
log(Avg_spod_mass_mg + 1):Species 17.527923 13 0.1762801

TAKEAWAY

  • Species and Mass*species are significant
  • Data are normally distributed
  • No species dropped from analysis
  • R2 = .94

2023 Herbivore count ~ Avg Spodoptera mass

ANOVA results of log Herbivore Count ~ log Average Spodoptera mass- 2023
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 6.932155 1 0.008466 0.71733290098759
Species 62.652181 11 0.000000

TAKEAWAY

  • Interactive effect is ignored as too many species are rank deficient.
  • Species and Mass significant by themselves
  • Residuals not normally distributed
  • 3 species dropped from analysis
  • R2 = .82

ANOVA Significance Table

ANOVA results for Herbivore Count ~ Spodoptera Mass * Species for each year
Year Spod_mass Species Spod_mass.Species
2021 ns *** **
2022 ns *** *
2023 ** *** n/a

4.3 Herbivore count ~ Percent Spodoptera Survival (separate years)

Since there is a species x year interaction and a 3-way interaction but there is not a survival x species this suggests that the 2-way survival x species interaction varies among years. I will try each year separately to see if this is indeed true.

2021 Herbivore count (negative binomial) ~ Percent survival

ANOVA results of log Herbivore count ~ Percent spodoptera survival- 2021
Chisq Df Pr(>Chisq) Conditional_R2
Spod_perc_survival 5.505291 1 0.018959 0.777558405218547
Species 304.038136 13 0.000000

TAKEAWAY

  • Model convergence problem when interactive effect included. Will ingnore
  • Survival and Species significant
  • Data are close to normally distributed
  • R2 = .77

2022 log Herbivore count ~ Percent Survival

ANOVA results of log Herbivore count ~ Percent spodoptera survival- 2022
Chisq Df Pr(>Chisq) Conditional_R2
Spod_perc_survival 0.0835764 1 0.7725081 0.56007832535803
Species 177.3418229 13 0.0000000
Spod_perc_survival:Species 13.8601115 13 0.3837689

TAKEAWAY

  • Species significant
  • Data normally distributed
  • R2 = .73

2023 Herbivore count (negative biniomial) ~ Percent survival

ANOVA results of log Herbivore count ~ Percent spodoptera survival- 2023
Chisq Df Pr(>Chisq) Conditional_R2
Spod_perc_survival 1.175064 1 0.2783633 0.60444676733722
Species 252.148852 13 0.0000000

TAKEAWAY

  • Similar to 2021, can’t compute interactive effect. Will ignore in model.
  • Species significant
  • Data close to normally distributed
  • R2 = .60

ANOVA Significance Table

ANOVA results for Herbivore Count ~ Spodoptera Survival * Species for each year
Year Percent_survival Species Percent_survival.Species
2021 ** *** n/a
2022 ns *** ns
2023 ns *** n/a

4.4 Conclusions

  • The overall Species*Mass effect on Herbivore count is driven primarily ACMGLA, ERIFAS, MALFAS and MALSAX in 2021 and 2022.

  • The overall Survival*Species effect on Herbivore count is driven primarily by 2021

  • Species differ widely in the effect of spodoptera mass on Herbivore counts across across years for most analyses. Duh.

5- Regression and LMER analysis of Spodoptera mass and survival across all species

I will condense species-level data to means, plot regression lines and run LMER analysis for the following equation:

5.1 Species means table

Species means for log Herbivore count, log Spodoptera Mass and Spodoptera Percent Survival across years
Species Year Mean_log_abundance_herb Mean_log_spodmass_mg Spod_perc_survival
Acmispon_glaber 2021 4.2467113 3.8869407 0.4333333
Artemisia_californica 2021 4.3153377 1.1209790 0.2125000
Diplacus_aurantiacus 2021 2.0140065 0.8086717 0.0875000
Encelia_californica 2021 2.5502693 1.1550286 0.0375000
Eriogonum_fasciculatum 2021 1.9144399 3.6982949 0.3750000
Grindelia_camporum 2021 4.5434721 1.8720904 0.2187500
Isocoma_menziesii 2021 4.4602608 1.3165042 0.1093750
Malacothamnus_fasciculatus 2021 5.4922620 1.8363750 0.0500000
Malacothrix_saxatilis 2021 1.9965735 2.0661566 0.0892857
Mirabilis_laevis 2021 3.3784141 1.5880088 0.1000000
Salvia_apiana 2021 4.3159992 2.6836486 0.1470588
Salvia_mellifera 2021 4.1603955 2.6349589 0.0555556
Sisyrinchium_bellum 2021 1.5891719 0.9345423 0.0735294
Stachys_ajugoides_var_rigida 2021 2.1385920 NA 0.0000000
Acmispon_glaber 2022 0.8476371 2.3177974 0.6590909
Artemisia_californica 2022 3.6260061 0.5743303 0.3250000
Diplacus_aurantiacus 2022 2.4177707 0.3724995 0.3157895
Encelia_californica 2022 0.8493928 0.7121052 0.1617647
Eriogonum_fasciculatum 2022 2.6050648 1.0815205 0.5666667
Grindelia_camporum 2022 1.9228855 1.5251025 0.6500000
Isocoma_menziesii 2022 2.1668175 0.4432413 0.1250000
Malacothamnus_fasciculatus 2022 1.8823814 1.0968257 0.1730769
Malacothrix_saxatilis 2022 1.3647077 0.5802378 0.4230769
Mirabilis_laevis 2022 1.9885712 0.4643244 0.0694444
Salvia_apiana 2022 2.7630305 0.7028663 0.3437500
Salvia_mellifera 2022 2.1963871 0.8529361 0.1911765
Sisyrinchium_bellum 2022 0.5201210 0.2725614 0.2692308
Stachys_ajugoides_var_rigida 2022 0.8089944 0.5606482 0.1250000
Acmispon_glaber 2023 6.2102007 NA 0.0000000
Artemisia_californica 2023 4.2371805 0.7945135 0.0500000
Diplacus_aurantiacus 2023 2.4194210 0.0953102 0.1666667
Encelia_californica 2023 4.1244557 0.8628385 0.0441176
Eriogonum_fasciculatum 2023 2.8052491 0.4762438 0.5333333
Grindelia_camporum 2023 3.4339872 1.7917595 1.0000000
Isocoma_menziesii 2023 2.8497770 0.0953102 0.3333333
Malacothamnus_fasciculatus 2023 3.7392930 0.4874980 0.1666667
Malacothrix_saxatilis 2023 1.7969941 0.4396952 0.1500000
Mirabilis_laevis 2023 2.1819190 0.0953102 0.0147059
Salvia_apiana 2023 2.8818480 0.6472535 0.6166667
Salvia_mellifera 2023 3.1927864 0.1551344 0.2166667
Sisyrinchium_bellum 2023 1.2569468 0.0953102 0.1875000
Stachys_ajugoides_var_rigida 2023 2.4596037 NA 0.0000000

5.2 Linear regression

log Herbivore ~ log Avg Spodoptera mass * year

TAKEAWAY

A slight positive trend between increased spodoptera mass and increased herbivore abundance exists across all years of the data regardless of species. While there may not be a strong correlation between Spodoptera biomass and Herbivore abundance, there does seem to be a general pattern of decreasing Spodoptera biomass over time. One way to interpret this is that as the plants aged, they became more resistant/lower quality to Spodoptera, however that this did not seem to have a concomitant effect on Herbivore abundance. I will run the lmer to determine this quantitatively.

log Herbivore ~ Percent Spodoptera survival * year

TAKEAWAY

Slight negative correlation between Spodoptera survival and Herbivore abundance. I don’t see the same interannual patterns as with Spodoptera mass. Will run LMER to quantify.

5.3 LMER

log Herbivore ~ log Avg Spodoptera mass * year

ANOVA results of Means for log Herbivore count ~ log Spodoptera mass (All Species)
Chisq Df Pr(>Chisq) Conditional_R2
Mean_log_spodmass_mg 0.9778221 1 0.3227371 0.380431274472226
Year 13.9452254 2 0.0009372
Mean_log_spodmass_mg:Year 2.6723995 2 0.2628426

TAKEAWAY

No significant effect of log Spodoptera biomass on log Herbivore counts across all species. Furthermore, there is no significant Herbivore ~ Year or interactive Herbivore ~ Spodoptera mass * Year effect. R2 = .18

log Herbivore ~ Spodoptera percent survival * year

ANOVA results of Means for log Herbivore count ~ Spodoptera percent survival (All Species)
Chisq Df Pr(>Chisq) Conditional_R2
Spod_perc_survival 0.0036584 1 0.9517696 0.279261525771854
Year 13.9042734 2 0.0009566
Spod_perc_survival:Year 0.6133591 2 0.7358864

TAKEAWAY

No significant effect of percent survival, year or the interactive effect on log herbivore count. Residuals look good. R2 = .11.

6- Percent Survival deeper dive

I want to look deeper into the Percent Survival data and understand what may be driving the negative relationship for some species- is it a particular year? It is odd to see this.

To investigate, I am going to look at a LMER analysis for the following equation:

log Herbivore Count ~ Percent Spodoptera survival + (1|Species) + (1|Year)

The reason for looking at Species and Year as random effects is to see to what extent particular species and/or particular years may be driving the negative relationship we are seeing. Accounting for the residuals due to Species and Year should allow us to see the significance of the main effect- percent survival.

6.1 Survival regression line

I will plot a linear regression line for all survival data combined to visualize patterns irrespective of species and year.

TAKEAWAY

There is no strong effect overall. Will run LMER to quantify

Survival LMER

ANOVA results of log Herbivore count ~ Percent spodoptera survival (Year and Species = Random Effects)
Chisq Df Pr(>Chisq)
Spod_perc_survival 0.00226 1 0.9620836

TAKEAWAY

Assuming I am doing this correctly, when Species and Year are random effects, there is very little to no effect of Percent Survival on Herbivore count with a p = .8.

7- Overall Herbivore ~ Biomass relationship

Since we have low statistical power for any one species due to small numbers of spodoptera surviving and/or low population of plants due to mortality, we will do an overall analysis similar to section 5 but with raw data instead. Specifically, we will do the following:

7.1 Herbivore ~ Spodoptera * Year with species as a random effect.

I will look at herbivore ~ spodoptera first with year as a fixed effect then move to each year by itself if year has a strong influence.

7.1.1 Log herbivore ~ log spodoptera * year with species as a random effect

ANOVA results of log Herbivore count ~ log Avg Spod mass (Species = random effect)
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 1.270524 1 0.2596684 0.497145412520696
Year 83.228207 2 0.0000000
log(Avg_spod_mass_mg + 1):Year 13.167246 2 0.0013828

TAKEAWAY: Interesting results when treating “Species” as a random factor. Strong significance across years as well as an interactive efffect of spodoptera * year. Also, the R2 is very high (.45) suggesting much of the variance in herbivores is explained by Spodoptera biomass * year.

Since there is a strong effect of year on the data set, I will rerun the analyses for each year alone.

7.1.2 Log herbivore ~ log spodoptera * year with species as a random effect (2021)

ANOVA results of 2021 log Herbivore count ~ log Avg Spod mass (Species = random effect)
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 0.065054 1 0.7986791 0.889078637201039

TAKEAWAY: No strong effect of spod mass on herbivore count for 2021 data.

7.1.3 Log herbivore ~ log spodoptera * year with species as a random effect (2022)

ANOVA results of 2022 log Herbivore count ~ log Avg Spod mass (Species = random effect)
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 0.1158667 1 0.7335614 0.699889578392609

TAKEAWAY: Same as 2021. No significant effect of spodoptera mass on herbivore count.

7.1.4 Log herbivore ~ log spodoptera * year with species as a random effect (2023)

ANOVA results of 2023 log Herbivore count ~ log Avg Spod mass (Species = random effect)
Chisq Df Pr(>Chisq) Conditional_R2
log(Avg_spod_mass_mg + 1) 6.575338 1 0.0103401 0.826320055870446

TAKEAWAY: Data are significant for 2023 with spodoptera showing a strong relationship to herbivore count when counting species as a random effect.

7.2 Z-score transformations for herbivore count and spodoptera biomass

Raw data was z-transformed by species by year. I will then run lmer analyses for all years and each year to see which if year is having an effect alone and, if so, which year(s) are having the strongest effect on herbivore counts.

7.2.1 Plot z-transformed regression lines for all years and individual years from 2021-2023

TAKEAWAY: Very interesting trends present. It seems that a slight negative relationship exists in 2021 but that this relationship turns more positive for 2022 and 2023. In other words, year seems to be having an effect. Will rum lmer on z-transformed data for all years and each year individually.

7.2.2 Run lmer for all years

ANOVA results of z-Herbivore count ~ z-Spod mass across species- All Years
Chisq Df Pr(>Chisq) Conditional_R2 Marginal_R2
Z_spodmass 0.0972815 1 0.7551166 0.0920802246145703 0.0920802165327786
Year 12.6960156 2 0.0017502
Z_spodmass:Year 3.2659229 2 0.1953502

TAKEAWAY: There is no effect for Spodoptera biomass or interactive effect for Spodoptera biomass * year but there IS a year effect. Marginal and conditional R2 very low.

7.2.3 Run lmer for 2021

ANOVA results of 2021 z-Herbivore count ~ z-Spod mass across species
Chisq Df Pr(>Chisq) Conditional_R2 Marginal_R2
Z_spodmass 1.727236 1 0.1887643 0.0234273823882526 0.0234273807816507

TAKEAWAY: No strong effect for 2021 data. Data normally distributed, Conditional and marginal R2 very low.

7.2.4 Run lmer for 2022

ANOVA results of 2022 z-Herbivore count ~ z-Spod mass across species
Chisq Df Pr(>Chisq) Conditional_R2 Marginal_R2
Z_spodmass 3.922081 1 0.0476559 0.0341281805317853 0.0341281768817394

TAKEAWAY: Not significant but closer to significance than 2021. Data are not entirely normally distributed so that is something to consider. R2 still very low.

7.2.5 Run lmer for 2023 data

ANOVA results of 2023 z-Herbivore count ~ z-Spod mass across species
Chisq Df Pr(>Chisq) Conditional_R2 Marginal_R2
Z_spodmass 2.498868 1 0.1139281 0.120052665510907 0.0514919549852756

TAKEAWAY: Similar to the log-transformed data with species as random factor, 2023 is the year that seems to be driving the year effect. p = .03. Interestingly the conditional R2 = .3 while the marginal R2 is quite low at .08.

7.3 Z-score transformation of log transformation for herbivore count and spodoptera biomass

I will do the same but instead I will calcualte z-scored based on log-transformed spodoptera biomass and herbivore counts.

7.3.1 Plot z-transformed of log-transformed regression lines for all years and individual years from 2021-2023

TAKEAWAY: Qualitatively similar to z-transformed plots with positive slopes for all years, 2022 and 2023 but negative slope for 2021.

7.3.2 Run lmer for all years

ANOVA results of z-logHerbivore count ~ z-logSpod mass across species- All Years
Chisq Df Pr(>Chisq) Conditional_R2 Marginal_R2
Z_log_spodmass 1.312883 1 0.2518734 0.239771387708757 0.23977138313791
Year 53.177378 2 0.0000000
Z_log_spodmass:Year 2.791168 2 0.2476883

TAKEAWAY: Similar to non-log transformed data in terms of significance, however the residuals look much better.

7.3.3 Run lmer for 2021

ANOVA results of 2021 z-logHerbivore count ~ z-logSpod mass across species
Chisq Df Pr(>Chisq) Conditional_R2 Marginal_R2
Z_log_spodmass 0.1915476 1 0.6616317 0.00265332469552767 0.00265332398255276

TAKEAWAY: Similar to non-log transformed data in that there is no signficant relationship with extremely low R2 values. Difference is that the p-value is much higher with log-transformation (p = .76)

7.3.4 Run lmer for 2022

ANOVA results of 2022 z-logHerbivore count ~ z-logSpod mass across species
Chisq Df Pr(>Chisq) Conditional_R2 Marginal_R2
Z_log_spodmass 2.006934 1 0.1565815 0.0177593900142746 0.0177593886318746

TAKEAWAY: Qualitatively similar results to non-log transformed 2022 data. Approaching significance but still not. Should be noted that p-value for log-transformed is higher (p = .18) compared to non log-transformed (p = .08). Low R2 values.

7.3.5 Run lmer for 2023

ANOVA results of 2023 z-logHerbivore count ~ z-logSpod mass across species
Chisq Df Pr(>Chisq) Conditional_R2 Marginal_R2
Z_log_spodmass 3.002466 1 0.0831379 0.261390204569805 0.059037987076671

TAKEAWAY: Same story as non log-transforemd data. Signficant relationship for 2023 with p-values nearly identical (.28 vs .30.) Similar R2 values as without log transformation.

8- Least Square values for Spodoptera mass and Survival

We now want to turn our attention to how species and year are affecting both spodoptera mass and survival since it appears that these values are changing year to year and, possibly, also differentially with respect to species. I will run the following model to assess this:

log spodoptera mass ~ species x year spodoptera survival ~ species x year

From these models, I will extract the marginal or least square values with errors for each year. This will provide me with the average each year accounting for species variations.

8.1 Least estimates for log spodoptera mass ~ species x year

Species Year lsmean SE df lower.CL upper.CL
Acmispon_glaber 2021 3.8869396 0.1449441 217 3.6012611 4.1726182
Artemisia_californica 2021 1.1209787 0.1513892 217 0.8225972 1.4193602
Diplacus_aurantiacus 2021 0.8086725 0.2049820 217 0.4046621 1.2126830
Encelia_californica 2021 1.1550288 0.3550392 217 0.4552621 1.8547955
Eriogonum_fasciculatum 2021 3.6982924 0.1775196 217 3.3484091 4.0481757
Grindelia_camporum 2021 1.8720895 0.1897764 217 1.4980484 2.2461305
Isocoma_menziesii 2021 1.3165052 0.2245465 217 0.8739339 1.7590765
Malacothamnus_fasciculatus 2021 1.8363756 0.3550392 217 1.1366089 2.5361423
Malacothrix_saxatilis 2021 2.0661578 0.2510506 217 1.5713480 2.5609675
Mirabilis_laevis 2021 1.5880090 0.2049820 217 1.1839985 1.9920195
Salvia_apiana 2021 2.6836501 0.1897764 217 2.3096091 3.0576912
Salvia_mellifera 2021 2.6349606 0.3550392 217 1.9351939 3.3347273
Sisyrinchium_bellum 2021 0.9345400 0.2898883 217 0.3631829 1.5058972
Stachys_ajugoides_var_rigida 2021 NA NA NA NA NA
Acmispon_glaber 2022 2.3177983 0.1587783 217 2.0048531 2.6307435
Artemisia_californica 2022 0.5743300 0.1513892 217 0.2759485 0.8727115
Diplacus_aurantiacus 2022 0.3724990 0.1587783 217 0.0595538 0.6854442
Encelia_californica 2022 0.7121046 0.2049820 217 0.3080941 1.1161151
Eriogonum_fasciculatum 2022 1.0815205 0.1341922 217 0.8170335 1.3460074
Grindelia_camporum 2022 1.5251029 0.1587783 217 1.2121577 1.8380481
Isocoma_menziesii 2022 0.4432422 0.2898883 217 -0.1281149 1.0145993
Malacothamnus_fasciculatus 2022 1.0968217 0.2510506 217 0.6020120 1.5916315
Malacothrix_saxatilis 2022 0.5802362 0.1587783 217 0.2672911 0.8931814
Mirabilis_laevis 2022 0.4643246 0.2510506 217 -0.0304852 0.9591343
Salvia_apiana 2022 0.7028661 0.1449441 217 0.4171876 0.9885447
Salvia_mellifera 2022 0.8529361 0.1897764 217 0.4788951 1.2269771
Sisyrinchium_bellum 2022 0.2725614 0.1673671 217 -0.0573117 0.6024346
Stachys_ajugoides_var_rigida 2022 0.5606523 0.2898883 217 -0.0107048 1.1320094
Acmispon_glaber 2023 NA NA NA NA NA
Artemisia_californica 2023 0.7945098 0.2510506 217 0.2997001 1.2893196
Diplacus_aurantiacus 2023 0.0953085 0.1897764 217 -0.2787325 0.4693496
Encelia_californica 2023 0.8628367 0.2898883 217 0.2914796 1.4341938
Eriogonum_fasciculatum 2023 0.4762447 0.1392578 217 0.2017736 0.7507158
Grindelia_camporum 2023 1.7917571 0.5021012 217 0.8021376 2.7813767
Isocoma_menziesii 2023 0.0953088 0.1775196 217 -0.2545746 0.4451921
Malacothamnus_fasciculatus 2023 0.4874997 0.2049820 217 0.0834892 0.8915102
Malacothrix_saxatilis 2023 0.4396962 0.2898883 217 -0.1316609 1.0110533
Mirabilis_laevis 2023 0.0953080 0.5021012 217 -0.8943116 1.0849275
Salvia_apiana 2023 0.6472534 0.1341922 217 0.3827665 0.9117404
Salvia_mellifera 2023 0.1551342 0.2049820 217 -0.2488763 0.5591447
Sisyrinchium_bellum 2023 0.0953120 0.2898883 217 -0.4760451 0.6666691
Stachys_ajugoides_var_rigida 2023 NA NA NA NA NA

TAKEAWAY:

Least Square means range differs for each year in the following manner:

  • 2021: .095 - 3.88
  • 2022: .095 - 2.31
  • 2023: .095 - 1.79

Confidence intervals range as follows:

  • 2021: 3.6 - 4.1
  • 2022: .276 - 2.6
  • 2023: .292 - 1.791

Overall, the results seem to suggest that (ignoring species differences) that Spodoptera mass is affected by year, decreasing over time as the plants get older and, in theory, they become more resistant.

8.2 Least square estimates for spodoptera survival ~ species * year

Species Year lsmean SE df lower.CL upper.CL
Acmispon_glaber 2021 0.4333334 0.0618162 572 0.3119189 0.5547478
Artemisia_californica 2021 0.2125000 0.0535344 572 0.1073520 0.3176480
Diplacus_aurantiacus 2021 0.0875000 0.0535344 572 -0.0176480 0.1926479
Encelia_californica 2021 0.0374999 0.0535344 572 -0.0676481 0.1426479
Eriogonum_fasciculatum 2021 0.3750001 0.0598533 572 0.2574411 0.4925591
Grindelia_camporum 2021 0.2187500 0.0598533 572 0.1011910 0.3363090
Isocoma_menziesii 2021 0.1093749 0.0598533 572 -0.0081841 0.2269339
Malacothamnus_fasciculatus 2021 0.0499999 0.0618162 572 -0.0714145 0.1714143
Malacothrix_saxatilis 2021 0.0892857 0.0639858 572 -0.0363901 0.2149616
Mirabilis_laevis 2021 0.1000000 0.0535344 572 -0.0051480 0.2051480
Salvia_apiana 2021 0.1470588 0.0580662 572 0.0330098 0.2611078
Salvia_mellifera 2021 0.0555555 0.0564302 572 -0.0552802 0.1663912
Sisyrinchium_bellum 2021 0.0735293 0.0580662 572 -0.0405196 0.1875783
Stachys_ajugoides_var_rigida 2021 0.0000000 0.0549251 572 -0.1078795 0.1078795
Acmispon_glaber 2022 0.6590910 0.0721858 572 0.5173095 0.8008725
Artemisia_californica 2022 0.3250001 0.0535344 572 0.2198522 0.4301481
Diplacus_aurantiacus 2022 0.3157895 0.0549251 572 0.2079100 0.4236691
Encelia_californica 2022 0.1617648 0.0580662 572 0.0477158 0.2758138
Eriogonum_fasciculatum 2022 0.5666668 0.0618162 572 0.4452524 0.6880812
Grindelia_camporum 2022 0.6500001 0.0757091 572 0.5012984 0.7987017
Isocoma_menziesii 2022 0.1249999 0.0639858 572 -0.0006759 0.2506758
Malacothamnus_fasciculatus 2022 0.1730770 0.0664012 572 0.0426570 0.3034970
Malacothrix_saxatilis 2022 0.4230770 0.0664012 572 0.2926570 0.5534970
Mirabilis_laevis 2022 0.0694444 0.0564302 572 -0.0413913 0.1802801
Salvia_apiana 2022 0.3437500 0.0598533 572 0.2261910 0.4613090
Salvia_mellifera 2022 0.1911765 0.0580662 572 0.0771276 0.3052255
Sisyrinchium_bellum 2022 0.2692309 0.0664012 572 0.1388109 0.3996509
Stachys_ajugoides_var_rigida 2022 0.1250002 0.0977400 572 -0.0669729 0.3169732
Acmispon_glaber 2023 0.0000000 0.1692906 572 -0.3325071 0.3325071
Artemisia_californica 2023 0.0500000 0.0535344 572 -0.0551479 0.1551480
Diplacus_aurantiacus 2023 0.1666666 0.0564302 572 0.0558310 0.2775023
Encelia_californica 2023 0.0441175 0.0580662 572 -0.0699314 0.1581665
Eriogonum_fasciculatum 2023 0.5333334 0.0618162 572 0.4119190 0.6547478
Grindelia_camporum 2023 0.9999997 0.2394131 572 0.5297637 1.4702357
Isocoma_menziesii 2023 0.3333333 0.0691126 572 0.1975879 0.4690788
Malacothamnus_fasciculatus 2023 0.1666668 0.0691126 572 0.0309213 0.3024122
Malacothrix_saxatilis 2023 0.1500000 0.0757091 572 0.0012983 0.2987016
Mirabilis_laevis 2023 0.0147058 0.0580662 572 -0.0993432 0.1287548
Salvia_apiana 2023 0.6166666 0.0618162 572 0.4952522 0.7380810
Salvia_mellifera 2023 0.2166666 0.0618162 572 0.0952522 0.3380810
Sisyrinchium_bellum 2023 0.1875001 0.0846453 572 0.0212466 0.3537537
Stachys_ajugoides_var_rigida 2023 0.0000000 0.0846453 572 -0.1662535 0.1662535

TAKEAWAY:

Least square mean ranges according to year:

  • 2021: 0 - .43
  • 2022: 0 - .99
  • 2023: 0 - .61

Confidence intervals according to year:

  • 2021: -.1 - .55
  • 2022: -.16 - .79
  • 2023: -.33 - 1.47

Overall, the results seem to suggest that (ignoring species differences) that survival does not change over time as the lsmean ranges remain fairly consistent. It should be noted, however, that the confidence intervals grow each year.

Nell, Colleen S., and Kailen A. Mooney. 2019. “Plant Structural Complexity Mediates Trade-Off in Direct and Indirect Plant Defense by Birds.” Ecology 100 (10). https://doi.org/10.1002/ecy.2853.
Pratt, Jessica D., Andrew Datu, Thi Tran, Daniel C. Sheng, and Kailen A. Mooney. 2017. “Genetically Based Latitudinal Clines in Artemisia Californica Drive Parallel Clines in Arthropod Communities.” Ecology 98 (1): 79–91. https://doi.org/10.1002/ecy.1620.
Rogers, Lee E., W. T. Hinds, and Ray L. Buschbom. 1976. “A General Weight Vs. Length Relationship for Insects1.” Annals of the Entomological Society of America 69 (2): 387–89. https://doi.org/10.1093/aesa/69.2.387.
You, Yanrong, Chunpeng An, and Chuanyou Li. 2020. “Insect Feeding Assays with Spodoptera Exigua on Arabidopsis Thaliana.” BIO-PROTOCOL 10 (5). https://doi.org/10.21769/BioProtoc.3538.