Experimental resource supplementation shifts ant-mediated defense on silver cholla

library(dplyr)


Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(ggplot2)

Warning: package 'ggplot2' was built under R version 4.3.3

library(performance)
library(glmmTMB)
library(emmeans)
library(tidyr)
library(vegan)

Loading required package: permute

Loading required package: lattice

Warning: package 'lattice' was built under R version 4.3.3

This is vegan 2.6-4

library(gridExtra)


Attaching package: 'gridExtra'

The following object is masked from 'package:dplyr':

    combine

Data Cleaning

#load in the visual survey data
arth <- read.csv("Raw Data/visual_surveys.csv")
arth$plant.id <- as.character(arth$plant.id)

#load plant measurements
plt <- read.csv("Raw Data/plant_data.csv")

#filter out the silver cholla
sc <- filter(plt, Species == "silvercholla") 

#18trap is a plant with a pitfall trap that was never surveys
sc <- sc %>%
  filter(PlantID != "18trap")

#join the survey data to the plant measurements
data <- left_join(sc, arth, by = c("PlantID" = "plant.id", "year")) 
data <- select(data, 1:13, 31:103)

#remove the half survey for simplicity
data <- data %>% filter(days.since.treatment != 13)

str(data)

'data.frame':   1180 obs. of  86 variables:
 $ PlantID                   : chr  "1" "1" "1" "1" ...
 $ Species                   : chr  "silvercholla" "silvercholla" "silvercholla" "silvercholla" ...
 $ date.x                    : num  4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 ...
 $ year                      : int  2023 2023 2023 2023 2023 2023 2023 2023 2023 2023 ...
 $ x                         : int  58 58 58 58 58 58 58 58 58 39 ...
 $ y                         : int  36 36 36 36 36 36 36 36 36 28 ...
 $ z                         : int  65 65 65 65 65 65 65 65 65 54 ...
 $ treatment                 : chr  "C" "C" "C" "C" ...
 $ growth.nodes              : int  110 110 110 110 110 110 110 110 110 97 ...
 $ new.segments              : int  41 41 41 41 41 41 41 41 41 27 ...
 $ buds                      : int  49 49 49 49 49 49 49 49 49 55 ...
 $ aborted.buds              : int  NA NA NA NA NA NA NA NA NA NA ...
 $ open.flowers              : int  NA NA NA NA NA NA NA NA NA NA ...
 $ date.y                    : chr  "4.9.2023" "4.10.2023" "4.12.2023" "4.15.2023" ...
 $ time.of.day               : chr  "morning" "morning" "morning" "afternoon" ...
 $ days.since.treatment      : int  0 1 3 6 9 12 17 22 24 0 ...
 $ CB.eggs                   : int  NA NA NA NA NA NA NA NA 3 NA ...
 $ adult.Chelinidea.vittiger : int  3 2 NA 2 2 NA 1 NA NA NA ...
 $ nymph.Chelinidea.vittiger : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Notogramma.purpuratum     : int  NA NA NA NA NA NA NA NA NA 1 ...
 $ Narnia                    : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Crematogaster.depilis     : int  1 NA NA 5 2 1 1 15 10 3 ...
 $ Crematogaster.depilis.tube: int  NA NA NA NA NA NA NA NA NA NA ...
 $ Endiodioctes              : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Endiodioctes.tube         : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Eremocystus               : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Monomorium.ergatogyna     : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Camponotus.sayi           : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Solenopsis.xyloni         : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Forelius.pruinosis        : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Myrmecosystus.yellow      : int  NA NA NA NA NA NA NA NA NA NA ...
 $ geocoris                  : int  NA NA NA NA NA NA NA NA NA NA ...
 $ orange.geocorid           : int  NA NA NA 1 NA NA NA NA NA NA ...
 $ arrhysus                  : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Largus.californicus       : int  NA NA NA NA NA NA NA NA NA NA ...
 $ imm.cal                   : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Leptoglossus.clypealus    : int  NA NA NA NA NA NA NA NA NA NA ...
 $ black.parasitoid.wasp     : int  NA NA NA NA NA NA NA NA NA NA ...
 $ green.stink.bug           : int  NA NA NA NA NA NA NA NA NA NA ...
 $ stink.nymph               : int  NA NA NA NA NA NA NA NA NA NA ...
 $ small.hemip               : int  NA NA NA NA NA NA NA NA NA NA ...
 $ wingless.wasp             : int  NA NA NA NA NA NA NA NA NA NA ...
 $ long.parasitoid.wasp      : int  NA NA NA NA NA NA NA NA NA NA ...
 $ aphid.outbreak            : int  NA NA NA NA NA NA NA NA NA NA ...
 $ notes                     : chr  "" "" "" "" ...
 $ tiny.wasp                 : int  NA NA NA NA NA NA NA NA NA NA ...
 $ salticidae                : int  NA NA NA NA NA NA NA NA NA NA ...
 $ aphid                     : int  NA NA NA NA NA NA NA NA NA NA ...
 $ braconid                  : int  NA NA NA NA NA NA NA NA NA NA ...
 $ lacewing                  : int  NA NA NA NA NA NA NA NA NA NA ...
 $ shiny.black.fly           : int  NA NA NA 1 NA NA NA NA NA NA ...
 $ trupanea.                 : int  NA NA NA NA NA NA NA NA NA NA ...
 $ red.and.black.hemipteran  : int  NA NA NA NA NA NA NA NA NA NA ...
 $ other.hemip               : int  NA NA NA NA NA NA NA NA NA NA ...
 $ eremochrysa               : int  NA NA NA NA NA NA NA NA NA NA ...
 $ ladybug                   : int  NA NA NA NA NA NA NA NA NA NA ...
 $ pompilidae                : int  NA NA NA NA NA NA NA NA NA NA ...
 $ syrphid                   : int  NA NA NA NA NA NA NA NA NA NA ...
 $ chrysomelidae             : int  NA NA NA NA NA NA NA NA NA NA ...
 $ orb.spider                : int  NA NA NA NA NA NA NA NA NA NA ...
 $ antracine.bombyliidae     : int  NA NA NA NA NA NA NA NA NA NA ...
 $ caterpillar               : int  NA NA NA NA NA NA NA NA NA NA ...
 $ lynx.spider               : int  NA NA NA NA NA NA NA NA NA NA ...
 $ green.crab.spider         : int  NA NA NA NA NA NA NA NA NA NA ...
 $ ichneumoidae              : int  NA NA NA NA NA NA NA NA NA NA ...
 $ fungus.gnat               : int  NA NA NA NA NA NA NA NA NA NA ...
 $ anthomyiidae              : int  NA NA NA NA NA NA NA NA NA NA ...
 $ sptted.wing.fly           : int  NA NA NA NA NA NA NA NA NA NA ...
 $ black.geocoris            : int  NA NA NA NA NA NA NA NA NA NA ...
 $ lasioglossum              : int  NA NA NA NA NA NA NA NA NA NA ...
 $ sml.black.coccinelid      : int  NA NA NA NA NA NA NA NA NA NA ...
 $ X                         : logi  NA NA NA NA NA NA ...
 $ cranefly                  : int  NA NA NA NA NA NA NA NA NA NA ...
 $ other.spider              : int  NA NA NA NA NA NA NA NA NA NA ...
 $ tetracanthus.assasin.bug  : chr  "" "" "" "" ...
 $ redblackmoth              : int  NA NA NA NA NA NA NA NA NA NA ...
 $ redblackimmaturehemip     : chr  "" "" "" "" ...
 $ black.specid              : int  NA NA NA NA NA NA NA NA NA NA ...
 $ sarcophagidae             : int  NA NA NA NA NA NA NA NA NA NA ...
 $ green.lygus               : int  NA NA NA NA NA NA NA NA NA NA ...
 $ chrysididae               : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Berytidae                 : int  NA NA NA NA NA NA NA NA NA NA ...
 $ philodromidae             : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Trichodes                 : int  NA NA NA NA NA NA NA NA NA NA ...
 $ brown.hemip               : chr  "" "" "" "" ...
 $ assasin.bug..apiomeris.   : chr  "" "" "" "" ...

data$brown.hemip <- as.integer(data$brown.hemip)

Warning: NAs introduced by coercion

data$assasin.bug..apiomeris. <- as.integer(data$assasin.bug..apiomeris.)

Warning: NAs introduced by coercion

data$redblackimmaturehemip <- as.integer(data$redblackimmaturehemip)

Warning: NAs introduced by coercion

data$tetracanthus.assasin.bug <- as.integer(data$tetracanthus.assasin.bug)

Warning: NAs introduced by coercion

#blanks are zeroes
data[is.na(data)] <- 0

#add together tubes and plant obs here
data <- mutate(data, Crematogaster.depilis = Crematogaster.depilis + Crematogaster.depilis.tube, Endiodioctes = Endiodioctes + Endiodioctes.tube) 
data <- select(data, -Crematogaster.depilis.tube, -Endiodioctes.tube)
data <- mutate(data, ants = Crematogaster.depilis  + Endiodioctes + 
                 Monomorium.ergatogyna + Camponotus.sayi + Solenopsis.xyloni + Forelius.pruinosis + Myrmecosystus.yellow + Eremocystus, herbs.all = adult.Chelinidea.vittiger +Notogramma.purpuratum + Narnia + arrhysus + green.stink.bug + Largus.californicus +  Leptoglossus.clypealus + chrysomelidae + caterpillar + green.lygus + Berytidae + Trichodes, preds = geocoris + orange.geocorid + black.parasitoid.wasp + wingless.wasp + long.parasitoid.wasp + tiny.wasp + salticidae + braconid + lacewing + eremochrysa + ladybug + pompilidae + orb.spider + lynx.spider + green.crab.spider + ichneumoidae + sml.black.coccinelid + other.spider + tetracanthus.assasin.bug + black.specid + philodromidae + assasin.bug..apiomeris.)



#add columns for presence/absence
data$pres <- ifelse(data$ants > 0, 1, 0)
data$solpres <- ifelse(data$Solenopsis.xyloni > 0, 1, 0)
data$crempres <- ifelse(data$Crematogaster.depilis > 0, 1, 0)

str(data)

'data.frame':   1180 obs. of  90 variables:
 $ PlantID                  : chr  "1" "1" "1" "1" ...
 $ Species                  : chr  "silvercholla" "silvercholla" "silvercholla" "silvercholla" ...
 $ date.x                   : num  4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 ...
 $ year                     : int  2023 2023 2023 2023 2023 2023 2023 2023 2023 2023 ...
 $ x                        : int  58 58 58 58 58 58 58 58 58 39 ...
 $ y                        : int  36 36 36 36 36 36 36 36 36 28 ...
 $ z                        : int  65 65 65 65 65 65 65 65 65 54 ...
 $ treatment                : chr  "C" "C" "C" "C" ...
 $ growth.nodes             : num  110 110 110 110 110 110 110 110 110 97 ...
 $ new.segments             : num  41 41 41 41 41 41 41 41 41 27 ...
 $ buds                     : int  49 49 49 49 49 49 49 49 49 55 ...
 $ aborted.buds             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ open.flowers             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ date.y                   : chr  "4.9.2023" "4.10.2023" "4.12.2023" "4.15.2023" ...
 $ time.of.day              : chr  "morning" "morning" "morning" "afternoon" ...
 $ days.since.treatment     : int  0 1 3 6 9 12 17 22 24 0 ...
 $ CB.eggs                  : num  0 0 0 0 0 0 0 0 3 0 ...
 $ adult.Chelinidea.vittiger: num  3 2 0 2 2 0 1 0 0 0 ...
 $ nymph.Chelinidea.vittiger: num  0 0 0 0 0 0 0 0 0 0 ...
 $ Notogramma.purpuratum    : num  0 0 0 0 0 0 0 0 0 1 ...
 $ Narnia                   : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Crematogaster.depilis    : num  1 0 0 5 2 1 1 15 10 3 ...
 $ Endiodioctes             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Eremocystus              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Monomorium.ergatogyna    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Camponotus.sayi          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Solenopsis.xyloni        : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Forelius.pruinosis       : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Myrmecosystus.yellow     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ geocoris                 : num  0 0 0 0 0 0 0 0 0 0 ...
 $ orange.geocorid          : num  0 0 0 1 0 0 0 0 0 0 ...
 $ arrhysus                 : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Largus.californicus      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ imm.cal                  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Leptoglossus.clypealus   : num  0 0 0 0 0 0 0 0 0 0 ...
 $ black.parasitoid.wasp    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ green.stink.bug          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stink.nymph              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ small.hemip              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ wingless.wasp            : num  0 0 0 0 0 0 0 0 0 0 ...
 $ long.parasitoid.wasp     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ aphid.outbreak           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ notes                    : chr  "" "" "" "" ...
 $ tiny.wasp                : num  0 0 0 0 0 0 0 0 0 0 ...
 $ salticidae               : num  0 0 0 0 0 0 0 0 0 0 ...
 $ aphid                    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ braconid                 : num  0 0 0 0 0 0 0 0 0 0 ...
 $ lacewing                 : num  0 0 0 0 0 0 0 0 0 0 ...
 $ shiny.black.fly          : num  0 0 0 1 0 0 0 0 0 0 ...
 $ trupanea.                : num  0 0 0 0 0 0 0 0 0 0 ...
 $ red.and.black.hemipteran : num  0 0 0 0 0 0 0 0 0 0 ...
 $ other.hemip              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ eremochrysa              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ladybug                  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ pompilidae               : num  0 0 0 0 0 0 0 0 0 0 ...
 $ syrphid                  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ chrysomelidae            : num  0 0 0 0 0 0 0 0 0 0 ...
 $ orb.spider               : num  0 0 0 0 0 0 0 0 0 0 ...
 $ antracine.bombyliidae    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ caterpillar              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ lynx.spider              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ green.crab.spider        : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ichneumoidae             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ fungus.gnat              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ anthomyiidae             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ sptted.wing.fly          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ black.geocoris           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ lasioglossum             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ sml.black.coccinelid     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ X                        : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cranefly                 : num  0 0 0 0 0 0 0 0 0 0 ...
 $ other.spider             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ tetracanthus.assasin.bug : num  0 0 0 0 0 0 0 0 0 0 ...
 $ redblackmoth             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ redblackimmaturehemip    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ black.specid             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ sarcophagidae            : num  0 0 0 0 0 0 0 0 0 0 ...
 $ green.lygus              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ chrysididae              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Berytidae                : num  0 0 0 0 0 0 0 0 0 0 ...
 $ philodromidae            : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Trichodes                : num  0 0 0 0 0 0 0 0 0 0 ...
 $ brown.hemip              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ assasin.bug..apiomeris.  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ants                     : num  1 0 0 5 2 1 1 15 10 3 ...
 $ herbs.all                : num  3 2 0 2 2 0 1 0 0 1 ...
 $ preds                    : num  0 0 0 1 0 0 0 0 0 0 ...
 $ pres                     : num  1 0 0 1 1 1 1 1 1 1 ...
 $ solpres                  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ crempres                 : num  1 0 0 1 1 1 1 1 1 1 ...

data$year <- as.factor(data$year)

unique(data$year)

[1] 2023 2024
Levels: 2023 2024

data$uniID <- paste(data$PlantID, data$year)

#create datasets from survey data after one week of treatment
alldata <- data
data <- filter(data, days.since.treatment >= 7)
data <- mutate(data, efndensity = new.segments/x)

#create a dataset aggregrated to the plant level (vs surveys)
visits <- data %>% group_by(PlantID, x, year, new.segments, buds, treatment,) %>% summarise(ants = sum(ants), bugs = sum(adult.Chelinidea.vittiger), herbs = sum(herbs.all), s.xyloni = sum(Solenopsis.xyloni), attend = sum(pres), crem = sum(Crematogaster.depilis), cremattend = sum(crempres), preds = sum(preds))

`summarise()` has grouped output by 'PlantID', 'x', 'year', 'new.segments',
'buds'. You can override using the `.groups` argument.

#create a wide dataset aggregrated to the plant level (vs surveys)



data$pres <- as.factor(data$pres)



#create a long format dataset with species identity
antslong <- data %>% pivot_longer(22:29, names_to = "ant.species", values_to = "abun")
antslong <- select(antslong, 1:16, 83, 85:86)


#add EFN neighbours
sc <- mutate(sc, efn.neighbor = CE + CA + acacia, shrub.neighbor = LT + AS + CR + EC + YS + ecoop + TM + K + BB + ambsal + HH + unknown + tetradynia)

neigh <- cbind(sc$PlantID, sc$efn.neighbor)

neigh <- as.data.frame(neigh)
colnames(neigh) <- c("uniID", "neigh")

data <- left_join(data, neigh, by = "uniID")
data$neigh <- as.numeric(data$neigh)

Simple EDA

ggplot(data, aes(year, ants, color = treatment)) + geom_boxplot()

ggplot(data, aes(treatment, ants)) + geom_violin()

ggplot(data, aes(treatment, preds)) + geom_boxplot()

ggplot(data, aes(treatment, herbs.all)) + geom_boxplot()

ggplot(data, aes(treatment, adult.Chelinidea.vittiger)) + geom_boxplot()

ggplot(data, aes(treatment, adult.Chelinidea.vittiger)) + geom_violin()

ggplot(data, aes(treatment, CB.eggs)) + geom_boxplot()

ggplot(data, aes(ants, adult.Chelinidea.vittiger)) + geom_smooth(method = "lm")

`geom_smooth()` using formula = 'y ~ x'

ggplot(data, aes(ants, CB.eggs)) + geom_smooth(method = "lm")

`geom_smooth()` using formula = 'y ~ x'

Summary Calculations

sum(data$adult.Chelinidea.vittiger)

[1] 643

sum(data$ants)

[1] 1337

sum(data$herbs.all)

[1] 843

dherb <- mutate(data, herbpres = ifelse(herbs.all > 0, 1, 0))

filter(data, ants >= 1) %>% count()

    n
1 241

sum(as.numeric(as.character(data$pres)))

[1] 241

sum(dherb$herbpres)

[1] 432

sum(visits$crempres)

Warning: Unknown or uninitialised column: `crempres`.

[1] 0

sum(visits$solpres)

Warning: Unknown or uninitialised column: `solpres`.

[1] 0

cor.test(sc$new.segments, sc$growth.nodes)


    Pearson's product-moment correlation

data:  sc$new.segments and sc$growth.nodes
t = 15.899, df = 57, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8420304 0.9415951
sample estimates:
      cor 
0.9033288

sums <- data %>% select(17, 18, 20:84)
sums <- select(sums, -notes)

sums <- colSums(sums[,])
sums <- as.data.frame(sums)

#write.csv(sums, "Clean Data/counts.csv")

Ant Defense

Herbivore Abundance

These models use survey data from 1 week post treatment

#all herbivores together

#ant presence
m1 <- glmmTMB(herbs.all ~ pres + treatment + days.since.treatment +  new.segments + preds+ (1|year) + (1|uniID), family = "poisson", data = data)
check_overdispersion(m1)

# Overdispersion test

       dispersion ratio =   1.016
  Pearson's Chi-Squared = 865.728
                p-value =   0.364

No overdispersion detected.

m1int <- glmmTMB(herbs.all ~ pres + treatment * days.since.treatment + new.segments + preds +   (1|year) + (1|uniID), family = "poisson", data = data)
check_overdispersion(m1int)

# Overdispersion test

       dispersion ratio =   1.015
  Pearson's Chi-Squared = 863.368
                p-value =   0.377

No overdispersion detected.

AIC(m1, m1int)

      df      AIC
m1     8 2276.935
m1int  9 2274.430

#keep interaction
summary(m1int)

 Family: poisson  ( log )
Formula:          
herbs.all ~ pres + treatment * days.since.treatment + new.segments +  
    preds + (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  2274.4   2317.2  -1128.2   2256.4      851 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 6.437e-10 2.537e-05
 uniID  (Intercept) 4.290e-01 6.550e-01
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -1.095186   0.225124  -4.865 1.15e-06 ***
pres1                           -0.249518   0.097203  -2.567   0.0103 *  
treatmentN                       0.333740   0.245341   1.360   0.1737    
days.since.treatment             0.016137   0.009026   1.788   0.0738 .  
new.segments                     0.017241   0.003287   5.246 1.55e-07 ***
preds                           -0.115288   0.073161  -1.576   0.1151    
treatmentN:days.since.treatment -0.027597   0.013008  -2.122   0.0339 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m1int, type = 3)

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: herbs.all
                                 Chisq Df Pr(>Chisq)    
(Intercept)                    23.6664  1  1.146e-06 ***
pres                            6.5894  1    0.01026 *  
treatment                       1.8504  1    0.17373    
days.since.treatment            3.1965  1    0.07380 .  
new.segments                   27.5199  1  1.555e-07 ***
preds                           2.4832  1    0.11507    
treatment:days.since.treatment  4.5009  1    0.03388 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# pres * treatment interaction is non-significant
# days since treatment * treatment IS significant, adding year interaction doesn't impact it

ggplot(data, aes(days.since.treatment, herbs.all, fill = treatment)) + geom_smooth(method = "lm") + facet_grid(~year)

`geom_smooth()` using formula = 'y ~ x'

test(emtrends(m1int, pairwise~treatment, var = "days.since.treatment", type = "response"))

$emtrends
 treatment days.since.treatment.trend      SE  df z.ratio p.value
 C                             0.0161 0.00903 Inf   1.788  0.0738
 N                            -0.0115 0.00948 Inf  -1.209  0.2268

Results are averaged over the levels of: pres 

$contrasts
 contrast estimate    SE  df z.ratio p.value
 C - N      0.0276 0.013 Inf   2.122  0.0339

Results are averaged over the levels of: pres

emmip(m1int, days.since.treatment ~ treatment, data = data)

emmip(m1int, treatment ~ days.since.treatment, cov.reduce = range, xlab = "Days Since Treatment")

# using ant abundance instead of presence

m1abunint <- glmmTMB(herbs.all ~ ants + treatment * days.since.treatment+  new.segments + preds + (1|year) + (1|uniID), family = "poisson", data = data)
check_overdispersion(m1abunint)

# Overdispersion test

       dispersion ratio =   1.021
  Pearson's Chi-Squared = 869.165
                p-value =   0.325

No overdispersion detected.

check_collinearity(m1abunint)

# Check for Multicollinearity

Low Correlation

                           Term  VIF    VIF 95% CI adj. VIF Tolerance
                           ants 1.07 [1.02,  1.21]     1.03      0.94
                      treatment 3.25 [2.92,  3.65]     1.80      0.31
           days.since.treatment 1.97 [1.80,  2.19]     1.41      0.51
                   new.segments 1.01 [1.00, 52.46]     1.00      0.99
                          preds 1.02 [1.00,  2.21]     1.01      0.98
 treatment:days.since.treatment 4.07 [3.63,  4.58]     2.02      0.25
 Tolerance 95% CI
     [0.83, 0.98]
     [0.27, 0.34]
     [0.46, 0.56]
     [0.02, 1.00]
     [0.45, 1.00]
     [0.22, 0.28]

summary(m1abunint)

 Family: poisson  ( log )
Formula:          
herbs.all ~ ants + treatment * days.since.treatment + new.segments +  
    preds + (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  2277.0   2319.8  -1129.5   2259.0      851 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 9.996e-10 3.162e-05
 uniID  (Intercept) 4.135e-01 6.430e-01
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -1.149965   0.223254  -5.151 2.59e-07 ***
ants                            -0.021964   0.011296  -1.944   0.0518 .  
treatmentN                       0.367759   0.243361   1.511   0.1307    
days.since.treatment             0.017388   0.009044   1.923   0.0545 .  
new.segments                     0.017261   0.003246   5.318 1.05e-07 ***
preds                           -0.117645   0.073014  -1.611   0.1071    
treatmentN:days.since.treatment -0.029342   0.012952  -2.265   0.0235 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m1abunint)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: herbs.all
                                 Chisq Df Pr(>Chisq)    
ants                            3.7809  1    0.05184 .  
treatment                       0.4483  1    0.50313    
days.since.treatment            0.2766  1    0.59893    
new.segments                   28.2812  1  1.049e-07 ***
preds                           2.5962  1    0.10712    
treatment:days.since.treatment  5.1320  1    0.02349 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Crematogaster depilis presence only

m2 <- glmmTMB(herbs.all ~ crempres + treatment * days.since.treatment + new.segments + preds+  (1|year) + (1|uniID), family = "poisson", data = data)
check_overdispersion(m2)

# Overdispersion test

       dispersion ratio =   1.014
  Pearson's Chi-Squared = 862.937
                p-value =   0.381

No overdispersion detected.

summary(m2)

 Family: poisson  ( log )
Formula:          
herbs.all ~ crempres + treatment * days.since.treatment + new.segments +  
    preds + (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  2272.6   2315.4  -1127.3   2254.6      851 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 6.779e-10 2.604e-05
 uniID  (Intercept) 4.360e-01 6.603e-01
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -1.087010   0.225619  -4.818 1.45e-06 ***
crempres                        -0.330387   0.114791  -2.878   0.0040 ** 
treatmentN                       0.369922   0.244095   1.515   0.1296    
days.since.treatment             0.015958   0.009018   1.770   0.0768 .  
new.segments                     0.017304   0.003307   5.232 1.67e-07 ***
preds                           -0.116386   0.073253  -1.589   0.1121    
treatmentN:days.since.treatment -0.031679   0.012865  -2.462   0.0138 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m2)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: herbs.all
                                 Chisq Df Pr(>Chisq)    
crempres                        8.2838  1     0.0040 ** 
treatment                       0.8313  1     0.3619    
days.since.treatment            0.0063  1     0.9365    
new.segments                   27.3760  1  1.675e-07 ***
preds                           2.5244  1     0.1121    
treatment:days.since.treatment  6.0633  1     0.0138 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Crematogaster presence has same impact


m3 <- glmmTMB(herbs.all ~ solpres + treatment * days.since.treatment + new.segments + preds + (1|year) + (1|uniID), family = "poisson", data = data)
check_overdispersion(m3)

# Overdispersion test

       dispersion ratio =   1.028
  Pearson's Chi-Squared = 874.856
                p-value =   0.278

No overdispersion detected.

summary(m3)

 Family: poisson  ( log )
Formula:          
herbs.all ~ solpres + treatment * days.since.treatment + new.segments +  
    preds + (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  2280.8   2323.6  -1131.4   2262.8      851 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 8.125e-10 0.0000285
 uniID  (Intercept) 4.195e-01 0.6476635
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -1.140123   0.224287  -5.083 3.71e-07 ***
solpres                         -0.120635   0.191906  -0.629   0.5296    
treatmentN                       0.388740   0.243062   1.599   0.1097    
days.since.treatment             0.016096   0.009028   1.783   0.0746 .  
new.segments                     0.016855   0.003258   5.173 2.30e-07 ***
preds                           -0.113400   0.073261  -1.548   0.1216    
treatmentN:days.since.treatment -0.030968   0.012930  -2.395   0.0166 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m3)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: herbs.all
                                 Chisq Df Pr(>Chisq)    
solpres                         0.3952  1    0.52960    
treatment                       0.4703  1    0.49286    
days.since.treatment            0.0546  1    0.81522    
new.segments                   26.7617  1  2.302e-07 ***
preds                           2.3960  1    0.12165    
treatment:days.since.treatment  5.7364  1    0.01662 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Solenopsis does not


ggplot(data, aes(as.factor(crempres), herbs.all)) + geom_boxplot()

ggplot(data, aes(ants, herbs.all)) + geom_smooth(method = "lm") + geom_point(position = "jitter") + ylim(0, 7)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 2 rows containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 234 rows containing missing values or values outside the scale range
(`geom_point()`).

Cactus bug Abundance

#just cactus bugs
#with ant presence

m1 <- glmmTMB(adult.Chelinidea.vittiger ~ pres + treatment + days.since.treatment + new.segments + preds + (1|year) + (1|uniID), family = "poisson", data = data)

check_overdispersion(m1)

# Overdispersion test

       dispersion ratio =   0.953
  Pearson's Chi-Squared = 811.733
                p-value =   0.835

No overdispersion detected.

summary(m1)

 Family: poisson  ( log )
Formula:          
adult.Chelinidea.vittiger ~ pres + treatment + days.since.treatment +  
    new.segments + preds + (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  1979.7   2017.8   -981.9   1963.7      852 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 6.600e-10 2.569e-05
 uniID  (Intercept) 5.417e-01 7.360e-01
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -1.246078   0.228360  -5.457 4.85e-08 ***
pres1                -0.308745   0.110111  -2.804  0.00505 ** 
treatmentN           -0.142179   0.155596  -0.914  0.36084    
days.since.treatment  0.006336   0.007501   0.845  0.39824    
new.segments          0.017671   0.003733   4.734 2.20e-06 ***
preds                -0.038893   0.080597  -0.483  0.62940    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m1int_pres <- glmmTMB(adult.Chelinidea.vittiger ~ treatment * days.since.treatment + pres + preds + new.segments + (1|year) + (1|uniID), family = "poisson", data = data)
check_overdispersion(m1int_pres)

# Overdispersion test

       dispersion ratio =   0.948
  Pearson's Chi-Squared = 807.071
                p-value =   0.857

No overdispersion detected.

AIC(m1, m1int_pres)

           df      AIC
m1          8 1979.728
m1int_pres  9 1976.271

#better with interaction 
summary(m1int_pres)

 Family: poisson  ( log )
Formula:          
adult.Chelinidea.vittiger ~ treatment * days.since.treatment +  
    pres + preds + new.segments + (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  1976.3   2019.1   -979.1   1958.3      851 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 8.859e-10 2.976e-05
 uniID  (Intercept) 5.338e-01 7.306e-01
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -1.507595   0.256249  -5.883 4.02e-09 ***
treatmentN                       0.406000   0.281012   1.445   0.1485    
days.since.treatment             0.022282   0.010188   2.187   0.0287 *  
pres1                           -0.275737   0.110968  -2.485   0.0130 *  
preds                           -0.038750   0.080781  -0.480   0.6314    
new.segments                     0.017616   0.003713   4.745 2.08e-06 ***
treatmentN:days.since.treatment -0.034774   0.014897  -2.334   0.0196 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m1int_pres)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: adult.Chelinidea.vittiger
                                 Chisq Df Pr(>Chisq)    
treatment                       0.8287  1    0.36265    
days.since.treatment            0.6909  1    0.40585    
pres                            6.1744  1    0.01296 *  
preds                           0.2301  1    0.63145    
new.segments                   22.5151  1  2.085e-06 ***
treatment:days.since.treatment  5.4488  1    0.01958 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

test(emtrends(m1int_pres, pairwise ~ treatment, var = "days.since.treatment"))

$emtrends
 treatment days.since.treatment.trend     SE  df z.ratio p.value
 C                             0.0223 0.0102 Inf   2.187  0.0287
 N                            -0.0125 0.0110 Inf  -1.135  0.2565

Results are averaged over the levels of: pres 

$contrasts
 contrast estimate     SE  df z.ratio p.value
 C - N      0.0348 0.0149 Inf   2.334  0.0196

Results are averaged over the levels of: pres

m1int_crempres <- glmmTMB(adult.Chelinidea.vittiger ~ treatment * days.since.treatment + crempres + preds + new.segments + (1|year) + (1|uniID), family = "poisson", data = data)
summary(m1int_crempres)

 Family: poisson  ( log )
Formula:          
adult.Chelinidea.vittiger ~ treatment * days.since.treatment +  
    crempres + preds + new.segments + (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  1975.5   2018.3   -978.7   1957.5      851 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 6.298e-10 0.0000251
 uniID  (Intercept) 5.455e-01 0.7385850
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -1.498798   0.257046  -5.831 5.51e-09 ***
treatmentN                       0.445762   0.279675   1.594  0.11097    
days.since.treatment             0.021867   0.010175   2.149  0.03163 *  
crempres                        -0.339852   0.129452  -2.625  0.00866 ** 
preds                           -0.039883   0.080991  -0.492  0.62241    
new.segments                     0.017636   0.003743   4.711 2.46e-06 ***
treatmentN:days.since.treatment -0.039194   0.014736  -2.660  0.00782 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m1int_crempres)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: adult.Chelinidea.vittiger
                                 Chisq Df Pr(>Chisq)    
treatment                       1.1883  1   0.275668    
days.since.treatment            0.2018  1   0.653267    
crempres                        6.8923  1   0.008657 ** 
preds                           0.2425  1   0.622411    
new.segments                   22.1970  1  2.461e-06 ***
treatment:days.since.treatment  7.0740  1   0.007821 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m1int_solpres <- glmmTMB(adult.Chelinidea.vittiger ~ treatment * days.since.treatment + solpres + preds + new.segments + (1|year) + (1|uniID), family = "poisson", data = data)
summary(m1int_solpres)

 Family: poisson  ( log )
Formula:          
adult.Chelinidea.vittiger ~ treatment * days.since.treatment +  
    solpres + preds + new.segments + (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  1982.5   2025.3   -982.3   1964.5      851 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 7.298e-10 2.702e-05
 uniID  (Intercept) 5.137e-01 7.167e-01
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -1.542919   0.254479  -6.063 1.34e-09 ***
treatmentN                       0.468987   0.277875   1.688  0.09146 .  
days.since.treatment             0.021711   0.010176   2.133  0.03289 *  
solpres                         -0.063042   0.212254  -0.297  0.76646    
preds                           -0.037448   0.080912  -0.463  0.64349    
new.segments                     0.017070   0.003659   4.665 3.09e-06 ***
treatmentN:days.since.treatment -0.038818   0.014820  -2.619  0.00881 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m1int_solpres)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: adult.Chelinidea.vittiger
                                 Chisq Df Pr(>Chisq)    
treatment                       0.8089  1    0.36846    
days.since.treatment            0.2972  1    0.58563    
solpres                         0.0882  1    0.76646    
preds                           0.2142  1    0.64349    
new.segments                   21.7623  1  3.086e-06 ***
treatment:days.since.treatment  6.8610  1    0.00881 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#with ant abundance 

m1int <- glmmTMB(adult.Chelinidea.vittiger ~ treatment * days.since.treatment + ants + preds + new.segments +  (1|year) + (1|uniID), family = "poisson", data = data)

m1int2 <- glmmTMB(adult.Chelinidea.vittiger ~ treatment * days.since.treatment +  preds + new.segments +  (1|year) + (1|uniID), family = "poisson", data = data)

AIC(m1int, m1int2)

       df      AIC
m1int   9 1978.067
m1int2  8 1980.594

#ants improve the model parsimony

check_overdispersion(m1int)

# Overdispersion test

       dispersion ratio =   0.953
  Pearson's Chi-Squared = 810.810
                p-value =   0.835

No overdispersion detected.

AIC(m1, m1int)

      df      AIC
m1     8 1979.728
m1int  9 1978.067

summary(m1int)

 Family: poisson  ( log )
Formula:          
adult.Chelinidea.vittiger ~ treatment * days.since.treatment +  
    ants + preds + new.segments + (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  1978.1   2020.9   -980.0   1960.1      851 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 7.991e-10 2.827e-05
 uniID  (Intercept) 5.145e-01 7.173e-01
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -1.571144   0.254483  -6.174 6.66e-10 ***
treatmentN                       0.442759   0.278888   1.588   0.1124    
days.since.treatment             0.023914   0.010214   2.341   0.0192 *  
ants                            -0.025774   0.012834  -2.008   0.0446 *  
preds                           -0.042242   0.080630  -0.524   0.6003    
new.segments                     0.017671   0.003667   4.818 1.45e-06 ***
treatmentN:days.since.treatment -0.036578   0.014837  -2.465   0.0137 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m1int)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: adult.Chelinidea.vittiger
                                 Chisq Df Pr(>Chisq)    
treatment                       0.7421  1    0.38899    
days.since.treatment            0.8751  1    0.34954    
ants                            4.0333  1    0.04461 *  
preds                           0.2745  1    0.60035    
new.segments                   23.2175  1  1.447e-06 ***
treatment:days.since.treatment  6.0779  1    0.01369 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

test(emtrends(m1int, pairwise ~ treatment, var = "days.since.treatment"))

$emtrends
 treatment days.since.treatment.trend     SE  df z.ratio p.value
 C                             0.0239 0.0102 Inf   2.341  0.0192
 N                            -0.0127 0.0111 Inf  -1.145  0.2521


$contrasts
 contrast estimate     SE  df z.ratio p.value
 C - N      0.0366 0.0148 Inf   2.465  0.0137

ggplot(data, aes(days.since.treatment, adult.Chelinidea.vittiger, fill = treatment, color = treatment)) +
    geom_smooth(method = "lm", size = 0.5)  +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) +
    ylab("Count of adult cactus bugs") +
    xlab("Number of days since treatment") +
  theme(legend.title = element_blank()) +
  scale_fill_manual(values = c("#F8766D", "darkcyan")) + facet_grid(~year)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

`geom_smooth()` using formula = 'y ~ x'

ggplot(data, aes(ants, herbs.all)) + 
  geom_smooth(method = "lm", color = "black") + geom_point() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  ylab("Herbivores detected per survey") +
  xlab("Number of ants per survey") +
  theme(legend.title = element_blank()) +
  scale_fill_manual(values = c("darkcyan", "#F8766D"))

`geom_smooth()` using formula = 'y ~ x'

ggplot(data, aes(pres, herbs.all)) + geom_boxplot() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) + stat_summary(fun.y=mean, geom="point", color = "red", size=4) + ylab("Herbivores detected per survey") + xlab("Ant presence")

Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
ℹ Please use the `fun` argument instead.

Ant Response

Ant Abundance

mp <- glmmTMB(ants ~ new.segments + treatment +days.since.treatment + x+  (1|year) + (1|uniID), family = "poisson", data = data)

check_overdispersion(mp)

# Overdispersion test

       dispersion ratio =    2.270
  Pearson's Chi-Squared = 1935.941
                p-value =  < 0.001

Overdispersion detected.

m4 <- glmmTMB(ants ~ new.segments + treatment +days.since.treatment + x+  (1|year) + (1|uniID), family = "nbinom2", data = data)

summary(m4)

 Family: nbinom2  ( log )
Formula:          ants ~ new.segments + treatment + days.since.treatment + x +  
    (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  2100.1   2138.2  -1042.1   2084.1      852 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 5.690e-09 7.543e-05
 uniID  (Intercept) 3.215e+00 1.793e+00
Number of obs: 860, groups:  year, 2; uniID, 140

Dispersion parameter for nbinom2 family (): 0.419 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -2.34142    0.77632  -3.016  0.00256 ** 
new.segments          0.03083    0.01003   3.073  0.00212 ** 
treatmentN            0.73149    0.36923   1.981  0.04758 *  
days.since.treatment  0.10958    0.01529   7.166 7.74e-13 ***
x                    -0.03613    0.01436  -2.515  0.01189 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(mp, m4)

   df      AIC
mp  7 2924.984
m4  8 2100.108

summary(m4)

 Family: nbinom2  ( log )
Formula:          ants ~ new.segments + treatment + days.since.treatment + x +  
    (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  2100.1   2138.2  -1042.1   2084.1      852 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 5.690e-09 7.543e-05
 uniID  (Intercept) 3.215e+00 1.793e+00
Number of obs: 860, groups:  year, 2; uniID, 140

Dispersion parameter for nbinom2 family (): 0.419 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -2.34142    0.77632  -3.016  0.00256 ** 
new.segments          0.03083    0.01003   3.073  0.00212 ** 
treatmentN            0.73149    0.36923   1.981  0.04758 *  
days.since.treatment  0.10958    0.01529   7.166 7.74e-13 ***
x                    -0.03613    0.01436  -2.515  0.01189 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m4)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: ants
                       Chisq Df Pr(>Chisq)    
new.segments          9.4437  1   0.002119 ** 
treatment             3.9249  1   0.047577 *  
days.since.treatment 51.3465  1  7.742e-13 ***
x                     6.3268  1   0.011893 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m4.int <- glmmTMB(ants ~ new.segments + treatment * days.since.treatment + x+  (1|year) + (1|uniID), family = "nbinom2", data = data) 

summary(m4.int)

 Family: nbinom2  ( log )
Formula:          ants ~ new.segments + treatment * days.since.treatment + x +  
    (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  2086.8   2129.6  -1034.4   2068.8      851 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 5.534e-09 7.439e-05
 uniID  (Intercept) 3.295e+00 1.815e+00
Number of obs: 860, groups:  year, 2; uniID, 140

Dispersion parameter for nbinom2 family (): 0.445 

Conditional model:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -1.48458    0.80781  -1.838  0.06609 .  
new.segments                     0.03153    0.01011   3.119  0.00182 ** 
treatmentN                      -1.25570    0.62854  -1.998  0.04574 *  
days.since.treatment             0.05221    0.02050   2.547  0.01088 *  
x                               -0.03561    0.01447  -2.461  0.01385 *  
treatmentN:days.since.treatment  0.11797    0.03017   3.910 9.21e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(m4, m4.int)

       df      AIC
m4      8 2100.108
m4.int  9 2086.771

#emtrends(m4.int, pairwise ~ treatment | days.since.treatment)
test(emtrends(m4.int, pairwise ~ treatment, var = "days.since.treatment"))

$emtrends
 treatment days.since.treatment.trend     SE  df z.ratio p.value
 C                             0.0522 0.0205 Inf   2.547  0.0109
 N                             0.1702 0.0222 Inf   7.680  <.0001


$contrasts
 contrast estimate     SE  df z.ratio p.value
 C - N      -0.118 0.0302 Inf  -3.910  0.0001

ggplot(data, aes(days.since.treatment, ants, fill = treatment)) + geom_smooth(method = "lm")

`geom_smooth()` using formula = 'y ~ x'

#m4 <- glmmTMB(ants ~ new.segments + treatment +days.since.treatment + x+ neigh + (1|year) + (1|uniID), family = "nbinom2", data = data)

Ant Presence

#m5 <- glmmTMB(pres ~ new.segments + x + treatment + days.since.treatment + neigh + (1|year) + (1|uniID), family = "binomial", data = data)
#summary(m5)
m5 <- glmmTMB(pres ~ new.segments + x + treatment + days.since.treatment +  (1|year) + (1|uniID), family = "binomial", data = data)
summary(m5)

 Family: binomial  ( logit )
Formula:          pres ~ new.segments + x + treatment + days.since.treatment +  
    (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
   870.6    903.9   -428.3    856.6      853 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 1.185e-09 3.442e-05
 uniID  (Intercept) 3.158e+00 1.777e+00
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -2.06004    0.80203  -2.569  0.01021 *  
new.segments          0.03006    0.01055   2.849  0.00439 ** 
x                    -0.03790    0.01479  -2.562  0.01041 *  
treatmentN            0.26986    0.37547   0.719  0.47232    
days.since.treatment  0.08332    0.01838   4.533 5.81e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_overdispersion(m5)

# Overdispersion test

 dispersion ratio = 1.009
          p-value = 0.976

No overdispersion detected.

m5.int <- glmmTMB(pres ~ new.segments + x + treatment * days.since.treatment +  (1|year) + (1|uniID), family = "binomial", data = data)

summary(m5.int)

 Family: binomial  ( logit )
Formula:          pres ~ new.segments + x + treatment * days.since.treatment +  
    (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
   858.2    896.3   -421.1    842.2      852 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 1.229e-09 3.506e-05
 uniID  (Intercept) 3.284e+00 1.812e+00
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -0.89992    0.85999  -1.046 0.295363    
new.segments                     0.03099    0.01075   2.881 0.003959 ** 
x                               -0.03892    0.01505  -2.585 0.009735 ** 
treatmentN                      -2.01644    0.72028  -2.800 0.005118 ** 
days.since.treatment             0.01176    0.02596   0.453 0.650397    
treatmentN:days.since.treatment  0.13931    0.03736   3.728 0.000193 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

test(emtrends(m5.int, pairwise ~ treatment, var = "days.since.treatment"))

$emtrends
 treatment days.since.treatment.trend     SE  df z.ratio p.value
 C                             0.0118 0.0260 Inf   0.453  0.6504
 N                             0.1511 0.0269 Inf   5.608  <.0001


$contrasts
 contrast estimate     SE  df z.ratio p.value
 C - N      -0.139 0.0374 Inf  -3.728  0.0002

ggplot(data, aes(x=days.since.treatment, y=as.numeric(data$pres) - 1, color = treatment)) + 
  geom_point(alpha=.5) +
  stat_smooth(method="glm", se=FALSE, method.args = list(family=binomial))

Warning: Use of `data$pres` is discouraged.
ℹ Use `pres` instead.
Use of `data$pres` is discouraged.
ℹ Use `pres` instead.

`geom_smooth()` using formula = 'y ~ x'

#crematogaster and solenopsis 



m4 <- glmmTMB(crempres ~ treatment + new.segments  +x + (1|year) + (1|uniID), family = "binomial", data = data)
check_overdispersion(m4)

# Overdispersion test

 dispersion ratio = 1.094
          p-value = 0.408

No overdispersion detected.

summary(m4)

 Family: binomial  ( logit )
Formula:          
crempres ~ treatment + new.segments + x + (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
   614.2    642.8   -301.1    602.2      854 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 1.683e-09 4.103e-05
 uniID  (Intercept) 1.318e+01 3.630e+00
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -1.70267    1.55405  -1.096  0.27324   
treatmentN   -1.01926    0.74321  -1.371  0.17024   
new.segments  0.06027    0.02241   2.689  0.00716 **
x            -0.07316    0.02977  -2.458  0.01399 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m5 <- glmmTMB(solpres ~ treatment + new.segments  + x + (1|year) + (1|uniID), family = "binomial", data = data)
check_overdispersion(m5)

# Overdispersion test

 dispersion ratio = 0.996
          p-value = 0.976

No overdispersion detected.

summary(m5)

 Family: binomial  ( logit )
Formula:          
solpres ~ treatment + new.segments + x + (1 | year) + (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
   297.5    326.1   -142.8    285.5      854 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 8.279e-10 2.877e-05
 uniID  (Intercept) 1.826e+00 1.351e+00
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -6.50812    1.14922  -5.663 1.49e-08 ***
treatmentN    1.69511    0.53009   3.198  0.00138 ** 
new.segments  0.02330    0.01140   2.044  0.04098 *  
x             0.01009    0.01624   0.621  0.53447    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ant models with species

#make a presence column

antslong <- mutate(antslong, presence = ifelse(antslong$abun > 0, 1, 0))


#abundances

m4p <- glmmTMB(abun ~ new.segments + treatment*ant.species + year + (1|uniID), family = "poisson", data = antslong)
check_overdispersion(m4p)

# Overdispersion test

       dispersion ratio =     4.298
  Pearson's Chi-Squared = 29489.761
                p-value =   < 0.001

Overdispersion detected.

m4 <- glmmTMB(abun ~ new.segments + treatment*ant.species + year + (1|uniID), family = "nbinom1", data = antslong)

summary(m4)

 Family: nbinom1  ( log )
Formula:          
abun ~ new.segments + treatment * ant.species + year + (1 | uniID)
Data: antslong

     AIC      BIC   logLik deviance df.resid 
  2759.1   2895.9  -1359.6   2719.1     6860 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 uniID  (Intercept) 1.22     1.105   
Number of obs: 6880, groups:  uniID, 140

Dispersion parameter for nbinom1 family ():  9.8 

Conditional model:
                                              Estimate Std. Error z value
(Intercept)                                 -4.980e+00  7.937e-01  -6.275
new.segments                                 1.118e-02  5.718e-03   1.956
treatmentN                                   9.078e-01  8.928e-01   1.017
ant.speciesCrematogaster.depilis             4.326e+00  7.141e-01   6.059
ant.speciesEndiodioctes                      9.227e-01  8.367e-01   1.103
ant.speciesEremocystus                      -6.930e-01  1.225e+00  -0.566
ant.speciesForelius.pruinosis               -6.934e-01  1.225e+00  -0.566
ant.speciesMonomorium.ergatogyna            -1.653e+01  2.743e+03  -0.006
ant.speciesMyrmecosystus.yellow             -1.748e+01  4.417e+03  -0.004
ant.speciesSolenopsis.xyloni                 1.265e+00  8.018e-01   1.577
year2024                                    -2.117e-01  2.466e-01  -0.858
treatmentN:ant.speciesCrematogaster.depilis -1.228e+00  8.799e-01  -1.396
treatmentN:ant.speciesEndiodioctes           3.442e-01  1.011e+00   0.341
treatmentN:ant.speciesEremocystus           -2.345e+01  8.700e+04   0.000
treatmentN:ant.speciesForelius.pruinosis     9.182e-01  1.396e+00   0.658
treatmentN:ant.speciesMonomorium.ergatogyna  1.583e+01  2.743e+03   0.006
treatmentN:ant.speciesMyrmecosystus.yellow   1.788e+01  4.417e+03   0.004
treatmentN:ant.speciesSolenopsis.xyloni      9.460e-01  9.609e-01   0.985
                                            Pr(>|z|)    
(Intercept)                                 3.49e-10 ***
new.segments                                  0.0505 .  
treatmentN                                    0.3092    
ant.speciesCrematogaster.depilis            1.37e-09 ***
ant.speciesEndiodioctes                       0.2701    
ant.speciesEremocystus                        0.5715    
ant.speciesForelius.pruinosis                 0.5713    
ant.speciesMonomorium.ergatogyna              0.9952    
ant.speciesMyrmecosystus.yellow               0.9968    
ant.speciesSolenopsis.xyloni                  0.1147    
year2024                                      0.3907    
treatmentN:ant.speciesCrematogaster.depilis   0.1627    
treatmentN:ant.speciesEndiodioctes            0.7334    
treatmentN:ant.speciesEremocystus             0.9998    
treatmentN:ant.speciesForelius.pruinosis      0.5108    
treatmentN:ant.speciesMonomorium.ergatogyna   0.9954    
treatmentN:ant.speciesMyrmecosystus.yellow    0.9968    
treatmentN:ant.speciesSolenopsis.xyloni       0.3248    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ct <- emmeans(m4, pairwise ~ treatment | ant.species)
xprt = print(pairs(ct), export = TRUE)
s <- do.call(rbind, xprt$summary)
s <- as.data.frame(s)
knitr::kable(s)

	contrast	estimate	SE	df	z.ratio	p.value
X	C - N	-0.908	0.89	Inf	-1.017	0.3092
X.1	C - N	0.320	0.26	Inf	1.242	0.2142
X.2	C - N	-1.252	0.56	Inf	-2.219	0.0265
X.3	C - N	22.538	87000.71	Inf	0.000	0.9998
X.4	C - N	-1.826	1.12	Inf	-1.635	0.1020
X.5	C - N	-16.739	2743.14	Inf	-0.006	0.9951
X.6	C - N	-18.791	4417.16	Inf	-0.004	0.9966
X.7	C - N	-1.854	0.47	Inf	-3.950	0.0001

m4 <- glmmTMB(presence ~ treatment*ant.species + (1|uniID), family = "binomial", data = antslong)

ct <- emmeans(m4, pairwise ~ treatment | ant.species)
ct <- emmeans(m4, pairwise ~ treatment | ant.species)
xprt = print(pairs(ct), export = TRUE)
s <- do.call(rbind, xprt$summary)
s <- as.data.frame(s)
knitr::kable(s)

	contrast	estimate	SE	df	z.ratio	p.value
X	C - N	-1.069	9.10e-01	Inf	-1.174	0.2403
X.1	C - N	0.445	2.90e-01	Inf	1.513	0.1303
X.2	C - N	-1.431	5.90e-01	Inf	-2.420	0.0155
X.3	C - N	23.426	1.47e+05	Inf	0.000	0.9999
X.4	C - N	-1.998	1.13e+00	Inf	-1.766	0.0774
X.5	C - N	-17.686	4.05e+03	Inf	-0.004	0.9965
X.6	C - N	-20.206	8.16e+03	Inf	-0.002	0.9980
X.7	C - N	-2.028	5.00e-01	Inf	-4.049	0.0001

library(mvabund)
antspp <- mvabund(data[, 22:29])
antpres <- antspp
antpres[antpres>1] = 1



plot(antspp ~ as.factor(data$treatment), cex.axis = 0.8, cex = 0.8, transformation="no")

Overlapping points were shifted along the y-axis to make them visible.


 PIPING TO 2nd MVFACTOR

mod1 <- manyglm(antspp ~ as.factor(data\(treatment) + data\)new.segments + data\(year + data\)x, family = “negative_binomial”) plot(mod1) anova(mod1, resamp = “case”, block = data\(uniID) summary(mod1, resamp = "case", block = data\)uniID) AIC(mod1)

Resource partitioning

#unchosen plants not in partition
nozer <- filter(antslong, abun > 0)

nozer <- filter(nozer, ant.species == "Solenopsis.xyloni" | ant.species == "Crematogaster.depilis" | ant.species == "Endiodioctes")



ggplot(nozer, aes(new.segments, fill = ant.species, alpha = 0.5)) + geom_density(position = "dodge") + facet_grid(~treatment)

Warning: Width not defined
ℹ Set with `position_dodge(width = ...)`

ggplot(nozer, aes(ant.species, alpha = 0.5)) + geom_bar() + facet_grid(~treatment)

ggplot(nozer, aes(treatment, alpha = 0.5)) + geom_bar() + facet_grid(~ant.species)

ggplot(nozer, aes(new.segments, fill = treatment, alpha = 0.5)) + geom_density(position = "dodge") + facet_grid(~ant.species)

Warning: Width not defined
ℹ Set with `position_dodge(width = ...)`

ggplot(nozer, aes(x, fill = treatment, alpha = 0.5)) + geom_density(position = "dodge") + facet_grid(~ant.species)

Warning: Width not defined
ℹ Set with `position_dodge(width = ...)`

#differences between species on control plants
cozer <- filter(nozer, treatment == "C")


ggplot(nozer, aes(ant.species, new.segments)) + geom_boxplot() + facet_grid(~treatment)

#ggplot(nozer, aes(ant.species, x)) + geom_boxplot() + facet_grid(~treatment)



#agregrate data at the plant level


sc$uniID <- paste(sc$PlantID, sc$year)

vis <- antslong %>% group_by(uniID, ant.species) %>% summarise(abun = sum(abun))

`summarise()` has grouped output by 'uniID'. You can override using the
`.groups` argument.

vis <- left_join(vis, sc, by = "uniID")
vis <- mutate(vis, pres = ifelse(abun > 0, 1, 0))

visno <- filter(vis, pres > 0)



nozer <- filter(visno, ant.species == "Solenopsis.xyloni" | ant.species == "Crematogaster.depilis" | ant.species == "Endiodioctes")

cozer <- filter(nozer, treatment == "C")
m1 <- glmmTMB(new.segments ~ ant.species + (1|year), family = "nbinom1", data = cozer)
#m2 <- glmmTMB(new.segments ~ ant.species + (1|year), family = "nbinom2", data = cozer)
#AIC(m1, m2)
#using nbinom1


summary(m1)

 Family: nbinom1  ( log )
Formula:          new.segments ~ ant.species + (1 | year)
Data: cozer

     AIC      BIC   logLik deviance df.resid 
   363.7    372.3   -176.9    353.7       36 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 year   (Intercept) 0.02284  0.1511  
Number of obs: 41, groups:  year, 2

Dispersion parameter for nbinom1 family (): 7.29 

Conditional model:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    3.7565     0.1337  28.087   <2e-16 ***
ant.speciesEndiodioctes       -0.6357     0.2928  -2.172   0.0299 *  
ant.speciesSolenopsis.xyloni   0.3737     0.1588   2.353   0.0186 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m1, type = 2)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: new.segments
            Chisq Df Pr(>Chisq)   
ant.species 11.31  2   0.003501 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(m1, pairwise~ ant.species)

$emmeans
 ant.species           emmean    SE  df asymp.LCL asymp.UCL
 Crematogaster.depilis   3.76 0.134 Inf      3.49      4.02
 Endiodioctes            3.12 0.301 Inf      2.53      3.71
 Solenopsis.xyloni       4.13 0.176 Inf      3.78      4.48

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast                                  estimate    SE  df z.ratio p.value
 Crematogaster.depilis - Endiodioctes         0.636 0.293 Inf   2.172  0.0761
 Crematogaster.depilis - Solenopsis.xyloni   -0.374 0.159 Inf  -2.353  0.0488
 Endiodioctes - Solenopsis.xyloni            -1.009 0.320 Inf  -3.157  0.0045

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates

m2 <- glmmTMB(new.segments ~ ant.species * treatment + (1|year), family = "nbinom1", data = nozer)
summary(m2)

 Family: nbinom1  ( log )
Formula:          new.segments ~ ant.species * treatment + (1 | year)
Data: nozer

     AIC      BIC   logLik deviance df.resid 
   882.7    903.3   -433.3    866.7       90 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 year   (Intercept) 0.01225  0.1107  
Number of obs: 98, groups:  year, 2

Dispersion parameter for nbinom1 family (): 9.14 

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)
(Intercept)                              3.77226    0.11646   32.39   <2e-16
ant.speciesEndiodioctes                 -0.58570    0.30414   -1.93   0.0541
ant.speciesSolenopsis.xyloni             0.36278    0.17206    2.11   0.0350
treatmentN                               0.05561    0.12567    0.44   0.6581
ant.speciesEndiodioctes:treatmentN       0.66789    0.33896    1.97   0.0488
ant.speciesSolenopsis.xyloni:treatmentN -0.34844    0.21939   -1.59   0.1122
                                           
(Intercept)                             ***
ant.speciesEndiodioctes                 .  
ant.speciesSolenopsis.xyloni            *  
treatmentN                                 
ant.speciesEndiodioctes:treatmentN      *  
ant.speciesSolenopsis.xyloni:treatmentN    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_overdispersion(m2)

# Overdispersion test

 dispersion ratio = 1.113
          p-value = 0.512

No overdispersion detected.

car::Anova(m2, type = 2)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: new.segments
                       Chisq Df Pr(>Chisq)  
ant.species           1.5004  2    0.47226  
treatment             0.0573  1    0.81083  
ant.species:treatment 8.1392  2    0.01708 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(m2, pairwise~ ant.species)

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 ant.species           emmean    SE  df asymp.LCL asymp.UCL
 Crematogaster.depilis   3.80 0.101 Inf      3.60      4.00
 Endiodioctes            3.55 0.177 Inf      3.20      3.90
 Solenopsis.xyloni       3.99 0.121 Inf      3.75      4.23

Results are averaged over the levels of: treatment 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast                                  estimate    SE  df z.ratio p.value
 Crematogaster.depilis - Endiodioctes         0.252 0.170 Inf   1.484  0.2988
 Crematogaster.depilis - Solenopsis.xyloni   -0.189 0.111 Inf  -1.702  0.2045
 Endiodioctes - Solenopsis.xyloni            -0.440 0.186 Inf  -2.373  0.0464

Results are averaged over the levels of: treatment 
Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates

emmeans(m2, pairwise~treatment|ant.species)

$emmeans
ant.species = Crematogaster.depilis:
 treatment emmean    SE  df asymp.LCL asymp.UCL
 C           3.77 0.116 Inf      3.54      4.00
 N           3.83 0.121 Inf      3.59      4.07

ant.species = Endiodioctes:
 treatment emmean    SE  df asymp.LCL asymp.UCL
 C           3.19 0.302 Inf      2.60      3.78
 N           3.91 0.143 Inf      3.63      4.19

ant.species = Solenopsis.xyloni:
 treatment emmean    SE  df asymp.LCL asymp.UCL
 C           4.14 0.170 Inf      3.80      4.47
 N           3.84 0.129 Inf      3.59      4.10

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
ant.species = Crematogaster.depilis:
 contrast estimate    SE  df z.ratio p.value
 C - N     -0.0556 0.126 Inf  -0.443  0.6581

ant.species = Endiodioctes:
 contrast estimate    SE  df z.ratio p.value
 C - N     -0.7235 0.312 Inf  -2.317  0.0205

ant.species = Solenopsis.xyloni:
 contrast estimate    SE  df z.ratio p.value
 C - N      0.2928 0.181 Inf   1.622  0.1048

Results are given on the log (not the response) scale.

emmeans(m2, pairwise~ant.species|treatment)

$emmeans
treatment = C:
 ant.species           emmean    SE  df asymp.LCL asymp.UCL
 Crematogaster.depilis   3.77 0.116 Inf      3.54      4.00
 Endiodioctes            3.19 0.302 Inf      2.60      3.78
 Solenopsis.xyloni       4.14 0.170 Inf      3.80      4.47

treatment = N:
 ant.species           emmean    SE  df asymp.LCL asymp.UCL
 Crematogaster.depilis   3.83 0.121 Inf      3.59      4.07
 Endiodioctes            3.91 0.143 Inf      3.63      4.19
 Solenopsis.xyloni       3.84 0.129 Inf      3.59      4.10

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
treatment = C:
 contrast                                  estimate    SE  df z.ratio p.value
 Crematogaster.depilis - Endiodioctes        0.5857 0.304 Inf   1.926  0.1314
 Crematogaster.depilis - Solenopsis.xyloni  -0.3628 0.172 Inf  -2.108  0.0881
 Endiodioctes - Solenopsis.xyloni           -0.9485 0.331 Inf  -2.868  0.0115

treatment = N:
 contrast                                  estimate    SE  df z.ratio p.value
 Crematogaster.depilis - Endiodioctes       -0.0822 0.150 Inf  -0.547  0.8478
 Crematogaster.depilis - Solenopsis.xyloni  -0.0143 0.138 Inf  -0.104  0.9941
 Endiodioctes - Solenopsis.xyloni            0.0679 0.158 Inf   0.430  0.9031

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates

ggplot(nozer, aes(ant.species, new.segments)) + geom_boxplot() + facet_grid(~treatment)

ggplot(nozer, aes(ant.species, x)) + geom_boxplot() + facet_grid(~treatment)

#ggplot(nozer, aes(ant.species, efndensity)) + geom_boxplot() + facet_grid(~treatment)

visits <- mutate(visits, solpres = ifelse(s.xyloni > 0, 1, 0))
visits <- mutate(visits, crempres = ifelse(crem > 0, 1, 0))
visits <- mutate(visits, herbpres = ifelse(herbs > 0, 1, 0))

Figures

Figure 3

ggplot(antslong, aes(x = reorder(ant.species, -abun), y = abun, fill = treatment)) + geom_bar(stat = "summary", position=position_dodge(), fun = sum) + scale_fill_manual(values = c("darkcyan", "#F8766D")) +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + ylab("Total abundance across all surveys") + xlab("Ant Species")

means <- antslong %>% group_by(ant.species, treatment) %>% summarize(mean = mean(abun), se = sd(abun)/sqrt(6880), sd = sd(abun), pressum = sum(presence))

`summarise()` has grouped output by 'ant.species'. You can override using the
`.groups` argument.

means

# A tibble: 16 × 6
# Groups:   ant.species [8]
   ant.species           treatment    mean      se     sd pressum
   <chr>                 <chr>       <dbl>   <dbl>  <dbl>   <dbl>
 1 Camponotus.sayi       C         0.00698 0.00130 0.108        2
 2 Camponotus.sayi       N         0.0233  0.00348 0.289        4
 3 Crematogaster.depilis C         1.07    0.0365  3.03       103
 4 Crematogaster.depilis N         0.793   0.0402  3.33        69
 5 Endiodioctes          C         0.0302  0.00422 0.350        5
 6 Endiodioctes          N         0.0907  0.00744 0.617       14
 7 Eremocystus           C         0.00698 0.00174 0.145        1
 8 Eremocystus           N         0       0       0            0
 9 Forelius.pruinosis    C         0.00465 0.00116 0.0964       1
10 Forelius.pruinosis    N         0.0558  0.00782 0.649        5
11 Monomorium.ergatogyna C         0       0       0            0
12 Monomorium.ergatogyna N         0.00930 0.00184 0.152        2
13 Myrmecosystus.yellow  C         0       0       0            0
14 Myrmecosystus.yellow  N         0.0163  0.00174 0.144        6
15 Solenopsis.xyloni     C         0.0837  0.0138  1.14         7
16 Solenopsis.xyloni     N         0.914   0.0501  4.15        33

str(antslong)

tibble [6,880 × 20] (S3: tbl_df/tbl/data.frame)
 $ PlantID             : chr [1:6880] "1" "1" "1" "1" ...
 $ Species             : chr [1:6880] "silvercholla" "silvercholla" "silvercholla" "silvercholla" ...
 $ date.x              : num [1:6880] 4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 ...
 $ year                : Factor w/ 2 levels "2023","2024": 1 1 1 1 1 1 1 1 1 1 ...
 $ x                   : int [1:6880] 58 58 58 58 58 58 58 58 58 58 ...
 $ y                   : int [1:6880] 36 36 36 36 36 36 36 36 36 36 ...
 $ z                   : int [1:6880] 65 65 65 65 65 65 65 65 65 65 ...
 $ treatment           : chr [1:6880] "C" "C" "C" "C" ...
 $ growth.nodes        : num [1:6880] 110 110 110 110 110 110 110 110 110 110 ...
 $ new.segments        : num [1:6880] 41 41 41 41 41 41 41 41 41 41 ...
 $ buds                : int [1:6880] 49 49 49 49 49 49 49 49 49 49 ...
 $ aborted.buds        : num [1:6880] 0 0 0 0 0 0 0 0 0 0 ...
 $ open.flowers        : num [1:6880] 0 0 0 0 0 0 0 0 0 0 ...
 $ date.y              : chr [1:6880] "4.18.2023" "4.18.2023" "4.18.2023" "4.18.2023" ...
 $ time.of.day         : chr [1:6880] "morning" "morning" "morning" "morning" ...
 $ days.since.treatment: int [1:6880] 9 9 9 9 9 9 9 9 12 12 ...
 $ uniID               : chr [1:6880] "1 2023" "1 2023" "1 2023" "1 2023" ...
 $ ant.species         : chr [1:6880] "Crematogaster.depilis" "Endiodioctes" "Eremocystus" "Monomorium.ergatogyna" ...
 $ abun                : num [1:6880] 2 0 0 0 0 0 0 0 1 0 ...
 $ presence            : num [1:6880] 1 0 0 0 0 0 0 0 1 0 ...

ggplot(antslong, aes(x = reorder(ant.species, -presence), y = presence, fill = treatment)) + geom_bar(stat = "summary", position=position_dodge(), fun = sum)

ggplot(means, aes(x = reorder(ant.species, -mean), y = mean, fill = treatment)) + geom_bar(stat = "identity", position=position_dodge()) + scale_fill_manual(values = c("darkcyan", "#F8766D")) +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + 
  ylab("Mean abundance per survey") + xlab("Ant Species") + 
  scale_x_discrete(labels = c("Crematogaster \n depilis", "Solenopsis \n xyloni", "Myrmecocystus \n flaviceps", "Forelius \n pruinosis", "Camponotus \n sayi", "Myrmecocystus \n mexicanus", "Monomorium \n ergatogyna", "Myrmecocystus \n (Eremocystus sp)")) + geom_errorbar(aes(ymin=mean, ymax=mean+se), width=.2,
                 position=position_dodge(.9))

means[means == 0] <- NA

a <- ggplot(means, aes(x = reorder(ant.species, -mean), y = mean, color = treatment)) + geom_point(position=position_dodge(.9), size = 4) + 
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + 
  ylab("Mean ant abundance\n (per plant, per survey)") + xlab("Ant Species") + 
  scale_x_discrete(labels = c("Crematogaster \n depilis", "Solenopsis \n xyloni", "Myrmecocystus \n flaviceps", "Forelius \n pruinosis", "Camponotus \n sayi", "Myrmecocystus \n mexicanus", "Monomorium \n ergatogyna", "Myrmecocystus \n (Eremocystus sp)")) + geom_errorbar(aes(ymin=ifelse(mean - se <= 0, 0, mean-se), ymax=mean+se), width=.2,
                 position=position_dodge(.9)) + scale_color_manual(name = "Treatment", labels = c("Control", "Nectar Supplementation"), values = c("darkcyan", "#F8766D")) +  theme(legend.position = c(0.85, 0.9), text=element_text(size=18), axis.text.x = element_text(face = "italic"))  + theme(axis.text.x = element_text(angle = 45, , hjust=1))

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

b <- ggplot(means, aes(x = reorder(ant.species, -pressum), y = pressum, fill = treatment)) + geom_bar(stat = "identity", position=position_dodge(.9))  +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + 
  ylab("Count of presences \nacross surveys") + xlab("Ant Species") + 
  scale_x_discrete(labels = c("Crematogaster \n depilis", "Solenopsis \n xyloni", "Myrmecocystus \n flaviceps", "Forelius \n pruinosis", "Camponotus \n sayi", "Myrmecocystus \n mexicanus", "Monomorium \n ergatogyna", "Myrmecocystus \n (Eremocystus sp)"))  + scale_fill_manual(name = "Treatment", labels = c("Control", "Nectar Supplementation"), values = c("darkcyan", "#F8766D")) +  theme(legend.position = c(0.85, 0.9), text=element_text(size=18), axis.text.x = element_text(face = "italic")) +  theme(axis.text.x = element_text(angle = 45, hjust=1))




plots <- list(a,b)
grobs <- list()
widths <- list()

for (i in 1:length(plots)){
    grobs[[i]] <- ggplotGrob(plots[[i]])
    widths[[i]] <- grobs[[i]]$widths[2:5]
}

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_bar()`).

maxwidth <- do.call(grid::unit.pmax, widths)
for (i in 1:length(grobs)){
     grobs[[i]]$widths[2:5] <- as.list(maxwidth)
}

p <- do.call("grid.arrange", c(grobs, ncol = 1))

TableGrob (2 x 1) "arrange": 2 grobs
  z     cells    name           grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (2-2,1-1) arrange gtable[layout]

Figure 4

cor.test(sc$new.segments, sc$x)


    Pearson's product-moment correlation

data:  sc$new.segments and sc$x
t = 6.6492, df = 137, p-value = 6.448e-10
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3567864 0.6102564
sample estimates:
      cor 
0.4939443

a <- ggplot(data, aes(new.segments, ants))  + geom_point(aes(x, ants), shape = 16, alpha = 0.4, position = "jitter") + ylab("Ants detected per survey")+ geom_smooth(aes(color = treatment, fill = treatment), method = "lm") + xlab("EFN Resources") +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + xlim(25, 110) + scale_color_manual(name = "Treatment", labels = c("Control", "Nectar \nSupplementation"), values = c("darkcyan", "#F8766D")) +  theme(legend.position = c(0.74, 0.9), text=element_text(size=18)) + scale_fill_manual(name = "Treatment", labels = c("Control", "Nectar \nSupplementation"), values = c("darkcyan", "#F8766D")) 
a

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 193 rows containing non-finite outside the scale range
(`stat_smooth()`).

b <- ggplot(data, aes(x, ants, fill = treatment)) + geom_smooth(aes(color = treatment), method = "lm", se = TRUE) + geom_point(aes(x, ants), shape = 16, alpha = 0.4, position = "jitter") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) + ylab("Ants detected per survey") + xlab("Cholla Width (cm)") + scale_color_manual(name = "Treatment", labels = c("Control", "Nectar \nSupplementation"), values = c("darkcyan", "#F8766D")) +  theme(legend.position = c(0.74, 0.9), text=element_text(size=18)) + scale_fill_manual(name = "Treatment", labels = c("Control", "Nectar \nSupplementation"), values = c("darkcyan", "#F8766D")) 
b

`geom_smooth()` using formula = 'y ~ x'

c <- ggplot(sc, aes(x, new.segments)) + geom_smooth(method = "lm", color = "black") + geom_point() + xlab("Cholla Width (cm)") + ylab("EFN Resources") +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text=element_text(size=18))


plots <- list(a,b,c)
grobs <- list()
widths <- list()

for (i in 1:length(plots)){
    grobs[[i]] <- ggplotGrob(plots[[i]])
    widths[[i]] <- grobs[[i]]$widths[2:5]
}

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 193 rows containing non-finite outside the scale range
(`stat_smooth()`).

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_point()`).

maxwidth <- do.call(grid::unit.pmax, widths)
for (i in 1:length(grobs)){
     grobs[[i]]$widths[2:5] <- as.list(maxwidth)
}

p <- do.call("grid.arrange", c(grobs, ncol = 3))

TableGrob (1 x 3) "arrange": 3 grobs
  z     cells    name           grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]
3 3 (1-1,3-3) arrange gtable[layout]

Figure 5

a <- ggplot(nozer, aes(new.segments, fill = ant.species, alpha = 0.5)) + geom_density(position = "dodge") + facet_grid(~treatment) + scale_fill_discrete(name = "", labels = c("Crematogaster depilis", "Myrmecocystus flaviceps", "Solenopsis xyloni")) + xlab("EFN Resources")  +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + scale_alpha(guide = "none") + theme(strip.background = element_blank(), strip.text.x = element_blank(), text=element_text(size=18)) +
  theme(legend.text = element_text(face = "italic"))
a

Warning: Width not defined
ℹ Set with `position_dodge(width = ...)`

treatment_labs <- c(C = "Control", N = "Nectar Supplementation")

b <- nozer %>% group_by(treatment, ant.species) %>%
  summarise(mean_value = mean(new.segments),
            se_value = sd(new.segments)/(sqrt(length(new.segments)))) %>%
  ggplot(aes(x = ant.species, y = mean_value, color = ant.species)) +
  geom_point(size=3) +
  geom_errorbar(aes(ymin = mean_value-se_value, ymax = mean_value+se_value, width=.2)) + facet_grid(~treatment, labeller = labeller(treatment = treatment_labs)) + scale_color_discrete(name = "", labels = c("Crematogaster depilis", "Myrmecocystus flaviceps", "Solenopsis xyloni")) +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + xlab("Ant Species") + ylab("EFN Resources") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), text=element_text(size=18)) +
  theme(legend.text = element_text(face = "italic"))

`summarise()` has grouped output by 'treatment'. You can override using the
`.groups` argument.

plots <- list(b, a)
grobs <- list()
widths <- list()

for (i in 1:length(plots)){
    grobs[[i]] <- ggplotGrob(plots[[i]])
    widths[[i]] <- grobs[[i]]$widths[2:5]
}

Warning: Width not defined
ℹ Set with `position_dodge(width = ...)`

maxwidth <- do.call(grid::unit.pmax, widths)
for (i in 1:length(grobs)){
     grobs[[i]]$widths[2:5] <- as.list(maxwidth)
}

p <- do.call("grid.arrange", c(grobs, ncol = 1))

TableGrob (2 x 1) "arrange": 2 grobs
  z     cells    name           grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (2-2,1-1) arrange gtable[layout]

SI Figure 1

a <- ggplot(sc, aes(x)) + geom_histogram() +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + ylab("Count") + xlab("Silver cholla width (cm)")

b <- ggplot(sc, aes(z)) + geom_histogram() +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + ylab("Count") + xlab("Silver cholla height (cm)")

b

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plots <- list(a, b)
grobs <- list()
widths <- list()

for (i in 1:length(plots)){
    grobs[[i]] <- ggplotGrob(plots[[i]])
    widths[[i]] <- grobs[[i]]$widths[2:5]
}

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

maxwidth <- do.call(grid::unit.pmax, widths)
for (i in 1:length(grobs)){
     grobs[[i]]$widths[2:5] <- as.list(maxwidth)
}

p <- do.call("grid.arrange", c(grobs, ncol = 1))

TableGrob (2 x 1) "arrange": 2 grobs
  z     cells    name           grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (2-2,1-1) arrange gtable[layout]

Species richness

str(data)

'data.frame':   860 obs. of  93 variables:
 $ PlantID                  : chr  "1" "1" "1" "1" ...
 $ Species                  : chr  "silvercholla" "silvercholla" "silvercholla" "silvercholla" ...
 $ date.x                   : num  4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 ...
 $ year                     : Factor w/ 2 levels "2023","2024": 1 1 1 1 1 1 1 1 1 1 ...
 $ x                        : int  58 58 58 58 58 39 39 39 39 39 ...
 $ y                        : int  36 36 36 36 36 28 28 28 28 28 ...
 $ z                        : int  65 65 65 65 65 54 54 54 54 54 ...
 $ treatment                : chr  "C" "C" "C" "C" ...
 $ growth.nodes             : num  110 110 110 110 110 97 97 97 97 97 ...
 $ new.segments             : num  41 41 41 41 41 27 27 27 27 27 ...
 $ buds                     : int  49 49 49 49 49 55 55 55 55 55 ...
 $ aborted.buds             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ open.flowers             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ date.y                   : chr  "4.18.2023" "4.21.2023" "4.26.2023" "5.1.2023" ...
 $ time.of.day              : chr  "morning" "morning" "morning" "midday" ...
 $ days.since.treatment     : int  9 12 17 22 24 9 12 17 22 24 ...
 $ CB.eggs                  : num  0 0 0 0 3 0 0 0 0 0 ...
 $ adult.Chelinidea.vittiger: num  2 0 1 0 0 1 1 0 0 0 ...
 $ nymph.Chelinidea.vittiger: num  0 0 0 0 0 0 0 0 0 0 ...
 $ Notogramma.purpuratum    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Narnia                   : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Crematogaster.depilis    : num  2 1 1 15 10 0 0 4 7 46 ...
 $ Endiodioctes             : num  0 0 0 0 0 1 0 0 0 0 ...
 $ Eremocystus              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Monomorium.ergatogyna    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Camponotus.sayi          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Solenopsis.xyloni        : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Forelius.pruinosis       : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Myrmecosystus.yellow     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ geocoris                 : num  0 0 0 0 0 0 0 0 0 0 ...
 $ orange.geocorid          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ arrhysus                 : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Largus.californicus      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ imm.cal                  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Leptoglossus.clypealus   : num  0 0 0 0 0 0 0 0 0 0 ...
 $ black.parasitoid.wasp    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ green.stink.bug          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stink.nymph              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ small.hemip              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ wingless.wasp            : num  0 0 0 0 0 0 0 0 0 0 ...
 $ long.parasitoid.wasp     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ aphid.outbreak           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ notes                    : chr  "" "" "" "" ...
 $ tiny.wasp                : num  0 0 0 0 0 0 0 0 0 0 ...
 $ salticidae               : num  0 0 0 0 0 0 0 0 0 0 ...
 $ aphid                    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ braconid                 : num  0 0 0 0 0 0 0 0 0 0 ...
 $ lacewing                 : num  0 0 0 0 0 0 0 0 0 0 ...
 $ shiny.black.fly          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ trupanea.                : num  0 0 0 0 0 0 0 0 0 0 ...
 $ red.and.black.hemipteran : num  0 0 0 0 0 0 0 0 0 0 ...
 $ other.hemip              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ eremochrysa              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ladybug                  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ pompilidae               : num  0 0 0 0 0 0 0 0 0 0 ...
 $ syrphid                  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ chrysomelidae            : num  0 0 0 0 0 0 0 0 0 0 ...
 $ orb.spider               : num  0 0 0 0 0 0 0 0 0 0 ...
 $ antracine.bombyliidae    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ caterpillar              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ lynx.spider              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ green.crab.spider        : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ichneumoidae             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ fungus.gnat              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ anthomyiidae             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ sptted.wing.fly          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ black.geocoris           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ lasioglossum             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ sml.black.coccinelid     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ X                        : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cranefly                 : num  0 0 0 0 0 0 0 0 0 0 ...
 $ other.spider             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ tetracanthus.assasin.bug : num  0 0 0 0 0 0 0 0 0 0 ...
 $ redblackmoth             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ redblackimmaturehemip    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ black.specid             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ sarcophagidae            : num  0 0 0 0 0 0 0 0 0 0 ...
 $ green.lygus              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ chrysididae              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Berytidae                : num  0 0 0 0 0 0 0 0 0 0 ...
 $ philodromidae            : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Trichodes                : num  0 0 0 0 0 0 0 0 0 0 ...
 $ brown.hemip              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ assasin.bug..apiomeris.  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ants                     : num  2 1 1 15 10 1 0 4 7 46 ...
 $ herbs.all                : num  2 0 1 0 0 1 1 0 0 0 ...
 $ preds                    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ pres                     : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 1 2 2 2 ...
 $ solpres                  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ crempres                 : num  1 1 1 1 1 0 0 1 1 1 ...
 $ uniID                    : chr  "1 2023" "1 2023" "1 2023" "1 2023" ...
 $ efndensity               : num  0.707 0.707 0.707 0.707 0.707 ...
 $ neigh                    : num  NA NA NA NA NA NA NA NA NA NA ...

ag <- data %>% select(PlantID, year, 22:29) 
ag <- ag %>% group_by(PlantID, year) %>% summarise(across(1:8, ~ sum(.x)))

`summarise()` has grouped output by 'PlantID'. You can override using the
`.groups` argument.

library(vegan)

sp <- specnumber(ag[,3:10])
ag$richness <- sp



ggplot(ag, aes(year, richness)) + geom_boxplot()

ag$uniID <- paste(ag$PlantID, ag$year)
sc$uniID <- paste(sc$PlantID, sc$year)
ag$year <- as.integer(ag$year)
ag <- left_join(ag, sc, by = "uniID")
str(ag)

tibble [140 × 44] (S3: tbl_df/tbl/data.frame)
 $ PlantID.x            : chr [1:140] "1" "1" "10" "10" ...
 $ year.x               : int [1:140] 1 2 1 2 1 2 1 2 1 2 ...
 $ Crematogaster.depilis: num [1:140] 29 15 0 0 47 0 0 0 13 2 ...
 $ Endiodioctes         : num [1:140] 0 0 0 0 0 0 0 0 0 0 ...
 $ Eremocystus          : num [1:140] 0 0 0 0 0 0 0 0 0 3 ...
 $ Monomorium.ergatogyna: num [1:140] 0 0 0 0 0 0 0 0 0 0 ...
 $ Camponotus.sayi      : num [1:140] 0 0 0 0 0 0 0 0 0 0 ...
 $ Solenopsis.xyloni    : num [1:140] 0 3 0 0 0 0 0 0 0 0 ...
 $ Forelius.pruinosis   : num [1:140] 0 0 0 0 0 0 0 0 0 0 ...
 $ Myrmecosystus.yellow : num [1:140] 0 0 0 0 0 0 0 0 0 0 ...
 $ richness             : int [1:140] 1 2 0 0 1 0 0 0 1 2 ...
 $ uniID                : chr [1:140] "1 2023" "1 2024" "10 2023" "10 2024" ...
 $ PlantID.y            : chr [1:140] "1" "1" "10" "10" ...
 $ Species              : chr [1:140] "silvercholla" "silvercholla" "silvercholla" "silvercholla" ...
 $ date                 : num [1:140] 4.25 NA 4.25 NA 4.25 NA 4.25 NA 4.25 NA ...
 $ year.y               : int [1:140] 2023 2024 2023 2024 2023 2024 2023 2024 2023 2024 ...
 $ x                    : int [1:140] 58 58 65 50 48 43 55 76 40 52 ...
 $ y                    : int [1:140] 36 35 38 31 40 33 50 70 38 29 ...
 $ z                    : int [1:140] 65 70 41 43 44 36 67 62 52 64 ...
 $ treatment            : chr [1:140] "C" "C" "N" "N" ...
 $ growth.nodes         : int [1:140] 110 NA 96 NA 98 NA 92 NA 166 NA ...
 $ new.segments         : int [1:140] 41 50 25 31 49 31 22 83 72 17 ...
 $ buds                 : int [1:140] 49 101 30 81 11 58 40 93 41 36 ...
 $ aborted.buds         : int [1:140] NA NA NA NA NA NA 1 NA 2 NA ...
 $ open.flowers         : int [1:140] NA NA NA NA NA NA NA NA NA NA ...
 $ LT                   : int [1:140] 2 NA 1 NA NA NA NA NA 2 NA ...
 $ AS                   : int [1:140] 1 NA NA NA NA NA 3 NA 2 NA ...
 $ dead                 : int [1:140] 2 NA 4 NA 5 NA 6 NA NA NA ...
 $ CR                   : int [1:140] NA NA 1 NA 1 NA NA NA NA NA ...
 $ EC                   : int [1:140] 1 NA NA NA 1 NA NA NA 2 NA ...
 $ YS                   : int [1:140] NA NA NA NA NA NA NA NA NA NA ...
 $ ecoop                : int [1:140] 1 NA NA NA NA NA NA NA NA NA ...
 $ TM                   : int [1:140] 1 NA NA NA NA NA NA NA NA NA ...
 $ K                    : int [1:140] 1 NA NA NA NA NA NA NA NA NA ...
 $ BB                   : int [1:140] NA NA NA NA NA NA NA NA NA NA ...
 $ CE                   : int [1:140] NA NA NA NA NA NA 1 NA NA NA ...
 $ CA                   : int [1:140] NA NA NA NA NA NA NA NA NA NA ...
 $ ambsal               : int [1:140] NA NA NA NA NA NA NA NA NA NA ...
 $ acacia               : int [1:140] NA NA NA NA NA NA NA NA NA NA ...
 $ HH                   : int [1:140] NA NA NA NA NA NA NA NA NA NA ...
 $ unknown              : int [1:140] NA NA NA NA NA NA NA NA NA NA ...
 $ tetradynia           : int [1:140] NA NA NA NA NA NA NA NA NA NA ...
 $ efn.neighbor         : int [1:140] NA NA NA NA NA NA NA NA NA NA ...
 $ shrub.neighbor       : int [1:140] NA NA NA NA NA NA NA NA NA NA ...

nozer1 <- filter(ag, richness > 0)

nozer1 %>% group_by(treatment) %>% summarize(mean = mean(richness))

# A tibble: 2 × 2
  treatment  mean
  <chr>     <dbl>
1 C          1.13
2 N          1.39

m1 <- glm(richness ~ treatment + new.segments + x, family = "poisson", data = nozer1)
check_overdispersion(m1)

# Overdispersion test

       dispersion ratio =  0.207
  Pearson's Chi-Squared = 17.764
                p-value =      1

No overdispersion detected.

summary(m1)


Call:
glm(formula = richness ~ treatment + new.segments + x, family = "poisson", 
    data = nozer1)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.064799   0.395392  -0.164    0.870
treatmentN    0.199520   0.192360   1.037    0.300
new.segments  0.003013   0.004796   0.628    0.530
x             0.001100   0.007584   0.145    0.885

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 17.123  on 89  degrees of freedom
Residual deviance: 15.268  on 86  degrees of freedom
AIC: 217.67

Number of Fisher Scoring iterations: 4

Predators

m5 <- glmmTMB(preds ~ new.segments + treatment  + solpres + x + (1|year) + (1|uniID), family = "poisson", data = data)
summary(m5)

 Family: poisson  ( log )
Formula:          
preds ~ new.segments + treatment + solpres + x + (1 | year) +      (1 | uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  1023.2   1056.5   -504.6   1009.2      853 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 1.346e-09 3.669e-05
 uniID  (Intercept) 7.558e-01 8.693e-01
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.812350   0.437683  -4.141 3.46e-05 ***
new.segments -0.006483   0.006339  -1.023    0.306    
treatmentN    0.201255   0.221673   0.908    0.364    
solpres       0.344011   0.289419   1.189    0.235    
x             0.001643   0.008154   0.202    0.840    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_overdispersion(m5)

# Overdispersion test

       dispersion ratio =   0.832
  Pearson's Chi-Squared = 710.010
                p-value =       1

No overdispersion detected.

m5 <- glmmTMB(preds ~ new.segments + treatment  + pres + x + (1|year) + (1|uniID), family = "poisson", data = data)
summary(m5)

 Family: poisson  ( log )
Formula:          
preds ~ new.segments + treatment + pres + x + (1 | year) + (1 |      uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  1023.7   1057.0   -504.8   1009.7      853 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 1.458e-09 3.819e-05
 uniID  (Intercept) 7.582e-01 8.707e-01
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.768558   0.443323  -3.989 6.63e-05 ***
new.segments -0.005646   0.006390  -0.883    0.377    
treatmentN    0.231290   0.220516   1.049    0.294    
pres1        -0.173819   0.193022  -0.901    0.368    
x             0.001105   0.008201   0.135    0.893    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m5 <- glmmTMB(preds ~ new.segments + treatment  + ants + x + (1|year) + (1|uniID), family = "poisson", data = data)
summary(m5)

 Family: poisson  ( log )
Formula:          
preds ~ new.segments + treatment + ants + x + (1 | year) + (1 |      uniID)
Data: data

     AIC      BIC   logLik deviance df.resid 
  1024.4   1057.7   -505.2   1010.4      853 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 year   (Intercept) 1.556e-09 3.945e-05
 uniID  (Intercept) 7.815e-01 8.840e-01
Number of obs: 860, groups:  year, 2; uniID, 140

Conditional model:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.822683   0.442948  -4.115 3.87e-05 ***
new.segments -0.005958   0.006458  -0.922    0.356    
treatmentN    0.232423   0.222726   1.044    0.297    
ants         -0.006143   0.018739  -0.328    0.743    
x             0.001517   0.008276   0.183    0.855    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pitfalls

pts <- read.csv("Raw Data/cacti_pitfalls.csv")
unique(pts$plant.species)

[1] "buckhorn"       "silver cholla"  "pencil"         "silver cholla "

pts$plant.species <- gsub(" ", "", pts$plant.species)

#clean silver cholla data
pts <- filter(pts, plant.species == "silvercholla")
pts$PlantID <- gsub("S", "", pts$PlantID)
pts$PlantID <- gsub("C", "", pts$PlantID)
unique(pts$samplingperiod)

[1] "summer"      "21-Apr"      "last spring" "spring"

pts <- filter(pts, samplingperiod == "last spring" | samplingperiod == "spring")

unique(pts$family)

 [1] "Bethylidae"              "cali plant bug"         
 [3] "chalcid"                 "Cicadellidae"           
 [5] "Formicidae"              "Geocoridae"             
 [7] "Ichneumoidae/Braconidae" "immature"               
 [9] "immature plant bug "     "Melyridae"              
[11] "Mutillidae"              "Oxyopidae"              
[13] "Phoridae"                "Platygastridae"         
[15] "pompilidae"              "Pompilidae"             
[17] "Salticidae"              "Sarcophagidae"          
[19] "shiny black"             "spotted eye?"           
[21] "Syrphidae"               "Tachinidae"             
[23] "Tettigonidae?"           "wingless "              
[25] ""

pts$uniID <- paste(pts$PlantID, pts$year)




pts <- filter(pts, family == "Formicidae")
pts$count <- as.numeric(pts$count)

unique(pts$genus)

[1] "Camponotus.fragilis" "Crematogaster"       "Endiodioctes"       
[4] "Forelius"            "Monomorium"          "Myrmecocystus"      
[7] "Pheidole.desertorum" "Solenopsis"

pts <- group_by(pts, uniID, family, genus) %>% summarise(sum = sum(count))

`summarise()` has grouped output by 'uniID', 'family'. You can override using
the `.groups` argument.

#make it wide
pts_wide <- pts %>% pivot_wider(names_from = "genus", values_from = "sum")
sc <- select(sc, 1:13)
sc$uniID <- paste(sc$PlantID, sc$year)
pts_wide <- left_join(sc, pts_wide, by = "uniID")
sc$uniID <- paste(sc$PlantID, sc$year)

pts_wide[,16:23][is.na(pts_wide[,16:23])] <- 0


#get rid of 18 2023 and others with issues

pts_wide <- filter(pts_wide, uniID != "18 2023")
pts_wide <- filter(pts_wide, uniID != "24 2023" & uniID != "59 2023")


nopits <- read.csv("Raw Data/no_pits_2024.csv")
nopits$uniID <- paste(nopits$PlantID, nopits$year)

pts_wide <- pts_wide %>% 
  filter(! uniID %in% nopits$uniID)






#clean other species data

pts <- read.csv("Raw Data/cacti_pitfalls.csv")
unique(pts$plant.species)

[1] "buckhorn"       "silver cholla"  "pencil"         "silver cholla "

pts$plant.species <- gsub(" ", "", pts$plant.species)


other <- filter(pts, plant.species != "silvercholla")
other$PlantID <- gsub("B", "", other$PlantID)
other$PlantID <- gsub("R", "", other$PlantID)
unique(other$samplingperiod)

[1] "later spring" "26-Apr"       "21-Apr"       "spring"       "summer"

other <- filter(other, samplingperiod == "later spring" | samplingperiod == "spring")
unique(other$family)

 [1] "Acrididae"     "Anthomyiidae"  "Bethylidae"    "Braconidae"   
 [5] "Cecidomyiidae" "Chrysididae"   "Chyrsididae"   "Cicadellidae" 
 [9] "Coccinelidae"  "Formicidae"    "Geocoridae"    "immature"     
[13] "Largidae"      "Melyridae"     "Phoridae"      "pompilidae"   
[17] "Pyrochroidae"  "Sarcophagidae" "small wasp"    "Syrphidae"    
[21] "tiny"          "Tipulidae"     "uliiidae"      ""

unique(other$genus)

 [1] ""                         "larva"                   
 [3] "Camponotus.fragilis"      "Camponotus.semitestaceus"
 [5] "Crematogaster"            "Eremocystus"             
 [7] "Forelius"                 "Myrmecocystus"           
 [9] "Endiodioctes"             "Pheidole.desertorum"     
[11] "Solenopsis"               "Geocoris"                
[13] "red"                      "Largus"                  
[15] "Dastyinae"

other <- filter(other, family == "Formicidae")
other$count <- as.numeric(other$count)
other$uniID <- paste(other$PlantID, other$year, other$plant.species)
other <- group_by(other, uniID, plant.species, family, genus) %>% summarise(sum = sum(count))

`summarise()` has grouped output by 'uniID', 'plant.species', 'family'. You can
override using the `.groups` argument.

other_wide <- other %>% pivot_wider(names_from = "genus", values_from = "sum")


othpt <- filter(plt, Species == "pencilcholla" | Species == "buckhorncholla")

othpt$Species <- gsub("cholla", "", othpt$Species)

othpt$uniID <- paste(othpt$PlantID, othpt$year, othpt$Species)

othpt <- select(othpt, 1:13, 31)

other_wide <- left_join(othpt, other_wide, by = "uniID")
other_wide <- select(other_wide, 1:14, 17:25)
other_wide[,15:23][is.na(other_wide[,15:23])] <- 0



#join plant info to others
#add the dataframes together, make wide and then make long again

#other_wide <- select(other_wide, -plant.species)
other_wide$uniID <- paste(other_wide$PlantID, other_wide$year, other_wide$Species)
#other_wide <- select(other_wide, -family)
pts_wide <- select(pts_wide, -family)

#make both wides long
other_long <- pivot_longer(other_wide, 15:23, names_to = "ant.species", values_to = "abun")

pts_long <- pivot_longer(pts_wide, 15:22, names_to = "ant.species", values_to = "abun")


all_long <- rbind(pts_long, other_long)


all_wide <- pivot_wider(all_long, names_from = "ant.species", values_from = "abun")

all_wide[,15:24][is.na(all_wide[,15:24])] <- 0


all_wide <- all_wide %>% rowwise() %>% mutate(abun = sum(c_across(15:24)))

sp <- specnumber(all_wide[,15:24])
all_wide$richness <- sp

all_wide <- mutate(all_wide, anypres = ifelse(abun > 0, 1, 0))
all_wide <- mutate(all_wide, crempres = ifelse(Crematogaster > 0, 1, 0))
all_wide <- mutate(all_wide, solpres = ifelse(Solenopsis > 0, 1, 0))

Analysis

#control chollas only
sc_wide <- filter(all_wide, Species == "silvercholla")
sc_wide <- filter(sc_wide, PlantID != "18trap")
sc_wide <- filter(sc_wide, treatment != "N")

#sc_wide$treatment <- gsub("", "C", sc_wide$treatment)


allcontrol <- filter(all_wide, treatment != "N")
control <- filter(all_wide, treatment != "N" & Species != "pencil")


#the three species together
m3p <- glmmTMB(abun ~ Species + x +(1|year), family = "poisson", data = allcontrol)
m3 <- glmmTMB(abun ~ Species + x + (1|year), family = "nbinom1", data = allcontrol)
m32 <- glmmTMB(abun ~ Species+ x + (1|year), family = "nbinom2", data = allcontrol)

AIC(m3p, m3, m32)

    df       AIC
m3p  5 1395.4221
m3   6  662.8772
m32  6  671.3076

#nbinom1 works best 

car::Anova(m3)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: abun
          Chisq Df Pr(>Chisq)    
Species 21.6822  2  1.958e-05 ***
x        2.3993  1     0.1214    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(m3, pairwise ~ Species)

$emmeans
 Species      emmean    SE  df asymp.LCL asymp.UCL
 buckhorn      1.169 0.236 Inf     0.706      1.63
 pencil       -0.259 0.397 Inf    -1.037      0.52
 silvercholla  1.858 0.244 Inf     1.380      2.34

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast                estimate    SE  df z.ratio p.value
 buckhorn - pencil          1.428 0.379 Inf   3.770  0.0005
 buckhorn - silvercholla   -0.688 0.302 Inf  -2.276  0.0592
 pencil - silvercholla     -2.116 0.456 Inf  -4.637  <.0001

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates

#buckhorn to silver only
m2 <- glmmTMB(richness ~ Species *  new.segments + (1|year), family = "poisson", data = control)
summary(m2)

 Family: poisson  ( log )
Formula:          richness ~ Species * new.segments + (1 | year)
Data: control

     AIC      BIC   logLik deviance df.resid 
   284.4    298.3   -137.2    274.4      113 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 year   (Intercept) 0.0561   0.2369  
Number of obs: 118, groups:  year, 2

Conditional model:
                                  Estimate Std. Error z value Pr(>|z|)  
(Intercept)                      -0.606235   0.306802  -1.976   0.0482 *
Speciessilvercholla               0.435207   0.399710   1.089   0.2762  
new.segments                      0.003751   0.003917   0.958   0.3383  
Speciessilvercholla:new.segments  0.001400   0.007704   0.182   0.8558  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(m2, pairwise ~ Species)

NOTE: Results may be misleading due to involvement in interactions

$emmeans
 Species       emmean    SE  df asymp.LCL asymp.UCL
 buckhorn     -0.4421 0.232 Inf    -0.897    0.0131
 silvercholla  0.0544 0.213 Inf    -0.362    0.4712

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast                estimate    SE  df z.ratio p.value
 buckhorn - silvercholla   -0.496 0.204 Inf  -2.438  0.0148

Results are given on the log (not the response) scale.

ggplot(allcontrol, aes(Species, abun)) + geom_boxplot()

ggplot(allcontrol, aes(Species, abun)) + geom_boxplot() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) + stat_summary(fun.y=mean, geom="point", color = "red", size=4) + ylab("Ant abundance within pitfall trap") + xlab("Cholla Species")

ggplot(allcontrol, aes(Species, richness)) + geom_boxplot() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) + stat_summary(fun.y=mean, geom="point", color = "red", size=4) + ylab("Ant species richness within pitfall trap") + xlab("Cholla Species")

c_long <- filter(all_long, treatment != "N")
c_means <- c_long %>% group_by(Species, ant.species) %>% summarize(means = mean(abun), se = sd(abun)/sqrt(1463))

`summarise()` has grouped output by 'Species'. You can override using the
`.groups` argument.

c_means <- filter(c_means, ant.species != "Eremocystus")
c_means$Species <- factor(c_means$Species, levels = c("silvercholla", "buckhorn", "pencil"))

Figure 6

c_means[c_means == 0] <- NA
c_means <- filter(c_means, ant.species != "Eremocystus")
ggplot(c_means, aes(reorder(ant.species, -means), means, fill = Species)) + geom_point(aes(color = Species), position = position_dodge(.9), size = 4) + geom_errorbar(aes(ymin=ifelse(means - se <= 0, 0, means-se), ymax=means+se, color = Species), width=.2,
                 position=position_dodge(.9)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), , text=element_text(size=16)) + xlab("Ant Species") + ylab("Mean ant abundance per pitfall trap") + scale_x_discrete(labels = c("Crematogaster \n depilis",
 "Myrmecocystus \n mexicanus","Solenopsis \n xyloni", "Pheidole \n desertorum",  "Camponotus \n fragilis", "Monomorium \n ergatogyna", "Camponotus \n semitestacaeus", "Myrmecocystus \n flaviceps",  "Forelius \n pruinosis"))  + scale_colour_brewer(palette = "Dark2", labels = c("Silver Cholla", "Buckhorn Cholla", "Pencil Cholla")) + labs(color="Plant Species", fill="Plant Species") + scale_fill_brewer(palette = "Dark2", labels = c("Silver Cholla", "Buckhorn Cholla", "Pencil Cholla")) + theme(legend.position = c(0.8, 0.8)) + theme(axis.text.x = element_text(angle = 45, hjust=1, face = "italic"))

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Plant fitness

Seeds per fruit

#read in plant fitness data
fit <- read.csv("Raw Data/cholla_seedset.csv")
fit$PlantID <- as.character(fit$PlantID)
unique(fit$PlantID)

 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15"
[16] "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "30"
[31] "32" "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45" "46"
[46] "47" "48" "49" "50" "51" "52" "53" "54" "55" "56" "57" "58" "59" "60"

fit <- filter(fit, total.seeds != "na" & total.seeds != " ")

fit$total.seeds <- as.numeric(fit$total.seeds)
#blanks mean no non-viable seeds
fit$non.viable[is.na(fit$non.viable)] <- 0
str(fit)

'data.frame':   191 obs. of  7 variables:
 $ PlantID    : chr  "1" "1" "1" "2" ...
 $ status     : chr  "full" "full" "full" "partial" ...
 $ rep        : int  1 2 3 1 2 3 1 2 3 1 ...
 $ total.seeds: num  46 39 45 2 3 12 31 44 17 53 ...
 $ non.viable : num  0 18 5 0 0 0 25 0 0 3 ...
 $ X          : chr  "" "21" "" "" ...
 $ X.1        : chr  "" "" "" "" ...

fit <- mutate(fit, viable = total.seeds - non.viable)



#read in plant fitness data
fit24 <- read.csv("Raw Data/cholla_seedset_2024.csv")
fit24$PlantID <- as.character(fit24$PlantID)
unique(fit24$PlantID)

 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "16"
[16] "15" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "30"
[31] "31" "32" "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45"
[46] "46" "47" "48" "49" "50" "51" "52" "53" "54" "55" "56" "57" "58" "59" "60"
[61] "61" "62" "63" "64" "65" "66" "67" "68" "69" "70" "71" "72" "73" "74" "75"
[76] "76" "77" "78" "79" "80"

fit24 <- filter(fit24, total.seeds != "na" & total.seeds != " ")

fit24$total.seeds <- as.numeric(fit24$total.seeds)
#blanks mean no non-viable seeds
fit24$non.viable[is.na(fit24$non.viable)] <- 0
str(fit24)

'data.frame':   213 obs. of  5 variables:
 $ PlantID    : chr  "1" "1" "1" "2" ...
 $ status     : chr  "full" "full" "full" "full" ...
 $ rep        : int  1 2 3 1 2 1 2 3 1 2 ...
 $ total.seeds: num  32 44 47 25 37 37 30 36 41 59 ...
 $ non.viable : num  0 0 1 6 3 4 6 0 1 0 ...

# calculate viable seeds
fit24 <- mutate(fit24, viable = total.seeds - non.viable)

fit <- select(fit, 1:5, 8)
fit$year <- 2023
fit24$year <- 2024
fit <- rbind(fit, fit24)

fit$uniID <- paste(fit$PlantID, fit$year)
sc$uniID <- paste(sc$PlantID, sc$year)
seed <- left_join(sc, fit, by = "uniID")

seed <- unique(seed)
seed <- filter(seed, status == "full")



# add in aggregation 
seed <- left_join(seed, ag, by = "uniID")


#data$pres <- as.numeric(data$pres)

visits$uniID <- paste(visits$PlantID, visits$year)

seed <- left_join(seed, visits, by = "uniID")


seed <- seed %>% select(1:11, 14:21, 24:32, 68:81)

# seed$uniID <- paste(seed$PlantID.x.x, seed$year.x.x)
# seed$crepres <- ifelse(seed$Crematogaster.depilis > 0, 1, 0)
# seed$solpres <- ifelse(seed$Solenopsis.xyloni > 0, 1, 0)
# 
seed$antpres <- ifelse(seed$ants > 0, 1, 0)


str(seed)

'data.frame':   377 obs. of  43 variables:
 $ PlantID.x.x          : chr  "1" "1" "1" "3" ...
 $ Species.x            : chr  "silvercholla" "silvercholla" "silvercholla" "silvercholla" ...
 $ date.x               : num  4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 4.25 ...
 $ year.x.x             : int  2023 2023 2023 2023 2023 2023 2023 2023 2023 2023 ...
 $ x.x                  : int  58 58 58 70 70 70 69 69 82 82 ...
 $ y.x                  : int  36 36 36 60 60 60 40 40 59 59 ...
 $ z.x                  : int  65 65 65 101 101 101 57 57 70 70 ...
 $ treatment.x          : chr  "C" "C" "C" "C" ...
 $ growth.nodes.x       : int  110 110 110 194 194 194 179 179 98 98 ...
 $ new.segments.x       : int  41 41 41 46 46 46 58 58 36 36 ...
 $ buds.x               : int  49 49 49 64 64 64 31 31 17 17 ...
 $ uniID                : chr  "1 2023" "1 2023" "1 2023" "3 2023" ...
 $ PlantID.y.x          : chr  "1" "1" "1" "3" ...
 $ status               : chr  "full" "full" "full" "full" ...
 $ rep                  : int  1 2 3 1 2 3 1 2 1 2 ...
 $ total.seeds          : num  46 39 45 31 44 17 53 62 31 33 ...
 $ non.viable           : num  0 18 5 25 0 0 3 2 0 0 ...
 $ viable               : num  46 21 40 6 44 17 50 60 31 33 ...
 $ year.y.x             : num  2023 2023 2023 2023 2023 ...
 $ Crematogaster.depilis: num  29 29 29 5 5 5 0 0 6 6 ...
 $ Endiodioctes         : num  0 0 0 0 0 0 3 3 0 0 ...
 $ Eremocystus          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Monomorium.ergatogyna: num  0 0 0 0 0 0 0 0 0 0 ...
 $ Camponotus.sayi      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Solenopsis.xyloni    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Forelius.pruinosis   : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Myrmecosystus.yellow : num  0 0 0 0 0 0 0 0 0 0 ...
 $ richness             : int  1 1 1 1 1 1 1 1 1 1 ...
 $ new.segments         : num  41 41 41 46 46 46 58 58 36 36 ...
 $ buds                 : int  49 49 49 64 64 64 31 31 17 17 ...
 $ treatment            : chr  "C" "C" "C" "C" ...
 $ ants                 : num  29 29 29 5 5 5 3 3 6 6 ...
 $ bugs                 : num  3 3 3 0 0 0 7 7 8 8 ...
 $ herbs                : num  3 3 3 0 0 0 9 9 8 8 ...
 $ s.xyloni             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ attend               : num  5 5 5 1 1 1 1 1 2 2 ...
 $ crem                 : num  29 29 29 5 5 5 0 0 6 6 ...
 $ cremattend           : num  5 5 5 1 1 1 0 0 2 2 ...
 $ preds                : num  0 0 0 0 0 0 1 1 1 1 ...
 $ solpres              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ crempres             : num  1 1 1 1 1 1 0 0 1 1 ...
 $ herbpres             : num  1 1 1 0 0 0 1 1 1 1 ...
 $ antpres              : num  1 1 1 1 1 1 1 1 1 1 ...

mean(seed$viable, na.rm = TRUE)

[1] 34.56757

sd(seed$viable, na.rm = TRUE)

[1] 16.91593

test <- seed
test$viable[is.na(test$viable)] <- 0

m1 <- glmmTMB(viable ~ as.factor(antpres)  + new.segments + treatment.x + x.x + (1|year.x.x) + (1|uniID), family = "nbinom1", data = seed)

summary(m1)

 Family: nbinom1  ( log )
Formula:          
viable ~ as.factor(antpres) + new.segments + treatment.x + x.x +  
    (1 | year.x.x) + (1 | uniID)
Data: seed

     AIC      BIC   logLik deviance df.resid 
  3183.7   3215.0  -1583.8   3167.7      362 

Random effects:

Conditional model:
 Groups   Name        Variance  Std.Dev. 
 year.x.x (Intercept) 6.550e-10 2.559e-05
 uniID    (Intercept) 2.094e-01 4.576e-01
Number of obs: 370, groups:  year.x.x, 2; uniID, 131

Dispersion parameter for nbinom1 family (): 6.45 

Conditional model:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)          3.7129085  0.2056535  18.054   <2e-16 ***
as.factor(antpres)1  0.1137809  0.1023549   1.112    0.266    
new.segments         0.0006048  0.0026036   0.232    0.816    
treatment.xN        -0.1536927  0.0959584  -1.602    0.109    
x.x                 -0.0051363  0.0036162  -1.420    0.156    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m1.nb2 <- glmmTMB(viable ~ as.factor(antpres)  +   new.segments + treatment.x + (1|year.x.x) + (1|uniID), family = "nbinom2", data = seed)

#nbinom1 is 100 AIC better than 2

AIC(m1, m1.nb2)

       df      AIC
m1      8 3183.692
m1.nb2  7 3282.309

summary(m1)

 Family: nbinom1  ( log )
Formula:          
viable ~ as.factor(antpres) + new.segments + treatment.x + x.x +  
    (1 | year.x.x) + (1 | uniID)
Data: seed

     AIC      BIC   logLik deviance df.resid 
  3183.7   3215.0  -1583.8   3167.7      362 

Random effects:

Conditional model:
 Groups   Name        Variance  Std.Dev. 
 year.x.x (Intercept) 6.550e-10 2.559e-05
 uniID    (Intercept) 2.094e-01 4.576e-01
Number of obs: 370, groups:  year.x.x, 2; uniID, 131

Dispersion parameter for nbinom1 family (): 6.45 

Conditional model:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)          3.7129085  0.2056535  18.054   <2e-16 ***
as.factor(antpres)1  0.1137809  0.1023549   1.112    0.266    
new.segments         0.0006048  0.0026036   0.232    0.816    
treatment.xN        -0.1536927  0.0959584  -1.602    0.109    
x.x                 -0.0051363  0.0036162  -1.420    0.156    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m1)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: viable
                    Chisq Df Pr(>Chisq)
as.factor(antpres) 1.2357  1     0.2663
new.segments       0.0540  1     0.8163
treatment.x        2.5653  1     0.1092
x.x                2.0174  1     0.1555

ggplot(seed, aes(treatment.x, viable, fill = as.factor(antpres))) + geom_boxplot()

Warning: Removed 7 rows containing non-finite outside the scale range
(`stat_boxplot()`).

m2 <- glmmTMB(viable ~ crempres  + treatment.x +  new.segments.x +  (1|year.x.x) + (1|uniID), family = "nbinom1", data = seed)

Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
NA/NaN function evaluation

summary(m2)

 Family: nbinom1  ( log )
Formula:          
viable ~ crempres + treatment.x + new.segments.x + (1 | year.x.x) +  
    (1 | uniID)
Data: seed

     AIC      BIC   logLik deviance df.resid 
  3156.7   3184.0  -1571.3   3142.7      360 

Random effects:

Conditional model:
 Groups   Name        Variance  Std.Dev. 
 year.x.x (Intercept) 5.150e-10 2.269e-05
 uniID    (Intercept) 2.063e-01 4.542e-01
Number of obs: 367, groups:  year.x.x, 2; uniID, 130

Dispersion parameter for nbinom1 family (): 6.52 

Conditional model:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)     3.468813   0.125150  27.717   <2e-16 ***
crempres        0.231155   0.097230   2.377   0.0174 *  
treatment.xN   -0.113751   0.095048  -1.197   0.2314    
new.segments.x -0.001439   0.002336  -0.616   0.5378    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m2)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: viable
                Chisq Df Pr(>Chisq)  
crempres       5.6521  1    0.01743 *
treatment.x    1.4323  1    0.23139  
new.segments.x 0.3797  1    0.53778  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m3 <- glmmTMB(viable ~ solpres + treatment.x +  new.segments.x +  (1|year.x.x) + (1|uniID), family = "nbinom1", data = seed)
summary(m3)

 Family: nbinom1  ( log )
Formula:          
viable ~ solpres + treatment.x + new.segments.x + (1 | year.x.x) +  
    (1 | uniID)
Data: seed

     AIC      BIC   logLik deviance df.resid 
  3161.5   3188.9  -1573.8   3147.5      360 

Random effects:

Conditional model:
 Groups   Name        Variance  Std.Dev. 
 year.x.x (Intercept) 3.859e-10 1.964e-05
 uniID    (Intercept) 2.179e-01 4.668e-01
Number of obs: 367, groups:  year.x.x, 2; uniID, 130

Dispersion parameter for nbinom1 family (): 6.52 

Conditional model:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     3.5389300  0.1232266  28.719   <2e-16 ***
solpres        -0.1074731  0.1281701  -0.839    0.402    
treatment.xN   -0.1228789  0.0984552  -1.248    0.212    
new.segments.x -0.0003312  0.0024437  -0.136    0.892    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m3)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: viable
                Chisq Df Pr(>Chisq)
solpres        0.7031  1     0.4017
treatment.x    1.5577  1     0.2120
new.segments.x 0.0184  1     0.8922

ggplot(seed, aes(as.factor(crempres), viable, fill = treatment.x)) + geom_boxplot()

Warning: Removed 7 rows containing non-finite outside the scale range
(`stat_boxplot()`).

ggplot(seed, aes(as.factor(crempres), viable)) + geom_boxplot() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) + stat_summary(fun.y=mean, geom="point", color = "red", size=4) + ylab("Seeds per fruit") + xlab("Crematogaster attendance")

Warning: Removed 7 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Warning: Removed 7 rows containing non-finite outside the scale range
(`stat_summary()`).

group_by(seed, crempres) %>% summarise(mean = mean(viable, na.rm = TRUE))

# A tibble: 2 × 2
  crempres  mean
     <dbl> <dbl>
1        0  32.5
2        1  37.8

t.test(seed$viable ~ seed$year.x.x)


    Welch Two Sample t-test

data:  seed$viable by seed$year.x.x
t = 1.212, df = 293.78, p-value = 0.2265
alternative hypothesis: true difference in means between group 2023 and group 2024 is not equal to 0
95 percent confidence interval:
 -1.386471  5.831912
sample estimates:
mean in group 2023 mean in group 2024 
          35.84713           33.62441

ggplot(seed, aes(ants, total.seeds, fill = treatment)) + geom_smooth(method = "lm")

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 7 rows containing non-finite outside the scale range
(`stat_smooth()`).

Fruit set

# read in fruiting data


fruit <- read.csv("Raw Data/cholla_fruitset.csv")
str(fruit)

'data.frame':   140 obs. of  5 variables:
 $ PlantID: int  1 2 3 4 5 6 7 8 9 10 ...
 $ year   : int  2023 2023 2023 2023 2023 2023 2023 2023 2023 2023 ...
 $ fruit  : int  3 NA 60 2 NA 7 12 17 8 10 ...
 $ partial: int  NA 6 NA 2 1 NA NA NA 5 NA ...
 $ empty  : int  NA 3 NA 3 NA 4 NA NA NA NA ...

fruit[is.na(fruit)] <- 0
fruit <- mutate(fruit, total.fruit = fruit + partial)

fruit$uniID <- paste(fruit$PlantID, fruit$year)


visits <- visits %>% ungroup

fruit <- left_join(visits, fruit, by = "uniID")

fruit <- left_join(fruit, ag, by = "uniID")


mean(fruit$total.fruit)

[1] 19.97857

m1 <- glm(cbind(total.fruit, buds.x) ~ as.factor(crempres) + year.x.x, family = "binomial", data = fruit)
summary(m1)


Call:
glm(formula = cbind(total.fruit, buds.x) ~ as.factor(crempres) + 
    year.x.x, family = "binomial", data = fruit)

Coefficients:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -1.05695    0.04040 -26.163  < 2e-16 ***
as.factor(crempres)1  0.13047    0.04558   2.862 0.004205 ** 
year.x.x2024          0.16140    0.04582   3.523 0.000427 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 953.21  on 136  degrees of freedom
Residual deviance: 932.87  on 134  degrees of freedom
AIC: 1489.4

Number of Fisher Scoring iterations: 4

#super dispersed
#check_overdispersion(m1)

m1 <- glm(cbind(total.fruit, buds.x) ~ as.factor(crempres) + year.x.x, family = "quasibinomial", data = fruit)
summary(m1)

Warning in summary.glm(m1): observations with zero weight not used for
calculating dispersion


Call:
glm(formula = cbind(total.fruit, buds.x) ~ as.factor(crempres) + 
    year.x.x, family = "quasibinomial", data = fruit)

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -1.0569     0.1034 -10.222   <2e-16 ***
as.factor(crempres)1   0.1305     0.1167   1.118    0.265    
year.x.x2024           0.1614     0.1173   1.376    0.171    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasibinomial family taken to be 6.551122)

    Null deviance: 953.21  on 136  degrees of freedom
Residual deviance: 932.87  on 134  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

fruit$crepres <- ifelse(fruit$Crematogaster.depilis > 0, 1, 0)
fruit$crepres <- as.factor(fruit$crepres)



m2 <- glmmTMB(total.fruit ~ as.factor(crempres) +buds.x + (1|year.x.x), family = "poisson", data = fruit)

Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
NA/NaN function evaluation

check_overdispersion(m2)

# Overdispersion test

       dispersion ratio =   11.357
  Pearson's Chi-Squared = 1544.509
                p-value =  < 0.001

Overdispersion detected.

m2.n <- glmmTMB(total.fruit ~ as.factor(crempres) +buds.x + (1|year.x.x), family = "nbinom1", data = fruit)

Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
NA/NaN function evaluation

summary(m2.n)

 Family: nbinom1  ( log )
Formula:          total.fruit ~ as.factor(crempres) + buds.x + (1 | year.x.x)
Data: fruit

     AIC      BIC   logLik deviance df.resid 
  1079.7   1094.4   -534.8   1069.7      135 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 year.x.x (Intercept) 0.02837  0.1684  
Number of obs: 140, groups:  year.x.x, 2

Dispersion parameter for nbinom1 family (): 10.3 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)          2.270976   0.170965  13.283   <2e-16 ***
as.factor(crempres)1 0.162201   0.115741   1.401    0.161    
buds.x               0.010799   0.001277   8.458   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m2.n2 <- glmmTMB(total.fruit ~ as.factor(crempres) +buds.x + (1|year.x.x), family = "nbinom2", data = fruit)
summary(m2.n2)

 Family: nbinom2  ( log )
Formula:          total.fruit ~ as.factor(crempres) + buds.x + (1 | year.x.x)
Data: fruit

     AIC      BIC   logLik deviance df.resid 
  1074.3   1089.0   -532.1   1064.3      135 

Random effects:

Conditional model:
 Groups   Name        Variance  Std.Dev. 
 year.x.x (Intercept) 9.104e-10 3.017e-05
Number of obs: 140, groups:  year.x.x, 2

Dispersion parameter for nbinom2 family (): 1.87 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)           2.02146    0.13169  15.350  < 2e-16 ***
as.factor(crempres)1  0.03315    0.13462   0.246    0.806    
buds.x                0.01634    0.00209   7.820 5.27e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m2.n2)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: total.fruit
                      Chisq Df Pr(>Chisq)    
as.factor(crempres)  0.0606  1     0.8055    
buds.x              61.1576  1  5.269e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(m2, m2.n, m2.n2)

      df      AIC
m2     4 2026.769
m2.n   5 1079.661
m2.n2  5 1074.262

fruit <- mutate(fruit, antpres = ifelse(ants > 0, 1, 0))
m3 <- glmmTMB(total.fruit ~ antpres +buds.x + (1|year.x.x), family = "nbinom2", data = fruit)

summary(m3)

 Family: nbinom2  ( log )
Formula:          total.fruit ~ antpres + buds.x + (1 | year.x.x)
Data: fruit

     AIC      BIC   logLik deviance df.resid 
  1074.0   1088.7   -532.0   1064.0      135 

Random effects:

Conditional model:
 Groups   Name        Variance  Std.Dev. 
 year.x.x (Intercept) 1.447e-09 3.804e-05
Number of obs: 140, groups:  year.x.x, 2

Dispersion parameter for nbinom2 family (): 1.88 

Conditional model:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.083858   0.154947  13.449  < 2e-16 ***
antpres     -0.077137   0.136150  -0.567    0.571    
buds.x       0.016330   0.002082   7.845 4.33e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m3)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: total.fruit
         Chisq Df Pr(>Chisq)    
antpres  0.321  1      0.571    
buds.x  61.542  1  4.334e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ggplot(fruit, aes(as.factor(crepres), total.fruit)) + geom_boxplot() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) + stat_summary(fun.y=mean, geom="point", color = "red", size=4) + ylab("Fruit per plant") + xlab("Crematogaster attendance")