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()
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)
# 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")
| 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")
# 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")
| 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")
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))
# 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")
| 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")
# 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")
| 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")
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
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")
| 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")
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")
| 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")
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))
# 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")
| 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")
# 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")
| 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")
ADD VALUE HERE
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
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")
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")
| 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")
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)")
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.
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