Elevation Results Summary

---

2012 Analyses

For all analyses of the data from 2012, I used GLMM hurdle models, with one model with a binomial distribution and another with a truncated-at-zero distribution. Helpful resources for understanding GLMMs and glmmADMB:

• Bolker et al. 2009

• http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#testing-hypotheses

• http://glmm.wikidot.com/faq

---

1. Does aphid abundance depend on elevation?

Binomial model: Does the presence of aphid colonies in a tree depend on elevation?

• 14.4% (48/333) of trees contained aphid colonies at high elevations.

• 20.6% (49/239) of trees contained aphids colonies at low elevations.

I excluded Tree #310 from this analysis because it contained 37 aphid colonies. Results are qualitatively the same either way.

zero_hurdle_aphids <- glmmadmb(formula = Col_pres ~ Elevation + Time + (1|Elevation:Valley) + (1|Valley:Plot), data = aphids1, family = "binomial") 

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-2.9414325	0.6648900	-4.423939	0.0000097
ElevationLow	0.9498695	0.8749200	1.085664	0.2776275
Time	0.0070657	0.0034825	2.028930	0.0424654

Zero-truncated count model: Does the number of aphid colonies per tree depend on elevation?

I again excluded Tree #310 from this analysis because it contained 37 aphid colonies. Either way there is no significant effect of elevation.

trunc1 <- glmmadmb(formula = Col_num ~ Elevation + Time + (1|Elevation:Valley) + (1|Valley:Plot), data = subaphids, family = "truncnbinom1") # truncated quasi-poisson distribution
trunc2 <- glmmadmb(formula = Col_num ~ Elevation + Time + (1|Elevation:Valley) + (1|Valley:Plot), data = subaphids, family = "truncnbinom") # truncated negative binomial distribution
trunc3 <- glmmadmb(formula = Col_num ~ Elevation + Time + (1|Elevation:Valley) + (1|Valley:Plot), data = subaphids, family = "truncpoiss") # truncated poisson distribution
AICtab(trunc1, trunc2, trunc3)

##        dAIC df
## trunc2  0.0 6 
## trunc1  5.1 6 
## trunc3 67.4 5

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-0.5446851	0.7766400	-0.7013354	0.4830937
ElevationLow	0.2375498	0.3349000	0.7093157	0.4781286
Time	0.0075180	0.0049087	1.5315662	0.1256295

Zero-truncated count model: Does the number of aphids per colony depend on elevation?

I again excluded Tree #310 from this analysis because it contained 37 aphid colonies. Either way the results are qualitatively the same.

# trunc1 <- glmmadmb(formula = Aphids ~ Elevation + Time + (1|Elevation:Valley) + (1|Valley:Plot) + (1|Plot:Tree_num), data = subaphids2, family = "truncnbinom1") # truncated quasi-poisson distribution
trunc2 <- glmmadmb(formula = Aphids ~ Elevation + Time + (1|Elevation:Valley) + (1|Valley:Plot) + (1|Plot:Tree_num), data = subaphids2, family = "truncnbinom") # truncated negative binomial distribution
trunc3 <- glmmadmb(formula = Aphids ~ Elevation + Time + (1|Elevation:Valley) + (1|Valley:Plot) + (1|Plot:Tree_num), data = subaphids2, family = "truncpoiss") # truncated poisson distribution
AICtab(trunc2, trunc3)

##        dAIC   df
## trunc2    0.0 7 
## trunc3 1116.2 6

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	2.4561373	0.1706100	14.396210	0.0000000
ElevationLow	0.4929275	0.1516000	3.251501	0.0011480
Time	-0.0034516	0.0019062	-1.810711	0.0701856

---

2. Does ant tending depend on elevation (colony level measure)?

Binomial model: presence of tending ants

• 83.6% (168/201) of colonies were ant-tended at low elevations.

• 78.2% (104/133) of colonies were ant-tended at high elevations.

hurdle_tended <- glmmadmb(formula = Ant_tended ~ Elevation + Aphids + (1|Elevation:Valley) + (1|Valley:Plot) + (1|Plot:Tree_num), family = "binomial")

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	1.1736601	0.510670	2.298275	0.0215461
ElevationLow	0.9638733	0.687950	1.401080	0.1611900
Aphids	0.0595468	0.018692	3.185682	0.0014441

Zero-truncated count model: number of ants per aphid

trunc1 <- glmmadmb(formula = Ants ~ Elevation*Aphids + (1|Elevation:Valley)+ (1|Valley:Plot) + (1|Plot:Tree_num), data = antsperaph, family = "truncnbinom1")
trunc2 <- glmmadmb(formula = Ants ~ Elevation*Aphids + (1|Elevation:Valley) + (1|Valley:Plot) + (1|Plot:Tree_num), data = antsperaph, family = "truncnbinom")
trunc3 <- glmmadmb(formula = Ants ~ Elevation*Aphids + (1|Elevation:Valley)+ (1|Valley:Plot) + (1|Plot:Tree_num), data = antsperaph, family = "truncpoiss") 
AICtab(trunc1, trunc2, trunc3)

##        dAIC  df
## trunc1   0.0 8 
## trunc3 120.8 7 
## trunc2 134.4 8

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	0.0558959	0.2608000	0.2143249	0.8302937
ElevationLow	0.1265049	0.3178000	0.3980644	0.6905827
Aphids	0.0304148	0.0073610	4.1318897	0.0000360
ElevationLow:Aphids	-0.0007066	0.0083661	-0.0844595	0.9326911

trunc1 <- glmmadmb(formula = Ants ~ Elevation + Aphids + (1|Elevation:Valley)+ (1|Valley:Plot) + (1|Plot:Tree_num), data = antsperaph, family = "truncnbinom1")

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	0.0638676	0.2427600	0.2630894	0.7924817
ElevationLow	0.1146748	0.2853100	0.4019305	0.6877352
Aphids	0.0299299	0.0045788	6.5366250	0.0000000

---

3. Does the abundance of non-tending ants depend on elevation (tree level measure)?

Binomial model:

• 26.4% (69/261) of trees contained non-tending ants at high elevations

• 48.5% (83/171) of trees contained non-tending ants at low elevations

ant_hurdle <- glmmadmb(formula = Ants_pres ~ Elevation + Time + (1|Elevation:Valley) + (1|Valley:Plot), data = antpres, family = "binomial")

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-1.6272803	0.3160600	-5.148644	0.0000003
ElevationLow	0.8124986	0.3876800	2.095797	0.0361002
Time	0.0112308	0.0032807	3.423283	0.0006187

Zero-truncated count model:

trunc1 <- glmmadmb(formula = Ants ~ Elevation + Time + (1|Elevation:Valley) + (1|Valley:Plot), data = antspertree, family = "truncnbinom1")
trunc2 <- glmmadmb(formula = Ants ~ Elevation + Time + (1|Valley:Plot), data = antspertree, family = "truncnbinom") # This model didn't work, so I removed valley.
trunc3 <- glmmadmb(formula = Ants ~ Elevation + Time + (1|Elevation:Valley) + (1|Valley:Plot), data = antspertree, family = "truncpoiss")
AICtab(trunc1, trunc2, trunc3)

##        dAIC df
## trunc1  0.0 6 
## trunc2 40.8 5 
## trunc3 71.3 5

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	0.1874435	0.2108400	0.8890319	0.3739860
ElevationLow	0.1570716	0.2545500	0.6170562	0.5371977
Time	0.0040654	0.0015987	2.5429461	0.0109922

---

2015 Analyses

---

1. Does per aphid colony capita growth rate depend on treatment and elevation (extinct colonies excluded)?

• For testing significance of fixed effects with lmer, prefer F tests over chi-squared tests (http://verde.esalq.usp.br/~jorge/cursos/modelos_longitudinais/Mixed%20Effects%20Models%20in%20S%20and%20S-Plus.pdf, see pg. 92)

• It is inappropriate to use LR tests if model is fit with REML, but you can do F or chi-squared. With F, Kenward-Roger should turn out about the same as chi-squared (http://socserv.socsci.mcmaster.ca/jfox/Courses/soc761/Appendix-Mixed-Models.pdf)

rnoextinctmodel1 <- lmer((Pop_growth+0.106)^2~Elevation*Treatment + (1|Elevation:Valley) + (1|Valley:Clone), data = rnoextinct)

	F	Df	Df.res	Pr(>F)
(Intercept)	43.2037922	1	11.39095	0.0000339
Elevation	0.0088957	1	10.57282	0.9266211
Treatment	3.9680778	2	83.09075	0.0225992
Elevation:Treatment	0.1367008	2	84.46368	0.8724238

rnoextinctmodel2 <- lmer((Pop_growth+0.106)^2~Elevation+Treatment + (1|Elevation:Valley) + (1|Valley:Clone), data = rnoextinct)

	F	Df	Df.res	Pr(>F)
(Intercept)	50.8752158	1	7.089604	0.0001772
Elevation	0.0831587	1	4.052823	0.7872100
Treatment	8.3317765	2	85.224267	0.0004955

summary(glht(rnoextinctmodel2, linfct = mcp(Treatment = "Tukey", interaction_average = T))) # Multiple comparisons between treatments

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = (Pop_growth + 0.106)^2 ~ Elevation + Treatment + 
##     (1 | Elevation:Valley) + (1 | Valley:Clone), data = rnoextinct)
## 
## Linear Hypotheses:
##                     Estimate Std. Error z value Pr(>|z|)    
## Three - Four == 0 -0.0008295  0.0025824  -0.321  0.94418    
## Two - Four == 0    0.0071279  0.0019814   3.597  < 0.001 ***
## Two - Three == 0   0.0079574  0.0025260   3.150  0.00456 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

---

2. Do the effects of predators on colony survival depend on elevation?

• 50% (8/16) of colonies survived at low elevations in the three trophic level treatment

• 67% (16/24) of colonies survived at high elevations in the three trophic level treatment

survival_mod <- glmmadmb(Survival ~ Elevation + (1|Elevation:Valley), family = "binomial", data = pred_effects)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	0.6931481	0.43302	1.60073	0.1094367
ElevationLow	-0.6931496	0.66144	-1.04794	0.2946661

---

3. Do the effects of ants on colony survival depend on elevation?

• 90% (27/30) of the aphid colonies survived in the four trophic level treatment at low elevations

• 73% (24/33) of the aphid colonies survived in the four trophic level treatment at high elevations

survivalmod <- glmmadmb(Survival ~ Elevation*Treatment +(1|Elevation:Valley), family = "binomial", data = ant_effects)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	0.9946933	0.42227	2.3555859	0.0184935
ElevationLow	1.2239975	0.75328	1.6248904	0.1041859
TreatmentThree	-0.2597703	0.59519	-0.4364493	0.6625108
ElevationLow:TreatmentThree	-1.9758773	1.01510	-1.9464854	0.0515965

Do ants increase the proportion of colonies that survived at low elevations?

At low elevations ants increased the proportion of colonies that survived by 40%.

ant_effects_low <- ant_effects %>% filter(Elevation == "Low")
Treatment <- as.factor(ant_effects_low$Corrected_treatment)
Survival <- as.factor(ant_effects_low$Survived)
Valley <- as.factor(ant_effects_low$Valley)
survivalmod1 <- glmmadmb(Survival ~ Treatment + (1|Elevation:Valley), family = "binomial", data = ant_effects_low)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	2.197205	0.60858	3.610380	0.0003057
TreatmentThree	-2.197210	0.78763	-2.789648	0.0052765

Do ants increase the proportion of colonies that survived at high elevations?

At high elevations ants increased the proportion of colonies that survived by 7%.

ant_effects_high <- ant_effects %>% filter(Elevation == "High")
Treatment <- as.factor(ant_effects_high$Corrected_treatment)
Survival <- as.factor(ant_effects_high$Survived)
Valley <- as.factor(ant_effects_high$Valley)
survivalmod2 <- glmmadmb(Survival ~ Treatment + (1|Elevation:Valley), family = "binomial", data = ant_effects_high)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	1.0564332	0.5424	1.9477013	0.0514507
TreatmentThree	-0.1989912	0.6078	-0.3273959	0.7433684

---

4. Does ant tending depend on elevation?

Binomial model: Does the probability of being tended depend on elevation?

• 45% of colonies were tended at high elevations.

• 90% of colonies were tended at low elevations.

The aphid variable included in the model is the average number of aphids during the entire time when the colony existed.

tendingmod2 <- glmmadmb(Tended ~ Elevation + Avg_aphids + (1|Elevation:Valley), family = "binomial", data = anttending)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-2.627808	0.934380	-2.812355	0.0049180
ElevationLow	2.190927	0.873590	2.507958	0.0121431
Avg_aphids	0.142169	0.049144	2.892908	0.0038169

Truncated part of the model: Does the per capita tending rate depend on elevation?

Variables are the maximum number of aphids out of all of the observations when ants were present and the number of ants during that observation.

truncmod1 <- glmmadmb(Ants ~ Elevation*Max_aphids + (1|Elevation:Valley), data = anttending, family = "truncnbinom1")
truncmod2 <- glmmadmb(Ants ~ Elevation*Max_aphids + (1|Elevation:Valley), data = anttending, family = "truncnbinom")
truncmod3 <- glmmadmb(Ants ~ Elevation*Max_aphids + (1|Elevation:Valley), data = anttending, family = "truncpoiss")
AICtab(truncmod1, truncmod2, truncmod3)

##           dAIC df
## truncmod1  0.0 6 
## truncmod2  7.5 6 
## truncmod3 17.8 5

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	1.3773616	0.2456800	5.6063236	0.0000000
ElevationLow	-1.0499038	0.3946200	-2.6605438	0.0078015
Max_aphids	0.0034868	0.0042141	0.8274116	0.4080038
ElevationLow:Max_aphids	0.0142462	0.0062882	2.2655417	0.0234795

5. Does the abundance of aphid predators depend on elevation?

This analysis again uses the average number of aphids across all observations when aphids were present.

Binomial model: Does the presence of predators at a colony depend on elevation?

• 38.6% of colonies had predators at high elevations.

• 39.1% of colonies had predators at low elevations.

hurdle <- glmmadmb(Pred_presence ~ Elevation + Avg_aphids + Treatment + Days + (1|Elevation:Valley) + (1|Valley:Clone), family = "binomial", data = preds1)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-2.6437467	1.28670	-2.0546722	0.0399107
ElevationLow	0.1183374	0.46353	0.2552962	0.7984944
Avg_aphids	-0.0285001	0.01783	-1.5984377	0.1099456
TreatmentThree	0.1548561	0.43668	0.3546214	0.7228733
Days	0.4757552	0.21984	2.1640974	0.0304569

Truncated model: Does the non-zero number of predators per colony depend on elevation?

trunc1 <- glmmadmb(Predators ~ Elevation*Avg_aphids + Treatment + Days + (1|Elevation:Valley) + (1|Valley:Clone), family = "truncnbinom", data = preds1)
trunc2 <- glmmadmb(Predators ~ Elevation*Avg_aphids + Treatment + Days + (1|Elevation:Valley) + (1|Valley:Clone), family = "truncnbinom1", data = preds1)
trunc3 <- glmmadmb(Predators ~ Elevation*Avg_aphids + Treatment + Days + (1|Elevation:Valley) + (1|Valley:Clone), family = "truncpoiss", data = preds1)
AICtab(trunc1, trunc2, trunc3)

##        dAIC df
## trunc2  0.0 9 
## trunc3 16.3 8 
## trunc1 18.3 9

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	0.7597229	1.032100	0.7360942	0.4616734
ElevationLow	0.6561158	0.864880	0.7586206	0.4480795
Avg_aphids	-0.0129305	0.020699	-0.6246928	0.5321727
TreatmentThree	0.1926551	0.316970	0.6078022	0.5433187
Days	-0.0456276	0.178130	-0.2561476	0.7978369
ElevationLow:Avg_aphids	-0.0546584	0.050313	-1.0863674	0.2773165

trunc2 <- glmmadmb(Predators ~ Elevation + Avg_aphids + Treatment + Days + (1|Elevation:Valley) + (1|Valley:Clone), family = "truncnbinom1", data = preds1)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	0.8158246	1.018000	0.8013994	0.4229005
ElevationLow	-0.2162406	0.334510	-0.6464399	0.5179945
Avg_aphids	-0.0248502	0.018111	-1.3721046	0.1700309
TreatmentThree	0.1612467	0.312280	0.5163530	0.6056079
Days	-0.0226291	0.164210	-0.1378061	0.8903937

6. Do ants discover aphid colonies more quickly at low elevations?

truncmod1 <- glmmadmb(Disc_days ~ Elevation + (1|Elevation:Valley), data = disc_data, family = "truncpoiss")
# truncmod2 <- glmmadmb(Disc_days ~ Elevation + (1|Elevation:Valley), data = disc_data, family = "truncnbinom") # model didn't converge
truncmod3 <- glmmadmb(Disc_days ~ Elevation + (1|Elevation:Valley), data = disc_data, family = "truncnbinom1")
AICtab(truncmod1, truncmod3)

##           dAIC df
## truncmod3  0.0 4 
## truncmod1 70.1 3

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	2.0936703	0.36383	5.7545289	0.0000000
ElevationLow	-0.4968188	0.50862	-0.9767977	0.3286693

7. Is ant fidelity to aphid colonies higher at low elevations?

This is a measure of whether once ants discovered an aphid colony they were still tending that colony on the subsequent observation.

• At low elevations fidelity was at 92.3% (24/26).

• At high elevations fidelity was at 78.6% (11/14).

fidelity_model <- glmmadmb(Fidelity ~ Elevation + (1|Elevation:Valley), family = "binomial", data = fidelity_data)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	1.299291	0.65156	1.994124	0.0461385
ElevationLow	1.185628	0.98287	1.206292	0.2277049

---

Ant abundance and resource use

---

1. Does ant abundance in pitfall traps depend on elevation (in 2012 and 2015)?

2015:

I used a GLMM, with the total number of ants per site as the response variable. To account for differences in sampling effort between aspen clones, I included the number of trap-days (combined variable accounting for differences in number of days and differences in number of traps) that the pitfall traps were in the field as a covariate.

antrate <- as.numeric(as.character(Ants/Trap_days))
mod <- lmer(log(antrate) ~ Elevation + (1|Elevation:Valley), data = pitfalls)
resids <- resid(mod)
shapiro.test(resids)

## 
##  Shapiro-Wilk normality test
## 
## data:  resids
## W = 0.97961, p-value = 0.3793

Anova(mod, test.statistic = "F", type = 3)

## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
## 
## Response: log(antrate)
##                  F Df Df.res  Pr(>F)  
## (Intercept) 6.9930  1 3.9371 0.05826 .
## Elevation   0.2149  1 3.9984 0.66706  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lsmeans <- lsmeans(mod, "Elevation")
lsmeans

##  Elevation   lsmean        SE   df    lower.CL upper.CL
##  High      1.051795 0.3976431 4.10 -0.04161929 2.145209
##  Low       1.313504 0.4006302 4.23  0.22450992 2.402498
## 
## Degrees-of-freedom method: satterthwaite 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

2012:

poismod1 <- glmmadmb(Ants ~ Elevation + Trap_days + (1|Elevation:Valley), family = "poisson", data = pitfalls2) # get error message stating that convergence failed
poismod2 <- glmmadmb(Ants ~ Elevation + Trap_days + (1|Elevation:Valley), family = "nbinom", data = pitfalls2)
poismod3 <- glmmadmb(Ants ~ Elevation + Trap_days + (1|Elevation:Valley), family = "nbinom1", data = pitfalls2)
AICtab(poismod1, poismod2, poismod3)

##          dAIC  df
## poismod3   0.0 5 
## poismod2   0.6 5 
## poismod1 200.2 4

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	4.3948440	1.086500	4.0449554	0.0000523
ElevationLow	0.1574721	0.582700	0.2702455	0.7869714
Trap_days	0.0177558	0.032498	0.5463665	0.5848140

---

2. Does ant community composition in pitfall traps depend on elevation?

• Bray-Curtis is the standard dissimilarity measure used in ecological research. It is unaffected by null values.

• I log + 1 transformed the ant count data because these analyses tend to be weighted towards more abundant species, which can be misleading if you happen to sample near foraging trails or nests (see Menke et al. 2014 and Gollan et al. 2010).

• Non-metric dimensional scaling (NMDS) is a method of unconstrained ordination that allows for any distance/dissimilarity measure (works with Bray-Curtis). I could have instead done a canonical analysis of principal coordinates (CAP), which is constrained. However, I believe that people often first look at NMDS for a general first step, and then if there are significant differences detected with the permanova, you could do a CAP.

• NMDS Models with a stress < 0.2 is acceptable.

• PERMANOVA is a non-parametric test.

2015:

set.seed(36)
# ants.permanova <- adonis(formula = log(ants2015+1) ~ Elevation, data = ants.env2015, permutations = 999, method = "bray") # strata is used to test hypotheses with nestedness
# ants.permanova 
nest.test.2015 <- nested.npmanova(log(ants2015+1) ~ Elevation + Valley, data = ants.env2015, permutations = 999, method = "bray") 

## Total sum of squares of distance matrix: 16.08441

nest.test.2015 # This is the correct test for the effect of elevation, with valley nested -- only need to use this package if there are multiple fixed effects

## Total sum of squares for non-parametric manova: 12.1694023358501 
## 
## Nested anova for Valley nested within Elevation 
## 
##           Df SumsofSquares      F N.Perm Pr(>F)    
## Elevation  1        0.1859 0.3508    999  1.000    
## Valley     4        2.1198 3.0625    999  0.001 ***
## Residuals 57        9.8636 0.1730                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Run 0 stress 0.1887176 
## Run 1 stress 0.2189534 
## Run 2 stress 0.2118809 
## Run 3 stress 0.1887198 
## ... Procrustes: rmse 0.00247347  max resid 0.01811251 
## Run 4 stress 0.1905581 
## Run 5 stress 0.2143029 
## Run 6 stress 0.1887157 
## ... New best solution
## ... Procrustes: rmse 0.001848861  max resid 0.01353706 
## Run 7 stress 0.1906784 
## Run 8 stress 0.1994604 
## Run 9 stress 0.1887151 
## ... New best solution
## ... Procrustes: rmse 0.0001510044  max resid 0.00107511 
## ... Similar to previous best
## Run 10 stress 0.215179 
## Run 11 stress 0.1905578 
## Run 12 stress 0.2036346 
## Run 13 stress 0.2049579 
## Run 14 stress 0.2092034 
## Run 15 stress 0.2112605 
## Run 16 stress 0.1896213 
## Run 17 stress 0.1887195 
## ... Procrustes: rmse 0.001990668  max resid 0.01457308 
## Run 18 stress 0.2041191 
## Run 19 stress 0.1938114 
## Run 20 stress 0.1896148 
## *** Solution reached

## 
## Call:
## metaMDS(comm = log(ants2015 + 1), distance = "bray", k = 2, autotransform = F,      noshare = F) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     log(ants2015 + 1) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1887151 
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'log(ants2015 + 1)'

2012:

# ants.permanova <- adonis(formula = log(ants2012+1) ~ Elevation, data = ants.env2012, permutations = 999, method ="bray", strata = ants.env2012$Valley)
# ants.permanova # test for the effect of trap-days
nest.test.2012 <- nested.npmanova(log(ants2012+1) ~ Elevation + Valley, data = ants.env2012, permutations = 999, method = "bray")

## Total sum of squares of distance matrix: 4.682083

nest.test.2012 # test for the effect of elevation, with valley nested

## Total sum of squares for non-parametric manova: 3.66036231929012 
## 
## Nested anova for Valley nested within Elevation 
## 
##           Df SumsofSquares      F N.Perm Pr(>F)
## Elevation  1       0.23231 0.8583    999  0.480
## Valley     4       1.08266 1.2694    999  0.253
## Residuals 11       2.34540 0.2132

## Run 0 stress 0.1298555 
## Run 1 stress 0.1353118 
## Run 2 stress 0.1425926 
## Run 3 stress 0.1353119 
## Run 4 stress 0.1298555 
## ... Procrustes: rmse 5.547594e-06  max resid 9.972617e-06 
## ... Similar to previous best
## Run 5 stress 0.1298555 
## ... Procrustes: rmse 2.110934e-06  max resid 5.994344e-06 
## ... Similar to previous best
## Run 6 stress 0.1298555 
## ... Procrustes: rmse 3.030815e-06  max resid 5.982492e-06 
## ... Similar to previous best
## Run 7 stress 0.1425793 
## Run 8 stress 0.1353118 
## Run 9 stress 0.1425793 
## Run 10 stress 0.1298555 
## ... Procrustes: rmse 2.735272e-06  max resid 4.691525e-06 
## ... Similar to previous best
## Run 11 stress 0.1353117 
## Run 12 stress 0.2121485 
## Run 13 stress 0.1425933 
## Run 14 stress 0.1298555 
## ... New best solution
## ... Procrustes: rmse 2.436696e-06  max resid 5.062304e-06 
## ... Similar to previous best
## Run 15 stress 0.1353118 
## Run 16 stress 0.1298555 
## ... New best solution
## ... Procrustes: rmse 1.578282e-06  max resid 2.830822e-06 
## ... Similar to previous best
## Run 17 stress 0.1916351 
## Run 18 stress 0.1353118 
## Run 19 stress 0.1353117 
## Run 20 stress 0.135312 
## *** Solution reached

## 
## Call:
## metaMDS(comm = log(ants2012 + 1), distance = "bray", k = 2, autotransform = F,      noshare = F) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     log(ants2012 + 1) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1298555 
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'log(ants2012 + 1)'

---

3. Does ant consumption of honey from baits depend on elevation (in 2015)?

Mixture model:

Helpful information about mixture models: http://tinyheero.github.io/2015/10/13/mixture-model.html

You can compute means and standard deviations for individual component distributions of a mixture model. If you are assuming a Gaussian/normal distribution, these models are called “Gaussian Mixture Models” (GMM).

setwd("~/Desktop")
baits <- read.xlsx('Final_Analyses.xlsx', 3) %>% filter(Year == "2015")
Consumption <- as.numeric(as.character(baits$Avg_consumption_rate..mg.h))

model <- normalmixEM(Consumption, mean.constr = c(0, NA)) # I set the mean of the first distribution at 0, since we are assuming that for this distribution no sugar was consumed. Does this sound appropriate?

## number of iterations= 37

summary(model)

## summary of normalmixEM object:
##          comp 1    comp 2
## lambda 0.349322  0.650678
## mu     0.000000 41.674530
## sigma  8.077495 30.171761
## loglik at estimate:  -296.5712

plot(model, density = T)

• 42.% (14/33) of the baits were discovered at high elevations

• 77% (23/30) of the baits were discovered at low elevations

binmod <- glmmadmb(formula = Discovered ~ Elevation + (1|Elevation:Valley), family = "binomial", data = baits)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-0.3007371	0.37157	-0.8093687	0.4183031
ElevationLow	1.4997766	0.58135	2.5798170	0.0098853

consumption_model <- lmer(Consumption ~ Elevation + (1|Elevation:Valley), data = baitspos)

	F	Df	Df.res	Pr(>F)
(Intercept)	11.7869774	1	4.757674	0.0200956
Elevation	0.8875598	1	4.021388	0.3992175