Foraging (Run1)
Andrés López-Sepulcre
10/25/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 |
In this vignette we want to reanalyze the data excluding any data from runs 2 and 3, and instead use data from every individual’s first run. What that means is that we need to substitute the Attempts data for individuals with more than 1 run to be 0, and Latency to be 600+
data_forage <- data_forage %>%
mutate(Attempts = ifelse(Run == 1, 0, Attempts),
Latency = ifelse(Run == 1, 600, Latency),
event = ifelse(Latency == 600, 0, 1),
time2peck = Surv(Latency, event = event))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) | -1.7630478 | 0.6250741 | -2.8205422 | 0.0047943 |
| PopulationAL | 1.4618845 | 0.5482685 | 2.6663662 | 0.0076676 |
| Tutor_popAL | 0.6822655 | 0.5488186 | 1.2431529 | 0.2138115 |
| Predpred+ | 1.0229243 | 0.5493438 | 1.8620839 | 0.0625913 |
| Waterpred+ | 0.3822435 | 0.5475571 | 0.6980889 | 0.4851216 |
(NOTE: We tested interactions up to 3-way, and removed non-significant ones)
There is only a significant effect of population, whereby LP fish eat less than HP fish. The interactions shown with the previous dataset are not there.
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 3 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) | -1.1374418 | 0.8286275 | -1.3726816 | 0.1698513 |
| PopulationAL | 0.9147975 | 0.5906286 | 1.5488541 | 0.1214168 |
| Predpred+ | -0.6440589 | 1.0440982 | -0.6168566 | 0.5373293 |
| Tutor_popAL | 2.2697796 | 0.9717505 | 2.3357637 | 0.0195036 |
| Waterpred+ | 1.4830185 | 0.9768728 | 1.5181285 | 0.1289820 |
| Predpred+:Tutor_popAL | -0.1919244 | 1.1901679 | -0.1612583 | 0.8718900 |
| Predpred+:Waterpred+ | -0.9839417 | 1.1901761 | -0.8267195 | 0.4083961 |
| Tutor_popAL:Waterpred+ | -2.4148880 | 1.1729384 | -2.0588362 | 0.0395099 |
The population result holds in males.
Now there is an interaction between Tutor population and testing water. Individuals with AH tutors react by having more attempts in predator water, but individuals raised with AL tutors don’t.
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## Warning: Removed 2 rows containing non-finite values (stat_boxplot).
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 | 1.3839659 | 3.990697 | 0.3659986 | 3.7813416 | 0.0001560 |
| Predpred+ | -0.1604963 | 0.851721 | 0.3167582 | -0.5066839 | 0.6123766 |
| Tutorsame | 0.5121996 | 1.668958 | 0.3234528 | 1.5835373 | 0.1132991 |
| Waterpred+ | 0.4719807 | 1.603166 | 0.3230156 | 1.4611700 | 0.1439688 |
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.6156348 | 1.8508311 | 0.4585247 | 1.3426427 | 0.1793877 |
| Predpred+ | -1.7836743 | 0.1680197 | 1.0691441 | -1.6683198 | 0.0952523 |
| Tutor_popAL | 0.5015750 | 1.6513201 | 0.5275220 | 0.9508134 | 0.3416991 |
| Waterpred+ | 0.7738679 | 2.1681361 | 0.5000699 | 1.5475192 | 0.1217381 |
| Tutor_popAL:Waterpred+ | -0.7537368 | 0.4706047 | 0.6873517 | -1.0965810 | 0.2728246 |
| PopulationAL:Predpred+ | 2.0012738 | 7.3984744 | 1.1343358 | 1.7642693 | 0.0776867 |
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 5 rows containing non-finite values (stat_summary).
## Warning: Removed 5 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 5 rows containing non-finite values (stat_summary).
## Warning: Removed 5 rows containing non-finite values (stat_summary).
## Warning: Removed 5 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'))## boundary (singular) fit: see ?isSingularsummary(model_pecks)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(0.0471) ( log )
## Formula: Attempts ~ (Population + Tutor_pop + Pred * Water) + (1 | Mom_ID)
## Data: subset(data_forage, Sex == "F")
##
## AIC BIC logLik deviance df.resid
## 515.1 545.0 -249.6 499.1 301
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.2161 -0.2132 -0.2087 -0.2000 10.1201
##
## Random effects:
## Groups Name Variance Std.Dev.
## Mom_ID (Intercept) 1.734e-12 1.317e-06
## Number of obs: 309, groups: Mom_ID, 24
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1381 0.7842 -2.726 0.00641 **
## PopulationAL 1.3795 0.5768 2.392 0.01677 *
## Tutor_popAL 0.8128 0.5948 1.367 0.17178
## Predpred+ 1.5931 0.7869 2.024 0.04292 *
## Waterpred+ 1.0213 0.8512 1.200 0.23018
## Predpred+:Waterpred+ -1.1648 1.1314 -1.029 0.30326
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PpltAL Ttr_AL Prdpr+ Wtrpr+
## PopulatinAL -0.370
## Tutor_popAL -0.539 -0.132
## Predpred+ -0.652 0.061 0.192
## Waterpred+ -0.679 -0.039 0.379 0.543
## Prdprd+:Wt+ 0.430 0.155 -0.244 -0.700 -0.735
## convergence code: 0
## boundary (singular) fit: see ?isSingularMother ID variance is very small. There is an effect of population (AL has longer latency), and an effect of Predator treatment while raising (higher latency than control)
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'))## boundary (singular) fit: see ?isSingularsummary(model_pecks)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(0.0485) ( log )
## Formula: Attempts ~ Population + (Pred + Tutor_pop + Water)^2 + (1 | Mom_ID)
## Data: subset(data_forage, Sex == "M")
##
## AIC BIC logLik deviance df.resid
## 466.0 501.5 -223.0 446.0 248
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.2195 -0.2184 -0.2138 -0.2031 5.5540
##
## Random effects:
## Groups Name Variance Std.Dev.
## Mom_ID (Intercept) 1.071e-12 1.035e-06
## Number of obs: 258, groups: Mom_ID, 25
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1374 0.7781 -1.462 0.1438
## PopulationAL 0.9148 0.8130 1.125 0.2605
## Predpred+ -0.6440 1.1765 -0.547 0.5841
## Tutor_popAL 2.2698 0.9833 2.308 0.0210 *
## Waterpred+ 1.4830 0.9804 1.513 0.1304
## Predpred+:Tutor_popAL -0.1919 1.2559 -0.153 0.8785
## Predpred+:Waterpred+ -0.9839 1.1818 -0.833 0.4051
## Tutor_popAL:Waterpred+ -2.4149 1.3161 -1.835 0.0665 .
## ---
## 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.284
## Predpred+ -0.372 -0.470
## Tutor_popAL -0.683 0.144 0.201
## Waterpred+ -0.643 0.011 0.259 0.375
## Prdpr+:T_AL 0.224 0.339 -0.652 -0.439 -0.003
## Prdprd+:Wt+ 0.332 0.025 -0.534 -0.013 -0.529 0.043
## Ttr_ppAL:W+ 0.483 -0.461 0.212 -0.604 -0.555 -0.144 0.016
## convergence code: 0
## boundary (singular) fit: see ?isSingularThe mother effect extremely small. The interaction between population and testing water has weakened, while there is now a clear effect of Tutor (fish with AL tutors take longer to attempt)
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 = 40, 309 (3 observations deleted due to missingness)
## Iterations= 2 12
## NULL Integrated Fitted
## Log-likelihood -226.6944 -216.3386 -216.3245
##
## Chisq df p AIC BIC
## Integrated loglik 20.71 6.00 0.00206680 8.71 -1.42
## Penalized loglik 20.74 5.01 0.00091796 10.71 2.24
##
## Model: time2peck ~ Population + Water + Pred * Tutor + (1 | Mom_ID)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## PopulationAL 1.38099187 3.9788462 0.3667215 3.77 0.00017
## Waterpred+ 0.46964268 1.5994226 0.3234093 1.45 0.15000
## Predpred+ -0.10358354 0.9016007 0.5010620 -0.21 0.84000
## Tutorsame 0.55764967 1.7465627 0.4494921 1.24 0.21000
## Predpred+:Tutorsame -0.09493848 0.9094289 0.6482168 -0.15 0.88000
##
## Random effects
## Group Variable Std Dev Variance
## Mom_ID Intercept 0.0199924335 0.0003996974| x | |
|---|---|
| PopulationAL | 1.3809919 |
| Waterpred+ | 0.4696427 |
| Predpred+ | -0.1035835 |
| Tutorsame | 0.5576497 |
| Predpred+:Tutorsame | -0.0949385 |
By using Run 1 only, we have lost any of the effects we saw before in this analysis.
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 = 36, 258 (2 observations deleted due to missingness)
## Iterations= 2 12
## NULL Integrated Fitted
## Log-likelihood -197.3443 -191.3689 -191.3581
##
## Chisq df p AIC BIC
## Integrated loglik 11.95 5.00 0.035469 1.95 -5.97
## Penalized loglik 11.97 4.01 0.017728 3.95 -2.41
##
## Model: time2peck ~ Population + Pred + Tutor_pop + Water + (1 | Mom_ID)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## PopulationAL 1.18946940 3.2853375 0.4020629 2.96 0.0031
## Predpred+ -0.14166785 0.8679095 0.3422821 -0.41 0.6800
## Tutor_popAL 0.06671849 1.0689945 0.3339240 0.20 0.8400
## Waterpred+ 0.39299454 1.4814103 0.3382896 1.16 0.2500
##
## Random effects
## Group Variable Std Dev Variance
## Mom_ID Intercept 0.0199745045 0.0003989808| x | |
|---|---|
| PopulationAL | 1.1894694 |
| Predpred+ | -0.1416678 |
| Tutor_popAL | 0.0667185 |
| Waterpred+ | 0.3929945 |
Only population has an effect on males, with AL taking less time to eat their first item.