This is an analysis of the ant-aphid fireweed data from 2010. There were several parts to this project: (1) ant feeders placed on Formica podzolica mounds in three high and three low elevation valleys, (2) aphid colonies established on potted fireweed placed next to F. podzolica mounds in these same locations, and (3) ant activity-abundance data from pitfall traps located within the same valleys in 2012 and 2015.

Field site map

Ant mound size and activity

1. Does ant mound size differ with elevation?

sizemodel <- lmer(log(mound.area) ~ elevation + (1|elevation:valley), data = mounds)

hist(resid(sizemodel))

shapiro.test(resid(sizemodel))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(sizemodel)
## W = 0.95682, p-value = 0.03296
qqnorm(resid(sizemodel))
qqline(resid(sizemodel))

## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: log(mound.area)
##                 Chisq Df Pr(>Chisq)    
## (Intercept) 2054.8527  1     <2e-16 ***
## elevation      1.1553  1     0.2824    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##   elevation mound.area
## 1      High   3830.669
## 2       Low   2582.650
##       SE       SE 
## 620.8133 212.1842

Removing the two unusually large colonies

mounds <- read.csv("MoundSites.csv", header = T) %>%
  filter(mound.area < 10000)

sizemodel <- lmer(log(mound.area) ~ elevation + (1|elevation:valley), data = mounds)

hist(resid(sizemodel))

shapiro.test(resid(sizemodel))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(sizemodel)
## W = 0.98012, p-value = 0.4569
qqnorm(resid(sizemodel))
qqline(resid(sizemodel))

## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: log(mound.area)
##                 Chisq Df Pr(>Chisq)    
## (Intercept) 2740.1055  1     <2e-16 ***
## elevation      0.6415  1     0.4232    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##   elevation mound.area
## 1      High   3057.154
## 2       Low   2582.650
##       SE       SE 
## 249.5129 212.1842

2. Does the number of ants on the mound surface depend on elevation?

Number of ants/observation on mound surface:

I excluded two ant mounds with volumes larger than 10,000 cm^2 (both from high elevations). Since ant mounds were observed 7, 8, or 9 times (depending on the valley they were in), I divided the total number of ants by the number of observations.

Area <- Area/10000 # convert to m^2, rather than cm^2
ant.model1 <- lmer(log(Mean_ants) ~ Elevation + (1|Elevation:Valley), data = ant_mounds2)
kable(Anova(ant.model1, type = 3))
Chisq Df Pr(>Chisq)
(Intercept) 258.048140 1 0.0000000
Elevation 3.193988 1 0.0739095 
hist(resid(ant.model1))

shapiro.test(resid(ant.model1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(ant.model1)
## W = 0.99191, p-value = 0.9661
qqnorm(resid(ant.model1))
qqline(resid(ant.model1))

Ant Feeders

1. Did rate of honey consumption (mg/h) depend on honey concentration and elevation?

Dates as repeated measures, setting all negative consumption rates to zero

feeder.model <- lmer(log(Consumption_rate + 1) ~ Elevation*Treatment + Date + (1|Elevation:Valley) + (1|Valley:Mound), data = feeders)

hist(resid(feeder.model))

shapiro.test(resid(feeder.model))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(feeder.model)
## W = 0.9971, p-value = 0.5281
qqnorm(resid(feeder.model))
qqline(resid(feeder.model))

Chisq Df Pr(>Chisq)
(Intercept) 21.884209 1 0.0000029
Elevation 4.852416 1 0.0276074
Treatment 60.572536 2 0.0000000
Date 68.023948 2 0.0000000
Elevation:Treatment 10.748480 2 0.0046344

2. Did the number of F. podzolica per tube depend on elevation and treatment?

Repeated measures

ant.num.model1 <- glmmadmb(AntNum ~ Elevation*Treatment + Date + (1|Elevation:Valley) + (1|Valley:Mound), data = ant_num, family = "gaussian")
# ant.num.model2 <- glmmadmb(AntNum ~ Elevation*Treatment + Date + (1|Elevation:Valley) + (1|Valley:Mound), data = ant_num, family = "poisson")
ant.num.model3 <- glmmadmb(AntNum ~ Elevation*Treatment + Date + (1|Elevation:Valley) + (1|Valley:Mound), data = ant_num, family = "nbinom")
ant.num.model4 <- glmmadmb(AntNum ~ Elevation*Treatment + Date + (1|Elevation:Valley) + (1|Valley:Mound), data = ant_num, family = "nbinom1")

AICtab(ant.num.model1, ant.num.model3, ant.num.model4)
##                dAIC  df
## ant.num.model4   0.0 11
## ant.num.model3   3.9 11
## ant.num.model1 867.8 11
Df Chisq Pr(>Chisq)
(Intercept) 1 24.171627 0.0000009
Elevation 1 2.022687 0.1549648
Treatment 2 29.897678 0.0000003
Date 2 2.166307 0.3385264
Elevation:Treatment 2 1.072164 0.5850358

Interaction removed:

ant.num.model4 <- glmmadmb(AntNum ~ Elevation + Treatment + Date + (1|Elevation:Valley) + (1|Valley:Mound), data = ant_num, family = "nbinom1")
Df Chisq Pr(>Chisq)
(Intercept) 1 35.652645 0.0000000
Elevation 1 2.071661 0.1500581
Treatment 2 60.592059 0.0000000
Date 2 2.138281 0.3433034

3. Did the consumption rate per ant depend on elevation?

To be able to divide consumption rates by the number of ants, we added 1 to each ant # to account for the 0 ant counts.

cons.ant.model <- lmer(log(ConsumptionPerAnt + 1) ~ Elevation*Treatment + Date + (1|Elevation:Valley) + (1|Valley:Mound), data = feeders)

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(cons.ant.model)
## W = 0.99192, p-value = 0.008454

Chisq Df Pr(>Chisq)
(Intercept) 22.863285 1 0.0000017
Elevation 4.817870 1 0.0281661
Treatment 37.010324 2 0.0000000
Date 63.656012 2 0.0000000
Elevation:Treatment 9.315474 2 0.0094879

Figure 1:

Aphid colonies on potted fireweed

1. Do ant tending rates depend on elevation?

  • In the ant tending treatment:
    • 83% (25/30) of all colonies were tended by >1 ant at low elevations.
    • 61% (17/28) of all colonies were tended by >1 ant at high elevations.
  • In the ant exclusion treatment (excluded from subsequent analyses):
    • 40% (12/30) of the colonies were actually tended by >1 ant at low elevations.
    • 29% (8/28) of the colonies were actually tended by >1 at high elevations.

Does the proportion of colonies that were ant tended (in the ant tending treatment) depend on elevation?

hurdlemod <- glmmadmb(Tended ~ elevation + Avg_aphs + (1|elevation:valley), data = tending, family = "binomial")
Df Chisq Pr(>Chisq)
(Intercept) 1 3.211162 0.0731375
elevation 1 5.848910 0.0155867
Avg_aphs 1 5.859612 0.0154922

Does the per capita tending rate (ants/aphid) for colonies that were discovered depend on elevation?

Mean number of ants vs. mean number of aphids across experiment

tendingmod <- lmer(log(Avg_ants) ~ Avg_aphs*elevation + (1|elevation:valley), data = tending)
Anova(tendingmod, type = 3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: log(Avg_ants)
##                     Chisq Df Pr(>Chisq)
## (Intercept)        0.0023  1     0.9618
## Avg_aphs           1.4485  1     0.2288
## elevation          0.1036  1     0.7476
## Avg_aphs:elevation 0.5635  1     0.4528
hist(resid(tendingmod))

shapiro.test(resid(tendingmod))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(tendingmod)
## W = 0.98245, p-value = 0.7555
qqnorm(resid(tendingmod))
qqline(resid(tendingmod))

tendingmod <- lmer(log(Avg_ants) ~ Avg_aphs + elevation + (1|elevation:valley), data = tending)
Anova(tendingmod, type = 3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: log(Avg_ants)
##              Chisq Df Pr(>Chisq)   
## (Intercept) 0.5928  1   0.441350   
## Avg_aphs    8.5079  1   0.003536 **
## elevation   2.6302  1   0.104850   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Negative binomial type I model:

2. Aphid colony survival and per capita growth rates

Survival of aphid colonies

  • Ant exclusion treatment (excluding when >1 ant breached exclusions):
    • 44% (8/18) survived at low elevations
    • 75% (15/20) survived at high elevations
  • Ant presence treatment:
    • 77% (23/30) survived at low elevations
    • 68% (19/28) survived at high elevations
survivmod <- glmmadmb(survived ~ elevation*treatment + (1|elevation:valley) + (1|valley:mound.id), family = "binomial", data = survival)
Df Chisq Pr(>Chisq)
(Intercept) 1 18.719162 0.0000151
elevation 1 3.985765 0.0458863
treatment 1 1.783680 0.1816981
elevation:treatment 1 7.610687 0.0058023

Per capita growth rates of the colonies that survived (excluding when ants snuck into the ant exclusion treatments):

popgrowthmod <- lmer((pop.growth + 0.25)^2 ~ elevation*treatment + (1|elevation:valley) + (1|valley:mound.id), data = survival)

hist(resid(popgrowthmod))

shapiro.test(resid(popgrowthmod))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(popgrowthmod)
## W = 0.98511, p-value = 0.6247
qqnorm(resid(popgrowthmod))
qqline(resid(popgrowthmod))

Chisq Df Pr(>Chisq)
(Intercept) 126.4008111 1 0.0000000
elevation 0.3113624 1 0.5768454
treatment 2.4378340 1 0.1184397
elevation:treatment 0.0283381 1 0.8663164 
popgrowthmod2 <- lmer((pop.growth+0.25)^2 ~ elevation + treatment + (1|elevation:valley) + (1|valley:mound.id), data = survival)
Chisq Df Pr(>Chisq)
(Intercept) 159.1746337 1 0.0000000
elevation 0.6549772 1 0.4183388
treatment 4.8380548 1 0.0278383

Pitfall trap data

2012

Antspertrapday <- ants12$tot.F.podz/ants12$trap.days

pitfall12_model5 <- lmer((Antspertrapday)^(1/3) ~ elevation + (1|elevation:valley), data = ants12)
Anova(pitfall12_model5, type = 3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: (Antspertrapday)^(1/3)
##              Chisq Df Pr(>Chisq)   
## (Intercept) 9.1367  1   0.002505 **
## elevation   0.0058  1   0.939045   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hist(resid(pitfall12_model5))

shapiro.test(resid(pitfall12_model5))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(pitfall12_model5)
## W = 0.90426, p-value = 0.08006
qqnorm(resid(pitfall12_model5))
qqline(resid(pitfall12_model5))

2015

Antspertrapday <- ants15$tot.F.podz/ants15$trap.days

pitfall15_model5 <- lmer((Antspertrapday)^(1/3) ~ elevation + (1|elevation:valley), data = ants15)
Anova(pitfall15_model5, type = 3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: (Antspertrapday)^(1/3)
##               Chisq Df Pr(>Chisq)    
## (Intercept) 15.9872  1  6.377e-05 ***
## elevation    0.2537  1     0.6145    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hist(resid(pitfall15_model5))

shapiro.test(resid(pitfall15_model5))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(pitfall15_model5)
## W = 0.98215, p-value = 0.4924
qqnorm(resid(pitfall15_model5))
qqline(resid(pitfall15_model5))