Foraging
Swanne Gordon
10/11/2021
Data entry
data_forage <- read.csv("../data_raw/ForagingAssayData_2021_10_24_9.07.csv")data_forage <- data_forage %>%
mutate(Tutor = factor(ifelse(Population == Tutor_pop, 'same', 'diff'))) %>%
filter(Sex != 'dead') %>%
filter(!(ID %in% c('ALF2_1O2O.2.6', 'AHF2_1O3Y.2.2', 'ALF2_1O3O.3.3', 'ALF2_1P2Y.2.4', 'ALF2_1P2Y.2.8', 'ALF2_3G4G.1.4', 'ALF2_1P3G.4.4', 'ALF2_1O3Y.3.7', 'ALF2_1P4G.3.4'))) %>%
mutate(Sex = droplevels(Sex, c('', 'dead'))) %>%
select(-X, -X.1, -X.2, -X.3)
kable(head(data_forage))| Population | ID | Sex | Mom_ID | Litter | Treatment | Pred_Tutor | Pred | Tutor_pop | Birth.Day | Day | Date | Time | Water | Run | Latency | Attempts | Notes | Tutor |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| AH | AHF2_1G4G.3.3 | F | AHF2_1G4G | 3 | pred-/same | pred-/AH | pred- | AH | 2 | 22-Mar | 3:47 PM | pred- | 1 | 16 | 17 | same | ||
| AH | AHF2_1G4G.3.3 | F | AHF2_1G4G | 3 | pred-/same | pred-/AH | pred- | AH | 1 | 21-Mar | 3:06 PM | pred+ | 1 | 73 | 99 | same | ||
| AH | AHF2_1G4G.6.2 | F | AHF2_1G4G | 6 | pred-/same | pred-/AH | pred- | AH | 1 | 29-Sep | 3:17 PM | pred- | 1 | 35 | 166 | first few in middle | same | |
| AH | AHF2_1G4G.6.2 | F | AHF2_1G4G | 6 | pred-/same | pred-/AH | pred- | AH | 2 | 30-Sep | 3:06 PM | pred+ | 1 | 58 | 7 | same | ||
| AH | AHF2_1G4P.2.1 | F | AHF2_1G4P | 2 | pred-/same | pred-/AH | pred- | AH | 1 | 2-Mar | 3:30 PM | pred+ | 1 | 39 | 55 | same | ||
| AH | AHF2_1G4P.2.1 | F | AHF2_1G4P | 2 | pred-/same | pred-/AH | pred- | AH | 2 | 3-Mar | 4:04 PM | pred- | 1 | 310 | 72 | same |
Attempts/Pecks
First we analyze only females.
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectmodel_pecks <- glm.nb(Attempts ~ Population + Tutor_pop + Pred * Water, subset(data_forage, Sex == 'F'))
kable(summary(model_pecks)$coefficients)| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 3.6732911 | 0.1560754 | 23.535368 | 0.0000000 |
| PopulationAL | -0.8938055 | 0.1289212 | -6.932959 | 0.0000000 |
| Tutor_popAL | -0.1728887 | 0.1286695 | -1.343665 | 0.1790568 |
| Predpred+ | 0.2545440 | 0.1801665 | 1.412826 | 0.1577068 |
| Waterpred+ | -0.3624840 | 0.1846607 | -1.962973 | 0.0496493 |
| Predpred+:Waterpred+ | 0.4207713 | 0.2567416 | 1.638890 | 0.1012362 |
(NOTE: We tested interactions up to 3-way, and removed non-significant ones)
There is a significant effect of population, whereby LP fish eat less than HP fish.
Testing water has an effect whereby fish raised in pred- water forage less in pred+ water, but that is not the case for fish raised in pred+ water.
There is no clear effect of tutor
We can now try to plot the results.
plot_forag <- data_forage %>%
filter(Sex == 'F') %>%
group_by(Population) %>%
ggplot(aes(y = Attempts, x = Pred, fill = Water)) +
geom_boxplot() + facet_wrap(~ Population) +
theme_bw()
plot_forag## Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Now we will repeat the same for males
model_pecks <- glm.nb(Attempts ~ Population + (Pred + Tutor_pop + Water)^2, subset(data_forage, Sex == 'M'))
kable(summary(model_pecks)$coefficients)| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 3.1323972 | 0.1933964 | 16.1967749 | 0.0000000 |
| PopulationAL | -0.5180359 | 0.1393579 | -3.7173050 | 0.0002014 |
| Predpred+ | 0.0521592 | 0.2420096 | 0.2155254 | 0.8293577 |
| Tutor_popAL | -0.0121646 | 0.2307346 | -0.0527211 | 0.9579541 |
| Waterpred+ | -0.1153480 | 0.2329695 | -0.4951207 | 0.6205149 |
| Predpred+:Tutor_popAL | 0.0400546 | 0.2804387 | 0.1428283 | 0.8864258 |
| Predpred+:Waterpred+ | -0.3835555 | 0.2804335 | -1.3677236 | 0.1713986 |
| Tutor_popAL:Waterpred+ | 0.4633094 | 0.2776182 | 1.6688728 | 0.0951426 |
- The population result holds in males. And that’s about it
plot_forag <- data_forage %>%
filter(Sex == 'M') %>%
group_by(Population) %>%
ggplot(aes(y = Attempts, x = Tutor_pop, fill = Water)) +
geom_boxplot() + facet_wrap(~ Population) +
theme_bw()
plot_forag
Latency
This is time-to-event data, and we will analyze it using Cox models.
First we need to convert the Latency variable into a Surv object. We will call it time2peck
library(coxme)
data_forage <- data_forage %>%
mutate(event = ifelse (Latency ==600, 0, 1),
time2peck = Surv(Latency, event = event))Now we can run the cox model
model_latency <- coxph(time2peck ~ Population + Pred + Tutor + Water, subset(data_forage, Sex == 'F'))
kable(summary(model_latency)$coefficients)| coef | exp(coef) | se(coef) | z | Pr(>|z|) | |
|---|---|---|---|---|---|
| PopulationAL | -0.5633603 | 0.5692928 | 0.1162317 | -4.8468724 | 0.0000013 |
| Predpred+ | 0.0882000 | 1.0922065 | 0.1151464 | 0.7659815 | 0.4436873 |
| Tutorsame | 0.2513069 | 1.2857046 | 0.1167362 | 2.1527751 | 0.0313364 |
| Waterpred+ | -0.0663850 | 0.9357705 | 0.1146896 | -0.5788230 | 0.5627086 |
Let’s plot the population differences crossed with water.
library(ggfortify)## Warning: package 'ggfortify' was built under R version 3.6.2data_forage <- data_forage %>%
mutate(popwater = factor(paste(Population, Water)))
autoplot(survfit(time2peck ~ popwater, subset(data_forage, Sex == 'F')))
And the males:
model_latency <- coxph(time2peck ~ Population + Pred + Tutor_pop + Water + Tutor_pop:Water + Population:Pred, subset(data_forage, Sex == 'M'))
kable(summary(model_latency)$coefficients)| coef | exp(coef) | se(coef) | z | Pr(>|z|) | |
|---|---|---|---|---|---|
| PopulationAL | -0.1385752 | 0.8705978 | 0.1669836 | -0.8298734 | 0.4066104 |
| Predpred+ | 0.4161952 | 1.5161817 | 0.1854655 | 2.2440569 | 0.0248287 |
| Tutor_popAL | -0.1666733 | 0.8464761 | 0.1764288 | -0.9447059 | 0.3448091 |
| Waterpred+ | -0.3727893 | 0.6888103 | 0.1780037 | -2.0942789 | 0.0362351 |
| Tutor_popAL:Waterpred+ | 0.3228353 | 1.3810378 | 0.2513577 | 1.2843659 | 0.1990139 |
| PopulationAL:Predpred+ | -0.3768484 | 0.6860201 | 0.2550545 | -1.4775210 | 0.1395360 |
FIGURES
Fig1. Comparison of foraging behaviors (attacks and latency) among the two populations of origin.
fig1a <- ggplot(data_forage, aes(Population, Attempts)) +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.07) +
theme_bw()
fig1a## Warning: Removed 2 rows containing non-finite values (stat_summary).
## Warning: Removed 2 rows containing non-finite values (stat_summary).
missing Fig2b: latency
Fig2. Effect of testing water on the feeding attempts (attacks) of fis from different origins, sexes, and treatments.
fig2 <- data_forage %>%
drop_na(Sex) %>%
mutate(Tutor_pop = recode(Tutor_pop, "AH" = "HP", "AL" = "LP"),
Population = recode(Population, "AH" = "HP", "AL" = "LP"),
Sex = recode(Sex, "M" = "Male", "F" = "Female")) %>%
group_by(Pred, Tutor) %>%
ggplot(aes(y = Latency, x = Water, group = Pred_Tutor, col = Pred, linetype = Tutor_pop)) +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun = mean, geom = "line") +
stat_summary(fun.data = "mean_se", geom = "errorbar", width = 0.1) +
facet_grid(Sex ~ Population, drop = T) +
theme_bw() + scale_colour_grey(start = 0, end = 0.7) +
labs(col = "Raising Environment", linetype ="Tutor Origin") +
theme(strip.background = element_blank(),
strip.text.x = element_text(size = 15), strip.text.y = element_text(size = 15),
axis.title.x = element_text(size = 15),
axis.title.y = element_text(size = 15),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12))
fig2## Warning: Removed 2 rows containing non-finite values (stat_summary).
## Warning: Removed 2 rows containing non-finite values (stat_summary).
## Warning: Removed 2 rows containing non-finite values (stat_summary).
Fig3 Environment (Pred) vs. Learning (Tutor) effect.
MOTHER EFFECTS
We here repeat the above analyses but with Mother ID as a random effect, in order to evaluate potential genetic effects on foraging behavior.
library(MASS)
model_pecks <- glmer.nb(Attempts ~ (Population + Tutor_pop + Pred * Water) + (1|Mom_ID), subset(data_forage, Sex == 'F'))
summary(model_pecks)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(0.8186) ( log )
## Formula: Attempts ~ (Population + Tutor_pop + Pred * Water) + (1 | Mom_ID)
## Data: subset(data_forage, Sex == "F")
##
## AIC BIC logLik deviance df.resid
## 2650.6 2680.5 -1317.3 2634.6 302
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.8963 -0.7299 -0.4263 0.3435 5.1113
##
## Random effects:
## Groups Name Variance Std.Dev.
## Mom_ID (Intercept) 0.01072 0.1035
## Number of obs: 310, groups: Mom_ID, 24
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.6530 0.1643 22.234 < 2e-16 ***
## PopulationAL -0.8767 0.1427 -6.143 8.09e-10 ***
## Tutor_popAL -0.1755 0.1307 -1.343 0.1794
## Predpred+ 0.2590 0.1815 1.427 0.1536
## Waterpred+ -0.3567 0.1847 -1.931 0.0535 .
## Predpred+:Waterpred+ 0.3944 0.2628 1.501 0.1334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PpltAL Ttr_AL Prdpr+ Wtrpr+
## PopulatinAL -0.420
## Tutor_popAL -0.395 -0.026
## Predpred+ -0.595 -0.035 0.114
## Waterpred+ -0.582 0.059 0.015 0.494
## Prdprd+:Wt+ 0.473 -0.065 -0.086 -0.695 -0.715Mother ID variance is quire small. The results remain, although weakened.
Let’s redo it for males.
model_pecks <- glmer.nb(Attempts ~ Population + (Pred + Tutor_pop + Water)^2 + (1|Mom_ID), subset(data_forage, Sex == 'M'))
summary(model_pecks)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(0.9692) ( log )
## Formula: Attempts ~ Population + (Pred + Tutor_pop + Water)^2 + (1 | Mom_ID)
## Data: subset(data_forage, Sex == "M")
##
## AIC BIC logLik deviance df.resid
## 2019.6 2055.2 -999.8 1999.6 250
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.9710 -0.7409 -0.4357 0.2597 4.7032
##
## Random effects:
## Groups Name Variance Std.Dev.
## Mom_ID (Intercept) 0.2315 0.4812
## Number of obs: 260, groups: Mom_ID, 25
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.11252 0.24945 12.477 <2e-16 ***
## PopulationAL -0.59416 0.25217 -2.356 0.0185 *
## Predpred+ 0.06039 0.23866 0.253 0.8002
## Tutor_popAL -0.03390 0.23989 -0.141 0.8876
## Waterpred+ -0.28546 0.22617 -1.262 0.2069
## Predpred+:Tutor_popAL -0.05635 0.29466 -0.191 0.8483
## Predpred+:Waterpred+ -0.38355 0.26969 -1.422 0.1550
## Tutor_popAL:Waterpred+ 0.62252 0.27084 2.298 0.0215 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PpltAL Prdpr+ Ttr_AL Wtrpr+ P+:T_A Pr+:W+
## PopulatinAL -0.561
## Predpred+ -0.445 0.003
## Tutor_popAL -0.563 0.084 0.358
## Waterpred+ -0.437 -0.010 0.322 0.355
## Prdpr+:T_AL 0.326 -0.032 -0.584 -0.580 -0.009
## Prdprd+:Wt+ 0.184 0.057 -0.553 0.016 -0.518 -0.034
## Ttr_ppAL:W+ 0.332 -0.063 -0.032 -0.589 -0.611 0.031 0.000The mother effect is larger for males. A new interaction emerges between population and testing water, such that males tutored by AL have more attempts in pred+ water (unlike males tutored by AH)
Let’s look at latency.
model_latency <- coxme(time2peck ~ Population + Water + Pred * Tutor + (1|Mom_ID), subset(data_forage, Sex == 'F'))
kable(summary(model_latency)$coefficients)## Cox mixed-effects model fit by maximum likelihood
## Data: subset(data_forage, Sex == "F")
## events, n = 307, 310 (2 observations deleted due to missingness)
## Iterations= 9 48
## NULL Integrated Fitted
## Log-likelihood -1470.333 -1451.505 -1440.823
##
## Chisq df p AIC BIC
## Integrated loglik 37.66 6.00 1.3114e-06 25.66 3.29
## Penalized loglik 59.02 13.24 9.5770e-08 32.53 -16.83
##
## Model: time2peck ~ Population + Water + Pred * Tutor + (1 | Mom_ID)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## PopulationAL -0.49042266 0.6123675 0.1694803 -2.89 0.0038
## Waterpred+ -0.04039032 0.9604145 0.1157616 -0.35 0.7300
## Predpred+ -0.27100863 0.7626099 0.1710449 -1.58 0.1100
## Tutorsame -0.04848794 0.9526688 0.1714653 -0.28 0.7800
## Predpred+:Tutorsame 0.54084819 1.7174630 0.2388624 2.26 0.0240
##
## Random effects
## Group Variable Std Dev Variance
## Mom_ID Intercept 0.25870315 0.06692732| x | |
|---|---|
| PopulationAL | -0.4904227 |
| Waterpred+ | -0.0403903 |
| Predpred+ | -0.2710086 |
| Tutorsame | -0.0484879 |
| Predpred+:Tutorsame | 0.5408482 |
We have an interaction between predation treatment and tutor, whereby latency is longer for predator-raised fish that are raised also with a high-predation fish.
model_latency <- coxme(time2peck ~ Population + Pred + Tutor_pop + Water + (1|Mom_ID), subset(data_forage, Sex == 'M'))
kable(summary(model_latency)$coefficients)## Cox mixed-effects model fit by maximum likelihood
## Data: subset(data_forage, Sex == "M")
## events, n = 258, 260
## Iterations= 2 10
## NULL Integrated Fitted
## Log-likelihood -1188.784 -1183.806 -1183.714
##
## Chisq df p AIC BIC
## Integrated loglik 9.96 5.00 0.076511 -0.04 -17.81
## Penalized loglik 10.14 4.09 0.040486 1.97 -12.55
##
## Model: time2peck ~ Population + Pred + Tutor_pop + Water + (1 | Mom_ID)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## PopulationAL -0.29816467 0.7421791 0.1263960 -2.36 0.018
## Predpred+ 0.21359411 1.2381200 0.1268287 1.68 0.092
## Tutor_popAL -0.01905357 0.9811268 0.1255251 -0.15 0.880
## Waterpred+ -0.19424451 0.8234565 0.1256642 -1.55 0.120
##
## Random effects
## Group Variable Std Dev Variance
## Mom_ID Intercept 0.0199699652 0.0003987995| x | |
|---|---|
| PopulationAL | -0.2981647 |
| Predpred+ | 0.2135941 |
| Tutor_popAL | -0.0190536 |
| Waterpred+ | -0.1942445 |
Only population has an effect on males, with AL taking less time to eat their first item.
Effect of RUN
First we will look at whether the average number of runs between populations, predation regime, tutor, and water
model_run <- glmer(Run ~ Population + Pred + Tutor_pop + Water + (1| Mom_ID), family = poisson, subset(data_forage, Sex == 'F'))## boundary (singular) fit: see ?isSingularkable(summary(model_run)$coefficients)| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -0.0571796 | 0.1207121 | -0.4736856 | 0.6357241 |
| PopulationAL | 0.1957950 | 0.1046988 | 1.8700793 | 0.0614728 |
| Predpred+ | 0.0345780 | 0.1049077 | 0.3296044 | 0.7416989 |
| Tutor_popAL | 0.1196548 | 0.1048202 | 1.1415241 | 0.2536519 |
| Waterpred+ | 0.0969107 | 0.1047895 | 0.9248136 | 0.3550629 |