Time series – aphid colony size and abundance within the plot

time <- demography %>% filter(Location == "Plot" & Surv != "NA" & Ants.t != "Myr") %>%
  group_by(t, Ants.t) %>%
  summarize(z = mean(z), n = length(Ants.t))

ggplot(data = time, aes(t, z)) + 
  geom_point(aes(color = Ants.t, size = n)) + 
  theme_bw()

time <- demography %>% filter(Location == "Plot" & Surv != "NA" & Ants.t != "Myr") %>%
  group_by(t) %>%
  summarize(z = mean(z), n = length(Ants.t))
ggplot(data = time, aes(t, n)) + 
  geom_point() + 
  theme_bw()

Survival

Overall

overall.surv <- demography %>% filter(Surv != "NA" & Ants.t != "Myr" & Ants.t != "Excl")
overall.surv$z <- as.numeric(as.character(overall.surv$z))
overall.surv$Surv <- as.numeric(as.character(overall.surv$Surv))

mod.surv.overall1 <- glm(Surv ~ poly(z, 2) + Ants.t, family = binomial, data = overall.surv)
mod.surv.overall2 <- glm(Surv ~ z + Ants.t, family = binomial, data = overall.surv)
AICtab(mod.surv.overall1, mod.surv.overall2)
##                   dAIC df
## mod.surv.overall1  0.0 6 
## mod.surv.overall2 12.1 5
Anova(mod.surv.overall2)
## Analysis of Deviance Table (Type II tests)
## 
## Response: Surv
##        LR Chisq Df Pr(>Chisq)    
## z       103.986  1    < 2e-16 ***
## Ants.t    8.986  3    0.02948 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glht(mod.surv.overall2, linfct = mcp(Ants.t = "Tukey")))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: glm(formula = Surv ~ z + Ants.t, family = binomial, data = overall.surv)
## 
## Linear Hypotheses:
##                         Estimate Std. Error z value Pr(>|z|)
## F.rufa - F.fusca == 0    0.28663    0.74255   0.386    0.979
## NA - F.fusca == 0       -0.83917    0.41953  -2.000    0.173
## Tapinoma - F.fusca == 0 -0.90905    0.45604  -1.993    0.176
## NA - F.rufa == 0        -1.12580    0.66421  -1.695    0.307
## Tapinoma - F.rufa == 0  -1.19569    0.63174  -1.893    0.215
## Tapinoma - NA == 0      -0.06988    0.30349  -0.230    0.995
## (Adjusted p values reported -- single-step method)

Split by ant tending

# create data files
ants <- c("NA", "Tapinoma", "F.rufa", "F.fusca")
for (i in seq_along(ants)) {
  assign(paste0("surv.", ants[i]), filter(demography, Ants.t == ants[i] & Surv != "NA"))
}

# define the structure of the variables
surv.NA$z <- as.numeric(as.character(surv.NA$z))
surv.NA$Surv <- as.numeric(as.character(surv.NA$Surv))
surv.F.rufa$z <- as.numeric(as.character(surv.F.rufa$z))
surv.F.rufa$Surv <- as.numeric(as.character(surv.F.rufa$Surv))
surv.F.fusca$z <- as.numeric(as.character(surv.F.fusca$z))
surv.F.fusca$Surv <- as.numeric(as.character(surv.F.fusca$Surv))
surv.Tapinoma$z <- as.numeric(as.character(surv.Tapinoma$z))
surv.Tapinoma$Surv <- as.numeric(as.character(surv.Tapinoma$Surv))

# GLMs
mod.surv.NA1 <- glm(Surv ~ poly(z, 2), family = binomial, data = surv.NA)
mod.surv.NA2 <- glm(Surv ~ z, family = binomial, data = surv.NA)
AICtab(mod.surv.NA1, mod.surv.NA2)
##              dAIC df
## mod.surv.NA1 0.0  3 
## mod.surv.NA2 0.7  2
mod.surv.F.rufa1 <- glm(Surv ~ poly(z, 2), family = binomial, data = surv.F.rufa)
mod.surv.F.rufa2 <- glm(Surv ~ z, family = binomial, data = surv.F.rufa)
AICtab(mod.surv.F.rufa1, mod.surv.F.rufa2)
##                  dAIC df
## mod.surv.F.rufa1  0.0 3 
## mod.surv.F.rufa2 13.7 2
mod.surv.F.fusca1 <- glm(Surv ~ poly(z, 2), family = binomial, data = surv.F.fusca)
mod.surv.F.fusca2 <- glm(Surv ~ z, family = binomial, data = surv.F.fusca)
AICtab(mod.surv.F.fusca1, mod.surv.F.fusca2)
##                   dAIC df
## mod.surv.F.fusca1 0.0  3 
## mod.surv.F.fusca2 0.1  2
mod.surv.Tapinoma1 <- glm(Surv ~ poly(z, 2), family = binomial, data = surv.Tapinoma)
mod.surv.Tapinoma2 <- glm(Surv ~ z, family = binomial, data = surv.Tapinoma)
AICtab(mod.surv.Tapinoma1, mod.surv.Tapinoma2)
##                    dAIC df
## mod.surv.Tapinoma1 0.0  3 
## mod.surv.Tapinoma2 0.4  2
surv.mod.df <- data.frame(mod.surv.NA1$coeff, mod.surv.F.rufa1$coeff, mod.surv.F.fusca1$coeff, 
                          mod.surv.Tapinoma1$coeff)
kable(surv.mod.df, col.names = c("No ants", "F. rufa", "F. fusca", "Tapinoma"), caption = "Model coefficients")
Model coefficients
No ants F. rufa F. fusca Tapinoma
(Intercept) 1.004112 1524.969 4.726515 3.868865
poly(z, 2)1 18.486377 15353.301 34.353316 43.161916
poly(z, 2)2 8.124536 6116.464 12.711934 8.869663
# graphs
par(mfrow = c(2,2))

#survNA1 <- transform(surv.NA, z.classes = cut(z, 20))
#survNA1 <- survNA1[order(survNA1$z), ]
#survNAps <- summaryBy(z + Surv ~ z.classes, data = survNA1, na.rm = TRUE)
#survNA1$Surv <- as.numeric(as.character(surv.NA$Surv))
#plot(Surv.mean ~ z.mean, data = survNAps, pch = 19, 
#     xlab = expression("Size t, " * italic(z)), ylab = "Probability of survival")

logi.hist.plot(surv.NA$z, surv.NA$Surv, boxp = F, type = "hist", col = "gray")
mtext(side = 3, line = 0.5, adj = 0, text = "a) No ants")
# lines(fitted(mod.Surv.no) ~ z, data = no.ants.surv, col = "red")

logi.hist.plot(surv.Tapinoma$z, surv.Tapinoma$Surv, boxp = F, type = "hist", col = "gray")
mtext(side = 3, line = 0.5, adj = 0, text = "b) Tapinoma")
# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

logi.hist.plot(surv.F.fusca$z, surv.F.fusca$Surv, boxp = F, type = "hist", col = "gray")
mtext(side = 3, line = 0.5, adj = 0, text = "c) F. fusca")
# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

logi.hist.plot(surv.F.rufa$z, surv.F.rufa$Surv, boxp = F, type = "hist", col = "gray")
mtext(side = 3, line = 0.5, adj = 0, text = "d) F. rufa")

# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

Split by ant tending – stricter ant treatment

# create data files
ants <- c("NA", "Tapinoma", "F.rufa", "F.fusca")
for (i in seq_along(ants)) {
  assign(paste0("surv.", ants[i]), filter(demography, Ants.t.detailed == ants[i] & Surv != "NA"))
}

# define the structure of the variables
surv.NA$z <- as.numeric(as.character(surv.NA$z))
surv.NA$Surv <- as.numeric(as.character(surv.NA$Surv))
surv.F.rufa$z <- as.numeric(as.character(surv.F.rufa$z))
surv.F.rufa$Surv <- as.numeric(as.character(surv.F.rufa$Surv))
surv.F.fusca$z <- as.numeric(as.character(surv.F.fusca$z))
surv.F.fusca$Surv <- as.numeric(as.character(surv.F.fusca$Surv))
surv.Tapinoma$z <- as.numeric(as.character(surv.Tapinoma$z))
surv.Tapinoma$Surv <- as.numeric(as.character(surv.Tapinoma$Surv))

# GLMs
mod.surv.NA <- glm(Surv ~ z, family = binomial, data = surv.NA)
mod.surv.F.rufa <- glm(Surv ~ z, family = binomial, data = surv.F.rufa)
mod.surv.F.fusca <- glm(Surv ~ z, family = binomial, data = surv.F.fusca)
mod.surv.Tapinoma <- glm(Surv ~ z, family = binomial, data = surv.Tapinoma)

surv.mod.df <- data.frame(mod.surv.NA$coeff, mod.surv.F.rufa$coeff, mod.surv.F.fusca$coeff, 
                          mod.surv.Tapinoma$coeff)
kable(surv.mod.df, col.names = c("No ants", "F. rufa", "F. fusca", "Tapinoma"), caption = "Model coefficients")
Model coefficients
No ants F. rufa F. fusca Tapinoma
(Intercept) -0.2378227 -0.4638057 0.6247694 -0.9790755
z 0.6300625 1.2030716 0.6513080 0.9758663
# graphs
par(mfrow = c(2,2))

survNA1 <- transform(surv.NA, z.classes = cut(z, 20))
survNA1 <- survNA1[order(survNA1$z), ]
survNAps <- summaryBy(z + Surv ~ z.classes, data = survNA1, na.rm = TRUE)
survNA1$Surv <- as.numeric(as.character(surv.NA$Surv))
plot(Surv.mean ~ z.mean, data = survNAps, pch = 19, 
     xlab = expression("Size t, " * italic(z)), ylab = "Probability of survival")
mtext(side = 3, line = 0.5, adj = 0, text = "a) No ants")
# lines(fitted(mod.Surv.no) ~ z, data = no.ants.surv, col = "red")

survTap1 <- transform(surv.Tapinoma, z.classes = cut(z, 20))
survTap1 <- survTap1[order(survTap1$z), ]
survTapps <- summaryBy(z + Surv ~ z.classes, data = survTap1, na.rm = TRUE)
survTap1$Surv <- as.numeric(as.character(surv.Tapinoma$Surv))
plot(Surv.mean ~ z.mean, data = survTapps, pch = 19, 
     xlab = expression("Size t, " * italic(z)), ylab = "Probability of survival")
mtext(side = 3, line = 0.5, adj = 0, text = "b) Tapinoma")
# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

survFusca1 <- transform(surv.F.fusca, z.classes = cut(z, 20))
survFusca1 <- survFusca1[order(survFusca1$z), ]
survFuscaps <- summaryBy(z + Surv ~ z.classes, data = survFusca1, na.rm = TRUE)
survFusca1$Surv <- as.numeric(as.character(survFusca1$Surv))
plot(Surv.mean ~ z.mean, data = survFuscaps, pch = 19, 
     xlab = expression("Size t, " * italic(z)), ylab = "Probability of survival")
mtext(side = 3, line = 0.5, adj = 0, text = "c) F. fusca")
# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

survRufa1 <- transform(surv.F.rufa, z.classes = cut(z, 20))
survRufa1 <- survRufa1[order(survRufa1$z), ]
survRufaps <- summaryBy(z + Surv ~ z.classes, data = survRufa1, na.rm = TRUE)
survRufa1$Surv <- as.numeric(as.character(surv.F.rufa$Surv))
plot(Surv.mean ~ z.mean, data = survRufaps, pch = 19, 
     xlab = expression("Size t, " * italic(z)), ylab = "Probability of survival")
mtext(side = 3, line = 0.5, adj = 0, text = "d) F. rufa")

# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

Growth

Overall

overall.growth <- demography %>% filter(Surv == 1 & Ants.t != "Myr" & Ants.t != "Excl")
overall.growth$z <- as.numeric(as.character(overall.growth$z))
overall.growth$Surv <- as.numeric(as.character(overall.growth$Surv))

mod.grow.overall1 <- lm(z1 ~ poly(z,2)*Ants.t, data = overall.growth)
mod.grow.overall2 <- lm(z1 ~ z*Ants.t, data = overall.growth)
Anova(mod.grow.overall2) # BUT residuals aren't normally distributed
## Anova Table (Type II tests)
## 
## Response: z1
##            Sum Sq  Df  F value    Pr(>F)    
## z         1200.02   1 1788.977 < 2.2e-16 ***
## Ants.t      28.64   3   14.231 4.577e-09 ***
## z:Ants.t    23.66   3   11.756 1.472e-07 ***
## Residuals  601.03 896                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#shapiro.test(resid(mod.grow.overall1))
#hist(resid(mod.grow.overall1))
#qqnorm(resid(mod.grow.overall1))
#qqline(resid(mod.grow.overall1))

Split by ant tending

# create data files
ants <- c("NA", "Tapinoma", "F.rufa", "F.fusca")
for (i in seq_along(ants)) {
  assign(paste0("grow.", ants[i]), filter(demography, Ants.t == ants[i] & Surv == 1))
}

# Linear models
mod.grow.NA1 <- lm(z1 ~ poly(z,2), data = grow.NA)
mod.grow.NA2 <- lm(z1 ~ z, data = grow.NA)
AICtab(mod.grow.NA1, mod.grow.NA2)
##              dAIC df
## mod.grow.NA1  0.0 4 
## mod.grow.NA2 19.8 3
mod.grow.F.rufa1 <- lm(z1 ~ poly(z,2), data = grow.F.rufa)
mod.grow.F.rufa2 <- lm(z1 ~ z, data = grow.F.rufa)
AICtab(mod.grow.F.rufa1, mod.grow.F.rufa2)
##                  dAIC df
## mod.grow.F.rufa1 0.0  4 
## mod.grow.F.rufa2 1.5  3
mod.grow.F.fusca1 <- lm(z1 ~ poly(z,2), data = grow.F.fusca)
mod.grow.F.fusca2 <- lm(z1 ~ z, data = grow.F.fusca)
AICtab(mod.grow.F.fusca1, mod.grow.F.fusca2)
##                   dAIC df
## mod.grow.F.fusca1 0.0  4 
## mod.grow.F.fusca2 7.1  3
mod.grow.Tapinoma1 <- lm(z1 ~ poly(z,2), data = grow.Tapinoma)
mod.grow.Tapinoma2 <- lm(z1 ~ z, data = grow.Tapinoma)
AICtab(mod.grow.Tapinoma1, mod.grow.Tapinoma2)
##                    dAIC df
## mod.grow.Tapinoma1  0.0 4 
## mod.grow.Tapinoma2 27.3 3
#shapiro.test(resid(mod.grow.NA))
#hist(resid(mod.grow.NA))
#qqnorm(resid(mod.grow.NA))
#qqline(resid(mod.grow.NA))
#shapiro.test(resid(mod.grow.F.rufa))
#hist(resid(mod.grow.F.rufa))
#qqnorm(resid(mod.grow.F.rufa))
#qqline(resid(mod.grow.F.rufa))
#shapiro.test(resid(mod.grow.F.fusca))
#hist(resid(mod.grow.F.fusca))
#qqnorm(resid(mod.grow.F.fusca))
#qqline(resid(mod.grow.F.fusca))
#shapiro.test(resid(mod.grow.Tapinoma))
#hist(resid(mod.grow.Tapinoma))
#qqnorm(resid(mod.grow.Tapinoma))
#qqline(resid(mod.grow.Tapinoma))

grow.mod.df <- data.frame(mod.grow.NA2$coeff, mod.grow.F.rufa2$coeff, mod.grow.F.fusca2$coeff, 
                          mod.grow.Tapinoma2$coeff)
kable(grow.mod.df, col.names = c("No ants", "F. rufa", "F. fusca", "Tapinoma"), caption = "Model coefficients")
Model coefficients
No ants F. rufa F. fusca Tapinoma
(Intercept) 0.9065345 0.6665212 1.196072 0.4524598
z 0.5869378 0.8804069 0.757278 0.8892087
# graphs
par(mfrow = c(2,2))

plot(z1 ~ z, data = grow.NA, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), ylab = expression("Size t+1, " * italic(z) * "'"))
mtext(side = 3, line = 0.5, adj = 0, text = "a) No ants")
abline(mod.grow.NA2, col = "Red")

plot(z1 ~ z, data = grow.Tapinoma, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), ylab = expression("Size t+1, " * italic(z) * "'"))
mtext(side = 3, line = 0.5, adj = 0, text = "b) Tapinoma")
abline(mod.grow.Tapinoma2, col = "Red")

plot(z1 ~ z, data = grow.F.fusca, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), ylab = expression("Size t+1, " * italic(z) * "'"))
mtext(side = 3, line = 0.5, adj = 0, text = "c) F. fusca")
abline(mod.grow.F.fusca2, col = "Red")

plot(z1 ~ z, data = grow.F.rufa, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), ylab = expression("Size t+1, " * italic(z) * "'"))
mtext(side = 3, line = 0.5, adj = 0, text = "d) F. rufa")
abline(mod.grow.F.rufa2, col = "Red")

Split by ant tending – stricter ant treatment

# create data files
ants <- c("NA", "Tapinoma", "F.rufa", "F.fusca")
for (i in seq_along(ants)) {
  assign(paste0("grow.", ants[i]), filter(demography, Ants.t.detailed == ants[i] & Surv == 1))
}

# Linear models
mod.grow.NA <- lm(z1 ~ z, data = grow.NA)
mod.grow.F.rufa <- lm(z1 ~ z, data = grow.F.rufa)
mod.grow.F.fusca <- lm(z1 ~ z, data = grow.F.fusca)
mod.grow.Tapinoma <- lm(z1 ~ z, data = grow.Tapinoma)

#shapiro.test(resid(mod.grow.NA))
#hist(resid(mod.grow.NA))
#qqnorm(resid(mod.grow.NA))
#qqline(resid(mod.grow.NA))
#shapiro.test(resid(mod.grow.F.rufa))
#hist(resid(mod.grow.F.rufa))
#qqnorm(resid(mod.grow.F.rufa))
#qqline(resid(mod.grow.F.rufa))
#shapiro.test(resid(mod.grow.F.fusca))
#hist(resid(mod.grow.F.fusca))
#qqnorm(resid(mod.grow.F.fusca))
#qqline(resid(mod.grow.F.fusca))
#shapiro.test(resid(mod.grow.Tapinoma))
#hist(resid(mod.grow.Tapinoma))
#qqnorm(resid(mod.grow.Tapinoma))
#qqline(resid(mod.grow.Tapinoma))

grow.mod.df <- data.frame(mod.grow.NA$coeff, mod.grow.F.rufa$coeff, mod.grow.F.fusca$coeff, 
                          mod.grow.Tapinoma$coeff)
kable(grow.mod.df, col.names = c("No ants", "F. rufa", "F. fusca", "Tapinoma"), caption = "Model coefficients")
Model coefficients
No ants F. rufa F. fusca Tapinoma
(Intercept) 0.9065345 0.6373004 1.1931460 0.4787014
z 0.5869378 0.8884791 0.7623144 0.8849949
# graphs
par(mfrow = c(2,2))

plot(z1 ~ z, data = grow.NA, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Size t+1, " * italic(z) * "'"))
mtext(side = 3, line = 0.5, adj = 0, text = "a) No ants")
abline(mod.grow.NA2, col = "Red")

plot(z1 ~ z, data = grow.Tapinoma, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Size t+1, " * italic(z) * "'"))
mtext(side = 3, line = 0.5, adj = 0, text = "b) Tapinoma")
abline(mod.grow.Tapinoma2, col = "Red")

plot(z1 ~ z, data = grow.F.fusca, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Size t+1, " * italic(z) * "'"))
mtext(side = 3, line = 0.5, adj = 0, text = "c) F. fusca")
abline(mod.grow.F.fusca2, col = "Red")

plot(z1 ~ z, data = grow.F.rufa, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Size t+1, " * italic(z) * "'"))
mtext(side = 3, line = 0.5, adj = 0, text = "d) F. rufa")
abline(mod.grow.F.rufa2, col = "Red")


Probability of reproducing

Overall

overall.repr <- demography %>% filter(Repr != "NA" & Ants.t != "Myr" & Ants.t != "Excl")
overall.repr$z <- as.numeric(as.character(overall.repr$z))
overall.repr$Repr <- as.numeric(as.character(overall.repr$Repr))

mod.repr.overall1 <- glm(Repr ~ poly(z,2) + Ants.t, family = binomial, data = overall.repr)
mod.repr.overall2 <- glm(Repr ~ z + Ants.t, family = binomial, data = overall.repr)
AICtab(mod.repr.overall1, mod.repr.overall2)
##                   dAIC df
## mod.repr.overall1  0.0 6 
## mod.repr.overall2 15.6 5
Anova(mod.repr.overall2)
## Analysis of Deviance Table (Type II tests)
## 
## Response: Repr
##        LR Chisq Df Pr(>Chisq)    
## z       101.077  1    < 2e-16 ***
## Ants.t    9.823  3    0.02014 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Split by ant tending

ants <- c("NA", "Tapinoma", "F.rufa", "F.fusca")
for (i in seq_along(ants)) {
  assign(paste0("repr1.", ants[i]), filter(demography, Ants.t == ants[i] & Repr != "NA"))
}

# define the structure of the variables
repr1.NA$z <- as.numeric(as.character(repr1.NA$z))
repr1.NA$Repr <- as.numeric(as.character(repr1.NA$Repr), na.rm = T)
repr1.F.rufa$z <- as.numeric(as.character(repr1.F.rufa$z))
repr1.F.rufa$Repr <- as.numeric(as.character(repr1.F.rufa$Repr))
repr1.F.fusca$z <- as.numeric(as.character(repr1.F.fusca$z))
repr1.F.fusca$Repr <- as.numeric(as.character(repr1.F.fusca$Repr))
repr1.Tapinoma$z <- as.numeric(as.character(repr1.Tapinoma$z))
repr1.Tapinoma$Repr <- as.numeric(as.character(repr1.Tapinoma$Repr))

# GLMs
mod.repr.NA1 <- glm(Repr ~ poly(z,2), family = binomial, data = repr1.NA)
mod.repr.NA2 <- glm(Repr ~ z, family = binomial, data = repr1.NA)
AICtab(mod.repr.NA1, mod.repr.NA2)
##              dAIC df
## mod.repr.NA2 0.0  2 
## mod.repr.NA1 0.6  3
mod.repr.F.rufa1 <- glm(Repr ~ poly(z,2), family = binomial, data = repr1.F.rufa)
mod.repr.F.rufa2 <- glm(Repr ~ z, family = binomial, data = repr1.F.rufa)
AICtab(mod.repr.F.rufa1, mod.repr.F.rufa2)
##                  dAIC df
## mod.repr.F.rufa1 0.0  3 
## mod.repr.F.rufa2 4.6  2
mod.repr.F.fusca1 <- glm(Repr ~ poly(z,2), family = binomial, data = repr1.F.fusca)
mod.repr.F.fusca2 <- glm(Repr ~ z, family = binomial, data = repr1.F.fusca)
AICtab(mod.repr.F.fusca1, mod.repr.F.fusca2)
##                   dAIC df
## mod.repr.F.fusca2 0.0  2 
## mod.repr.F.fusca1 1.1  3
mod.repr.Tapinoma1 <- glm(Repr ~ poly(z,2), family = binomial, data = repr1.Tapinoma)
mod.repr.Tapinoma2 <- glm(Repr ~ z, family = binomial, data = repr1.Tapinoma)
AICtab(mod.repr.Tapinoma1, mod.repr.Tapinoma2)
##                    dAIC df
## mod.repr.Tapinoma1  0.0 3 
## mod.repr.Tapinoma2 10.7 2
repr.mod.df <- data.frame(mod.repr.NA2$coeff, mod.repr.F.rufa2$coeff, mod.repr.F.fusca2$coeff, 
                          mod.repr.Tapinoma2$coeff)
kable(repr.mod.df, col.names = c("No ants", "F. rufa", "F. fusca", "Tapinoma"), caption = "Model coefficients")
Model coefficients
No ants F. rufa F. fusca Tapinoma
(Intercept) -1.1381661 -2.1176986 -0.8405134 -1.5443277
z 0.3721693 0.6592488 0.4417070 0.4636725
# graphs
par(mfrow = c(2,2))

#reprNA1 <- transform(repr1.NA, z.classes = cut(z, 20))
#reprNA1 <- reprNA1[order(reprNA1$z), ]
#reprNAps <- summaryBy(z + Repr ~ z.classes, data = reprNA1, na.rm = TRUE)
#reprNA1$Repr <- as.numeric(as.character(repr1.NA$Repr))
#plot(Repr.mean ~ z.mean, data = reprNAps, pch = 19, 
#     xlab = expression("Size t, " * italic(z)), ylab = "Probability of producing\nalates")
repr1.NA$Surv <- as.numeric(as.character(repr1.NA$Surv))
repr1.Tapinoma$Surv <- as.numeric(as.character(repr1.Tapinoma$Surv))
repr1.F.fusca$Surv <- as.numeric(as.character(repr1.F.fusca$Surv))
repr1.F.rufa$Surv <- as.numeric(as.character(repr1.F.rufa$Surv))

logi.hist.plot(repr1.NA$z, repr1.NA$Surv, boxp = F, type = "hist", col = "gray")
mtext(side = 3, line = 0.5, adj = 0, text = "a) No ants")
# lines(fitted(mod.Surv.no) ~ z, data = no.ants.surv, col = "red")

logi.hist.plot(repr1.Tapinoma$z, repr1.Tapinoma$Surv, boxp = F, type = "hist", col = "gray")
mtext(side = 3, line = 0.5, adj = 0, text = "b) Tapinoma")
# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

logi.hist.plot(repr1.F.fusca$z, repr1.F.fusca$Surv, boxp = F, type = "hist", col = "gray")
mtext(side = 3, line = 0.5, adj = 0, text = "c) F. fusca")
# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

logi.hist.plot(repr1.F.rufa$z, repr1.F.rufa$Surv, boxp = F, type = "hist", col = "gray")
mtext(side = 3, line = 0.5, adj = 0, text = "d) F. rufa")

# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

Split by ant tending – stricter ant treatment

ants <- c("NA", "Tapinoma", "F.rufa", "F.fusca")
for (i in seq_along(ants)) {
  assign(paste0("repr.", ants[i]), filter(demography, Ants.t.detailed == ants[i] & Repr != "NA"))
}

# define the structure of the variables
repr.NA$z <- as.numeric(as.character(repr.NA$z))
repr.NA$Repr <- as.numeric(as.character(repr.NA$Repr))
repr.F.rufa$z <- as.numeric(as.character(repr.F.rufa$z))
repr.F.rufa$Repr <- as.numeric(as.character(repr.F.rufa$Repr))
repr.F.fusca$z <- as.numeric(as.character(repr.F.fusca$z))
repr.F.fusca$Repr <- as.numeric(as.character(repr.F.fusca$Repr))
repr.Tapinoma$z <- as.numeric(as.character(repr.Tapinoma$z))
repr.Tapinoma$Repr <- as.numeric(as.character(repr.Tapinoma$Repr))

# GLMs
mod.repr.NA <- glm(Repr ~ z, family = binomial, data = repr.NA)
mod.repr.F.rufa <- glm(Repr ~ z, family = binomial, data = repr.F.rufa)
mod.repr.F.fusca <- glm(Repr ~ z, family = binomial, data = repr.F.fusca)
mod.repr.Tapinoma <- glm(Repr ~ z, family = binomial, data = repr.Tapinoma)

repr.mod.df <- data.frame(mod.repr.NA$coeff, mod.repr.F.rufa$coeff, mod.repr.F.fusca$coeff, 
                          mod.repr.Tapinoma$coeff)
kable(repr.mod.df, col.names = c("No ants", "F. rufa", "F. fusca", "Tapinoma"), caption = "Model coefficients")
Model coefficients
No ants F. rufa F. fusca Tapinoma
(Intercept) -1.1381661 -2.4253803 -0.9230556 -1.4990878
z 0.3721693 0.7402979 0.4918847 0.4750348
# graphs
par(mfrow = c(2,2))

reprNA1 <- transform(repr.NA, z.classes = cut(z, 20))
reprNA1 <- reprNA1[order(reprNA1$z), ]
reprNAps <- summaryBy(z + Repr ~ z.classes, data = reprNA1, na.rm = TRUE)
reprNA1$Repr <- as.numeric(as.character(repr.NA$Repr))
plot(Repr.mean ~ z.mean, data = reprNAps, pch = 19, 
     xlab = expression("Size t, " * italic(z)), ylab = "Probability of producing\nalates")
mtext(side = 3, line = 0.5, adj = 0, text = "a) No ants")
# lines(fitted(mod.Surv.no) ~ z, data = no.ants.surv, col = "red")

reprTap1 <- transform(repr.Tapinoma, z.classes = cut(z, 20))
reprTap1 <- reprTap1[order(reprTap1$z), ]
reprTapps <- summaryBy(z + Repr ~ z.classes, data = reprTap1, na.rm = TRUE)
reprTap1$Repr <- as.numeric(as.character(repr.Tapinoma$Repr))
plot(Repr.mean ~ z.mean, data = reprTapps, pch = 19, 
     xlab = expression("Size t, " * italic(z)), ylab = "Probability of producing\nalates")
mtext(side = 3, line = 0.5, adj = 0, text = "b) Tapinoma")
# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

reprFusca1 <- transform(repr.F.fusca, z.classes = cut(z, 20))
reprFusca1 <- reprFusca1[order(reprFusca1$z), ]
reprFuscaps <- summaryBy(z + Repr ~ z.classes, data = reprFusca1, na.rm = TRUE)
reprFusca1$Repr <- as.numeric(as.character(repr.F.fusca$Repr))
plot(Repr.mean ~ z.mean, data = reprFuscaps, pch = 19, 
     xlab = expression("Size t, " * italic(z)), ylab = "Probability of producing\nalates")
mtext(side = 3, line = 0.5, adj = 0, text = "c) F. fusca")
# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

reprRufa1 <- transform(repr.F.rufa, z.classes = cut(z, 20))
reprRufa1 <- reprRufa1[order(reprRufa1$z), ]
reprRufaps <- summaryBy(z + Repr ~ z.classes, data = reprRufa1, na.rm = TRUE)
reprRufa1$Repr <- as.numeric(as.character(repr.F.rufa$Repr))
plot(Repr.mean ~ z.mean, data = reprRufaps, pch = 19, 
     xlab = expression("Size t, " * italic(z)), ylab = "Probability of producing\nalates")
mtext(side = 3, line = 0.5, adj = 0, text = "d) F. rufa")

# lines(fitted(mod.Surv.tap) ~ z, data = tap.surv, col = "red")

Alate production

Overall

overall.alates <- demography %>% filter(Repr == 1 & Ants.t != "Myr" & Ants.t != "Excl")
overall.alates$z <- as.numeric(as.character(overall.alates$z))
overall.alates$Alates.t...1 <- as.numeric(as.character(overall.alates$Alates.t...1))

mod.alates.overall1 <- lm(log(Alates.t...1) ~ poly(z,2) + Ants.t, data = overall.alates)
mod.alates.overall2 <- lm(log(Alates.t...1) ~ z*Ants.t, data = overall.alates)
AICtab(mod.alates.overall1, mod.alates.overall2)
##                     dAIC df
## mod.alates.overall1  0.0 7 
## mod.alates.overall2 54.6 9
Anova(mod.alates.overall2)
## Anova Table (Type II tests)
## 
## Response: log(Alates.t...1)
##           Sum Sq  Df  F value    Pr(>F)    
## z         176.27   1 270.5184 < 2.2e-16 ***
## Ants.t      4.67   3   2.3909   0.06762 .  
## z:Ants.t   16.93   3   8.6617 1.218e-05 ***
## Residuals 418.32 642                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#shapiro.test(resid(mod.alates.overall))
#hist(resid(mod.alates.overall))
#qqnorm(resid(mod.alates.overall))
#qqline(resid(mod.alates.overall))

Split by ant tending

# create data files
ants <- c("NA", "Tapinoma", "F.rufa", "F.fusca")
for (i in seq_along(ants)) {
  assign(paste0("alate.", ants[i]), filter(demography, Ants.t == ants[i] & Repr == 1))
}

# Linear models
mod.alate.NA1 <- lm(log(Alates.t...1) ~ poly(z,2), data = alate.NA)
mod.alate.NA2 <- lm(log(Alates.t...1) ~ z, data = alate.NA)
AICtab(mod.alate.NA1, mod.alate.NA2)
##               dAIC df
## mod.alate.NA1 0.0  4 
## mod.alate.NA2 2.9  3
mod.alate.F.rufa1 <- lm(log(Alates.t...1) ~ poly(z,2), data = alate.F.rufa)
mod.alate.F.rufa2 <- lm(log(Alates.t...1) ~ z, data = alate.F.rufa)
AICtab(mod.alate.F.rufa1, mod.alate.F.rufa2)
##                   dAIC df
## mod.alate.F.rufa1  0.0 4 
## mod.alate.F.rufa2 11.8 3
mod.alate.F.fusca1 <- lm(log(Alates.t...1) ~ poly(z,2), data = alate.F.fusca)
mod.alate.F.fusca2 <- lm(log(Alates.t...1) ~ z, data = alate.F.fusca)
AICtab(mod.alate.F.fusca1, mod.alate.F.fusca2)
##                    dAIC df
## mod.alate.F.fusca1 0.0  4 
## mod.alate.F.fusca2 7.2  3
mod.alate.Tapinoma1 <- lm(log(Alates.t...1) ~ poly(z,2), data = alate.Tapinoma)
mod.alate.Tapinoma2 <- lm(log(Alates.t...1) ~ z, data = alate.Tapinoma)
AICtab(mod.alate.Tapinoma1, mod.alate.Tapinoma2)
##                     dAIC df
## mod.alate.Tapinoma1  0.0 4 
## mod.alate.Tapinoma2 30.7 3
#shapiro.test(resid(mod.alate.NA))
#hist(resid(mod.alate.NA))
#qqnorm(resid(mod.alate.NA))
#qqline(resid(mod.alate.NA))
#shapiro.test(resid(mod.alate.F.rufa))
#hist(resid(mod.alate.F.rufa))
#qqnorm(resid(mod.alate.F.rufa))
#qqline(resid(mod.alate.F.rufa))
#shapiro.test(resid(mod.alate.F.fusca))
#hist(resid(mod.alate.F.fusca))
#qqnorm(resid(mod.alate.F.fusca))
#qqline(resid(mod.alate.F.fusca))
#shapiro.test(resid(mod.alate.Tapinoma))
#hist(resid(mod.alate.Tapinoma))
#qqnorm(resid(mod.alate.Tapinoma))
#qqline(resid(mod.alate.Tapinoma))

alate.mod.df <- data.frame(mod.alate.NA2$coeff, mod.alate.F.rufa2$coeff, mod.alate.F.fusca2$coeff, 
                          mod.alate.Tapinoma2$coeff)
kable(alate.mod.df, col.names = c("No ants", "F. rufa", "F. fusca", "Tapinoma"), caption = "Model coefficients")
Model coefficients
No ants F. rufa F. fusca Tapinoma
(Intercept) 0.1917109 -1.2658844 -0.1426351 -0.5841869
z 0.1536506 0.5082736 0.3286193 0.3783253
# graphs
par(mfrow = c(2,2))

plot(log(Alates.t...1) ~ z, data = alate.NA, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Alates t+1"))
mtext(side = 3, line = 0.5, adj = 0, text = "a) No ants")
abline(mod.alate.NA2, col = "Red")

plot(log(Alates.t...1) ~ z, data = alate.Tapinoma, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Alates t+1"))
mtext(side = 3, line = 0.5, adj = 0, text = "b) Tapinoma")
abline(mod.alate.Tapinoma2, col = "Red")

plot(log(Alates.t...1) ~ z, data = alate.F.fusca, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Alates t+1"))
mtext(side = 3, line = 0.5, adj = 0, text = "c) F. fusca")
abline(mod.alate.F.fusca2, col = "Red")

plot(log(Alates.t...1) ~ z, data = alate.F.rufa, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Alates t+1"))
mtext(side = 3, line = 0.5, adj = 0, text = "d) F. rufa")
abline(mod.alate.F.rufa2, col = "Red")

Split by ant tending – stricter ant treatment

# create data files
ants <- c("NA", "Tapinoma", "F.rufa", "F.fusca")
for (i in seq_along(ants)) {
  assign(paste0("alate.", ants[i]), filter(demography, Ants.t.detailed == ants[i] & Repr == 1))
}

# Linear models
mod.alate.NA <- lm(log(Alates.t...1) ~ z, data = alate.NA)
mod.alate.F.rufa <- lm(log(Alates.t...1) ~ z, data = alate.F.rufa)
mod.alate.F.fusca <- lm(log(Alates.t...1) ~ z, data = alate.F.fusca)
mod.alate.Tapinoma <- lm(log(Alates.t...1) ~ z, data = alate.Tapinoma)

#shapiro.test(resid(mod.alate.NA))
#hist(resid(mod.alate.NA))
#qqnorm(resid(mod.alate.NA))
#qqline(resid(mod.alate.NA))
#shapiro.test(resid(mod.alate.F.rufa))
#hist(resid(mod.alate.F.rufa))
#qqnorm(resid(mod.alate.F.rufa))
#qqline(resid(mod.alate.F.rufa))
#shapiro.test(resid(mod.alate.F.fusca))
#hist(resid(mod.alate.F.fusca))
#qqnorm(resid(mod.alate.F.fusca))
#qqline(resid(mod.alate.F.fusca))
#shapiro.test(resid(mod.alate.Tapinoma))
#hist(resid(mod.alate.Tapinoma))
#qqnorm(resid(mod.alate.Tapinoma))
#qqline(resid(mod.alate.Tapinoma))

alate.mod.df <- data.frame(mod.alate.NA$coeff, mod.alate.F.rufa$coeff, mod.alate.F.fusca$coeff, 
                          mod.alate.Tapinoma$coeff)
kable(alate.mod.df, col.names = c("No ants", "F. rufa", "F. fusca", "Tapinoma"), caption = "Model coefficients")
Model coefficients
No ants F. rufa F. fusca Tapinoma
(Intercept) 0.1917109 -1.2499976 -0.1436534 -0.5767551
z 0.1536506 0.5153622 0.3305182 0.3784476
# graphs
par(mfrow = c(2,2))

plot(log(Alates.t...1) ~ z, data = alate.NA, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Alates t+1"))
mtext(side = 3, line = 0.5, adj = 0, text = "a) No ants")
abline(mod.alate.NA, col = "Red")

plot(log(Alates.t...1) ~ z, data = alate.Tapinoma, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Alates t+1"))
mtext(side = 3, line = 0.5, adj = 0, text = "b) Tapinoma")
abline(mod.alate.Tapinoma, col = "Red")

plot(log(Alates.t...1) ~ z, data = alate.F.fusca, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Alates t+1"))
mtext(side = 3, line = 0.5, adj = 0, text = "c) F. fusca")
abline(mod.alate.F.fusca, col = "Red")

plot(log(Alates.t...1) ~ z, data = alate.F.rufa, pch = 16, cex = 0.25,
     xlab = expression("Size t, " * italic(z)), 
     ylab = expression("Alates t+1"))
mtext(side = 3, line = 0.5, adj = 0, text = "d) F. rufa")
abline(mod.alate.F.rufa, col = "Red")

Alate success rate

Value from the literature

ADD VALUE HERE

New colonies in plot

For aphid colonies in the main plot, this is the number of new colonies divided by the number of alates produced in the previous timestep.

success <- demography %>% filter(Location == "Plot")
sum(success$Alates.t) #number of total alates
## [1] 1533
119 #number of colonies founded
## [1] 119
119/2135 #percentage of alates successfully founding colonies
## [1] 0.0557377

Size of new colonies

Maybe the colony with 150 aphids is just an accident and was missed during the previous timestep?? Should I exclude colonies above a certain threshold that we don’t think is possible for a 3 day span?

success <- demography %>% filter(Surv == "NA" & Location == "Plot" & z1 != "-Inf")
hist(exp(success$z1), xlab = "Number of aphids", main = "Histogram")

Ant tending

Per capita tending rate

tending <- demography %>% filter(Ants.t != "Excl" & Ants.t != "Myr" & Ants.t != "NA" & z != "-Inf")
tending$Ant.num.t <- as.numeric(as.character(tending$Ant.num.t))

tend.model <- lm(Ant.num.t ~ Ants.t*z, data = tending)
Anova(tend.model)
## Anova Table (Type II tests)
## 
## Response: Ant.num.t
##           Sum Sq  Df  F value    Pr(>F)    
## Ants.t       638   2   5.2397  0.005487 ** 
## z          17060   1 280.1363 < 2.2e-16 ***
## Ants.t:z    3802   2  31.2154 8.898e-14 ***
## Residuals  48413 795                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tend.ruf.model <- lm(Ant.num.t ~ z, data = subset(tending, Ants.t == "F.rufa"))
tend.tap.model <- lm(Ant.num.t ~ z, data = subset(tending, Ants.t == "Tapinoma"))
tend.fus.model <- lm(Ant.num.t ~ z, data = subset(tending, Ants.t == "F.fusca"))

tend.mod.df <- data.frame(tend.ruf.model$coeff, tend.tap.model$coeff, tend.fus.model$coeff)
kable(tend.mod.df, col.names = c("F. rufa", "F. fusca", "Tapinoma"), caption = "Model coefficients")
Model coefficients
F. rufa F. fusca Tapinoma
(Intercept) -18.448793 -4.935493 -0.750506
z 5.083794 2.959663 1.106511
ggplot(tending, aes(x = z, y = log(as.numeric(as.character(Ant.num.t))), color = Ants.t)) +
  geom_point() +
  stat_smooth(method = "lm", aes(color = Ants.t), se = F) + 
  ylab("Ants")

Colony size

For colonies within the plot, I just picked a random date (3rd observation). Should maybe average things across dates.

size <- demography %>% filter(Surv != "NA" & Location == "Plot" & Ants.t != "Myr")
size$Ants.t <- as.factor(as.character(size$Ants.t))
size$z <- as.numeric(as.character(size$z))

Anova(lm(z ~ Ants.t, data = size)) # should probably be averaged so that I'm not counting the same colonies 7 times
## Anova Table (Type II tests)
## 
## Response: z
##            Sum Sq  Df F value    Pr(>F)    
## Ants.t     912.79   3  155.21 < 2.2e-16 ***
## Residuals 1101.73 562                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glht(lm(z ~ Ants.t, data = size), linfct = mcp(Ants.t = "Tukey")))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = z ~ Ants.t, data = size)
## 
## Linear Hypotheses:
##                         Estimate Std. Error t value Pr(>|t|)    
## F.rufa - F.fusca == 0     2.9214     0.2304  12.682   <0.001 ***
## NA - F.fusca == 0        -0.5497     0.2092  -2.628   0.0425 *  
## Tapinoma - F.fusca == 0   1.1592     0.2108   5.499   <0.001 ***
## NA - F.rufa == 0         -3.4711     0.1675 -20.726   <0.001 ***
## Tapinoma - F.rufa == 0   -1.7622     0.1695 -10.398   <0.001 ***
## Tapinoma - NA == 0        1.7090     0.1394  12.259   <0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
size$Ants.t <- ordered(size$Ants.t, levels = c("NA", "F.fusca", "Tapinoma", "F.rufa"))

ggplot(data = size, aes(Ants.t, z)) +
  geom_boxplot(alpha = 0.8, outlier.colour = "white") + 
  theme_bw() + 
  geom_point(pch = 21, position = position_dodge(width = 0.75), alpha = 0.8) +
  theme(text = element_text(size = 15), panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank()) +
  scale_x_discrete("Ant tending") + 
  scale_y_continuous("Aphid colony size (ln-transformed)")

Transitions

10% (48/502) of colonies transitioned away from being tended by F. rufa. 16% (30/192) of colonies transitioned away from being tended by F. fusca. 9% (45/489) of colonies transitioned away from being tended by Tapinoma.

43% of all transitions were to F. fusca, 30% were to F. rufa, and 24% were to Tapinoma.

Abandonment

6% (29/502) of colonies tended by F. rufa were abandoned. 28% (53/192) of colonies tended by F. fusca were abandoned. 19% (92/489) of colonies tended by Tapinoma were abandoned.

abandoned <- demography %>% filter(Ants.t != "NA" & Ants.t != "Excl" & Ants.t != "Myr" & Ants.t...1 != "Excl")

aggregate(Plant.ID ~ Ants.t, data = abandoned, FUN = "length")
##     Ants.t Plant.ID
## 1  F.fusca      119
## 2   F.rufa      351
## 3 Tapinoma      351
aggregate(Plant.ID ~ Ants.t + Ants.t...1, data = abandoned, FUN = "length")
##      Ants.t Ants.t...1 Plant.ID
## 1   F.fusca    F.fusca       55
## 2    F.rufa    F.fusca       20
## 3  Tapinoma    F.fusca       13
## 4   F.fusca     F.rufa       10
## 5    F.rufa     F.rufa      304
## 6  Tapinoma     F.rufa       16
## 7   F.fusca        Myr        2
## 8  Tapinoma        Myr        1
## 9   F.fusca         NA       42
## 10   F.rufa         NA       17
## 11 Tapinoma         NA       70
## 12  F.fusca   Tapinoma       10
## 13   F.rufa   Tapinoma       10
## 14 Tapinoma   Tapinoma      251