Data Analysis

Load data: Individuals that didn’t eat in first trial, and thus didn’t have a second trial were removed. Individuals that did not eat during the first run (but did eat on the 2nd or 3rd run) had their attempts changed to 0 and latency changed to 600+ to more appropriately compare to fish that ate in Run 1. So all data represents behaviors in the 1st run of the fish for each water treatment.

load(file = here::here('data/data_forage.Rda'))

library(knitr)

kable(head(data_forage))
## Warning in `[<-.data.frame`(`*tmp*`, , isn, value = structure(list(Litter =
## structure(c("3", : provided 9 variables to replace 8 variables
Population ID Sex Mom_ID Litter Treatment Pred_Tutor Pred Tutor_pop Birth.Day Day Date Time Water Weight Run Latency Attempts Notes Tutor event time2peck
AH AHF2_1G4G.3.3 F AHF2_1G4G 3 pred-/same pred-/AH pred- AH 2-Feb 2 22-Mar 15:47 pred- 0.133 1 16 17 same 1 16
AH AHF2_1G4G.3.3 F AHF2_1G4G 3 pred-/same pred-/AH pred- AH 2-Feb 1 21-Mar 15:06 pred+ 0.133 1 73 99 same 1 73
AH AHF2_1G4G.6.2 F AHF2_1G4G 6 pred-/same pred-/AH pred- AH 14-Aug 1 29-Sep 15:17 pred- 0.144 1 35 166 first few in middle same 1 35
AH AHF2_1G4G.6.2 F AHF2_1G4G 6 pred-/same pred-/AH pred- AH 14-Aug 2 30-Sep 15:06 pred+ 0.144 1 58 7 same 1 58
AH AHF2_1G4P.2.1 F AHF2_1G4P 2 pred-/same pred-/AH pred- AH 15-Jan 1 2-Mar 15:30 pred+ 0.128 1 39 55 same 1 39
AH AHF2_1G4P.2.1 F AHF2_1G4P 2 pred-/same pred-/AH pred- AH 15-Jan 2 3-Mar 16:04 pred- 0.128 1 310 72 same 1 310

Attempts/Pecks

First we analyze only females

library(MASS)
library(knitr)
library(lme4)

# 'Simplest' model looking at all four variables independently
model_pecks_simple <- glm.nb(Attempts ~ Population + Tutor_pop + Pred + Water, subset(data_forage, Sex == 'F'))

# 'Full' model looking at all four variables independently and with all possible two way interactions
model_pecks_full <- glm.nb(Attempts ~ (Population + Tutor_pop + Pred + Water)^2, subset(data_forage, Sex == 'F'))

# 'Trimmed' model looking at all the signficant or trending on significant interactions in the 'Full' model

model_pecks_trim <- glm.nb(Attempts ~ Population + Tutor_pop + Pred * Water, subset(data_forage, Sex == 'F'))


kable(summary(model_pecks_simple)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.5987752 0.1744225 20.632515 0.0000000
PopulationAL -1.0215722 0.1586079 -6.440864 0.0000000
Tutor_popAL -0.2311055 0.1582042 -1.460805 0.1440691
Predpred+ 0.4573438 0.1581425 2.891972 0.0038283
Waterpred+ -0.1732993 0.1578257 -1.098042 0.2721862
kable(summary(model_pecks_full)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.8661537 0.2545692 15.1870442 0.0000000
PopulationAL -1.1241208 0.3173216 -3.5425279 0.0003963
Tutor_popAL -0.3916102 0.3070727 -1.2753014 0.2022026
Predpred+ 0.0066733 0.3033832 0.0219964 0.9824508
Waterpred+ -0.5599422 0.3085409 -1.8148066 0.0695536
PopulationAL:Tutor_popAL -0.0250555 0.3154572 -0.0794260 0.9366938
PopulationAL:Predpred+ 0.4069069 0.3154219 1.2900401 0.1970367
PopulationAL:Waterpred+ -0.2122775 0.3147673 -0.6743951 0.5000602
Tutor_popAL:Predpred+ -0.0733019 0.3140665 -0.2333960 0.8154539
Tutor_popAL:Waterpred+ 0.3497376 0.3138509 1.1143432 0.2651320
Predpred+:Waterpred+ 0.5918235 0.3137729 1.8861523 0.0592744
kable(summary(model_pecks_trim)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.7460114 0.1915669 19.5545886 0.0000000
PopulationAL -1.0384631 0.1580132 -6.5720009 0.0000000
Tutor_popAL -0.2555276 0.1575930 -1.6214397 0.1049234
Predpred+ 0.1893597 0.2218712 0.8534665 0.3934006
Waterpred+ -0.4524229 0.2264190 -1.9981663 0.0456986
Predpred+:Waterpred+ 0.5412919 0.3146490 1.7203042 0.0853772
# Need to take into account random effects 
model_pecks_simple_id <- glmer.nb(Attempts ~ Population + Tutor_pop + Pred + Water + (1|ID), subset(data_forage, Sex == 'F'))

# model_pecks_full <- glmer.nb(Attempts ~ (Population + Tutor_pop + Pred + Water)^2 + (1|ID), subset(data_forage, Sex == 'F'))

model_pecks_trim_id <- glmer.nb(Attempts ~ Population + Tutor_pop + Pred * Water + (1|ID), subset(data_forage, Sex == 'F'))

kable(summary(model_pecks_simple_id)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.5523136 0.1816699 19.553670 0.0000000
PopulationAL -1.0781444 0.1778470 -6.062201 0.0000000
Tutor_popAL -0.2252319 0.1678341 -1.341991 0.1795988
Predpred+ 0.4404657 0.1700072 2.590865 0.0095735
Waterpred+ -0.1888288 0.1616081 -1.168436 0.2426308
# kable(summary(model_pecks_full)$coefficients)
kable(summary(model_pecks_trim_id)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.6908451 0.2059306 17.9227588 0.0000000
PopulationAL -1.1078580 0.1793936 -6.1755715 0.0000000
Tutor_popAL -0.2456015 0.1698581 -1.4459218 0.1481991
Predpred+ 0.1539084 0.2330627 0.6603736 0.5090141
Waterpred+ -0.4846872 0.2287561 -2.1187948 0.0341078
Predpred+:Waterpred+ 0.5719984 0.3190408 1.7928692 0.0729938
# the results with random effects included are the same as when they were not included

(NOTE: We tested interactions up to 2-way, and removed non-signficant interactions) There are only significant effects of population, whereby LP fish eat less than HP fish, and rearing water, where fish reared in pred + water ate more than fish reared in pred- water.

Now we will do the same for males

# 'Simplest' model looking at all four variables independently
model_pecks_simple <- glm.nb(Attempts ~ Population + Tutor_pop + Pred + Water, subset(data_forage, Sex == 'M'))

# 'Full' model looking at all four variables independently and with all possible two way interactions
model_pecks_full <- glm.nb(Attempts ~ (Population + Tutor_pop + Pred + Water)^2, subset(data_forage, Sex == 'M'))

# 'Trimmed' model looking at all the signficant or trending on significant interactions in the 'Full' model

model_pecks_trim <- glm.nb(Attempts ~ Population + Pred + Water * Tutor_pop, subset(data_forage, Sex == 'M'))


kable(summary(model_pecks_simple)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.0155154 0.1974789 15.2700635 0.0000000
PopulationAL -0.5316341 0.1748043 -3.0413101 0.0023555
Tutor_popAL 0.1555127 0.1744521 0.8914350 0.3726958
Predpred+ -0.0438475 0.1756602 -0.2496152 0.8028849
Waterpred+ -0.0391066 0.1739754 -0.2247824 0.8221485
kable(summary(model_pecks_full)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.8530170 0.2927564 9.7453624 0.0000000
PopulationAL -0.1461726 0.3378470 -0.4326590 0.6652625
Tutor_popAL -0.0544962 0.3436567 -0.1585773 0.8740019
Predpred+ 0.3749826 0.3574287 1.0491116 0.2941268
Waterpred+ -0.0460824 0.3506592 -0.1314165 0.8954459
PopulationAL:Tutor_popAL -0.1675185 0.3457170 -0.4845537 0.6279930
PopulationAL:Predpred+ -0.5016977 0.3486900 -1.4388071 0.1502052
PopulationAL:Waterpred+ -0.2122664 0.3455806 -0.6142314 0.5390624
Tutor_popAL:Predpred+ 0.0714742 0.3482628 0.2052306 0.8373920
Tutor_popAL:Waterpred+ 0.5841628 0.3448415 1.6940035 0.0902646
Predpred+:Waterpred+ -0.4303792 0.3473643 -1.2389850 0.2153511
kable(summary(model_pecks_trim)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.1764835 0.2140655 14.8388379 0.0000000
PopulationAL -0.5447368 0.1738456 -3.1334522 0.0017276
Predpred+ -0.0599176 0.1747094 -0.3429557 0.7316318
Waterpred+ -0.3566677 0.2441707 -1.4607312 0.1440892
Tutor_popAL -0.1629372 0.2449072 -0.6653016 0.5058576
Waterpred+:Tutor_popAL 0.6341565 0.3460624 1.8324921 0.0668781
# random effects
model_pecks_simple_id <- glmer.nb(Attempts ~ Population + Tutor_pop + Pred + Water + (1|ID), subset(data_forage, Sex == 'M'))
model_pecks_trim_id <- glmer.nb(Attempts ~ Population + Pred + Water * Tutor_pop + (1|ID), subset(data_forage, Sex == 'M'))

kable(summary(model_pecks_simple_id)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.7007875 0.2281677 11.8368535 0.0000000
PopulationAL -0.6055147 0.2065172 -2.9320308 0.0033675
Tutor_popAL 0.1382724 0.2064848 0.6696491 0.5030815
Predpred+ 0.0556513 0.2071810 0.2686119 0.7882283
Waterpred+ -0.1561811 0.1703630 -0.9167548 0.3592711
kable(summary(model_pecks_trim_id)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.8783926 0.2454925 11.7249726 0.0000000
PopulationAL -0.6417831 0.2078371 -3.0879140 0.0020157
Predpred+ 0.0618702 0.2081572 0.2972282 0.7662923
Waterpred+ -0.5862757 0.2354779 -2.4897275 0.0127841
Tutor_popAL -0.2751185 0.2620306 -1.0499482 0.2937419
Waterpred+:Tutor_popAL 0.8372672 0.3301435 2.5360705 0.0112104
# The trimmed model shows different results compared to the model with out random effects - water, population, and water*turor interactions are signficant

(NOTE: We tested interactions up to 2-way, and removed non-signficant interactions) There is only a significant effect of population, whereby LP fish eat less than HP fish. But a trend with Tutor x Water.

Latency

This is time-to-event data, and we will analyze it using Cox models.

When creating the dataframe, we converted the Latency into a Surve object that we called time2peck

library(coxme)

# 'Simplest' model looking at all four variables independently
model_latency_simple <- coxph(time2peck ~ Population + Pred + Tutor + Water, subset(data_forage, Sex == 'F'))

# 'Full' model looking at all four variables independently and with all possible two way interactions
model_latency_full <- coxph(time2peck ~ (Population + Pred + Tutor + Water)^2, subset(data_forage, Sex == 'F'))

# 'Trimmed' model looking at all the signficant or trending on significant interactions in the 'Full' model
model_latency_trim <- coxph(time2peck ~ Population + Pred + Tutor + Water + Population:Tutor + Pred:Tutor + Population:Water, subset(data_forage, Sex == 'F'))

kable(summary(model_latency_simple)$coefficients)
coef exp(coef) se(coef) z Pr(>|z|)
PopulationAL -0.7567577 0.4691852 0.1266876 -5.9734177 0.0000000
Predpred+ 0.0416949 1.0425764 0.1231710 0.3385127 0.7349769
Tutorsame 0.0559924 1.0575896 0.1235532 0.4531843 0.6504160
Waterpred+ -0.1297774 0.8782909 0.1231942 -1.0534378 0.2921404
kable(summary(model_latency_full)$coefficients)
coef exp(coef) se(coef) z Pr(>|z|)
PopulationAL -0.2847564 0.7521975 0.2452942 -1.1608769 0.2456920
Predpred+ -0.4077995 0.6651122 0.2444727 -1.6680781 0.0953002
Tutorsame 0.2014658 1.2231944 0.2351556 0.8567338 0.3915920
Waterpred+ 0.0229861 1.0232523 0.2387353 0.0962827 0.9232961
PopulationAL:Predpred+ 0.2161835 1.2413302 0.2523266 0.8567606 0.3915772
PopulationAL:Tutorsame -0.7497399 0.4724894 0.2535380 -2.9571105 0.0031054
PopulationAL:Waterpred+ -0.4635441 0.6290503 0.2528427 -1.8333297 0.0667536
Predpred+:Tutorsame 0.5048453 1.6567293 0.2493055 2.0250068 0.0428667
Predpred+:Waterpred+ 0.2265151 1.2542216 0.2476186 0.9147742 0.3603102
Tutorsame:Waterpred+ -0.2222171 0.8007415 0.2493463 -0.8911986 0.3728226
kable(summary(model_latency_trim)$coefficients)
coef exp(coef) se(coef) z Pr(>|z|)
PopulationAL -0.1953424 0.8225529 0.2080349 -0.9389889 0.3477364
Predpred+ -0.2071615 0.8128884 0.1732515 -1.1957268 0.2318032
Tutorsame 0.1030820 1.1085824 0.2055906 0.5013947 0.6160934
Waterpred+ 0.0203545 1.0205630 0.1599622 0.1272455 0.8987461
PopulationAL:Tutorsame -0.7231762 0.4852087 0.2523996 -2.8652030 0.0041674
Predpred+:Tutorsame 0.4798674 1.6158601 0.2475397 1.9385475 0.0525565
PopulationAL:Waterpred+ -0.4380980 0.6452625 0.2510058 -1.7453704 0.0809204
library(coxme)

# 'Simplest' model looking at all four variables independently
model_latency_simple <- coxph(time2peck ~ Population + Pred + Tutor + Water, subset(data_forage, Sex == 'M'))

# 'Full' model looking at all four variables independently and with all possible two way interactions
model_latency_full <- coxph(time2peck ~ (Population + Pred + Tutor + Water)^2, subset(data_forage, Sex == 'M'))

# 'Trimmed' model looking at all the signficant or trending on significant interactions in the 'Full' model
model_latency_trim <- coxph(time2peck ~ Population * Pred + Tutor + Water, subset(data_forage, Sex == 'M'))

kable(summary(model_latency_simple)$coefficients)
coef exp(coef) se(coef) z Pr(>|z|)
PopulationAL -0.5297270 0.5887657 0.1364147 -3.883211 0.0001031
Predpred+ 0.2583134 1.2947445 0.1364153 1.893581 0.0582807
Tutorsame 0.0237455 1.0240297 0.1352951 0.175509 0.8606797
Waterpred+ -0.2227339 0.8003278 0.1354096 -1.644889 0.0999927
kable(summary(model_latency_full)$coefficients)
coef exp(coef) se(coef) z Pr(>|z|)
PopulationAL -0.0713313 0.9311533 0.2567880 -0.2777829 0.7811790
Predpred+ 0.5277458 1.6951069 0.2725409 1.9363909 0.0528198
Tutorsame 0.1791638 1.1962167 0.2681231 0.6682147 0.5039965
Waterpred+ -0.0384037 0.9623243 0.2653448 -0.1447315 0.8849229
PopulationAL:Predpred+ -0.6002284 0.5486863 0.2749177 -2.1833026 0.0290135
PopulationAL:Tutorsame -0.0976370 0.9069781 0.2731098 -0.3575009 0.7207169
PopulationAL:Waterpred+ -0.3299970 0.7189259 0.2728491 -1.2094487 0.2264905
Predpred+:Tutorsame -0.0780707 0.9248991 0.2731662 -0.2857992 0.7750319
Predpred+:Waterpred+ 0.1196902 1.1271476 0.2741843 0.4365318 0.6624510
Tutorsame:Waterpred+ -0.2046928 0.8148976 0.2730779 -0.7495766 0.4535097
kable(summary(model_latency_trim)$coefficients)
coef exp(coef) se(coef) z Pr(>|z|)
PopulationAL -0.2649116 0.7672738 0.1818483 -1.4567728 0.1451791
Predpred+ 0.5589274 1.7487958 0.1917097 2.9154886 0.0035513
Tutorsame 0.0146556 1.0147635 0.1353744 0.1082594 0.9137899
Waterpred+ -0.2400396 0.7865967 0.1357832 -1.7678153 0.0770918
PopulationAL:Predpred+ -0.6018387 0.5478035 0.2742801 -2.1942482 0.0282176

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(lme4)

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.5385)  ( log )
## Formula: Attempts ~ (Population + Tutor_pop + Pred + Water) + (1 | Mom_ID)
##    Data: subset(data_forage, Sex == "F")
## 
##      AIC      BIC   logLik deviance df.resid 
##   2535.0   2561.1  -1260.5   2521.0      302 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.7303 -0.6464 -0.3908  0.3039  5.2147 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Mom_ID (Intercept) 0.01928  0.1388  
## Number of obs: 309, groups:  Mom_ID, 24
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.5868     0.1764  20.331  < 2e-16 ***
## PopulationAL  -1.0159     0.1752  -5.797 6.75e-09 ***
## Tutor_popAL   -0.2398     0.1611  -1.489  0.13653    
## Predpred+      0.4445     0.1628   2.730  0.00634 ** 
## Waterpred+    -0.1872     0.1617  -1.158  0.24697    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PpltAL Ttr_AL Prdpr+
## PopulatinAL -0.449                     
## Tutor_popAL -0.423  0.004              
## Predpred+   -0.378 -0.102  0.097       
## Waterpred+  -0.352  0.049 -0.060 -0.113

Mother ID variance is very small. There is an effect of population (AL has more attempts), and an effect of rearing water (pred+ pecked more than pred-)

Let’s redo it for males.

library(lme4)

model_pecks <- glmer.nb(Attempts ~ (Population + Tutor_pop + Pred + Water) + (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.5845)  ( log )
## Formula: Attempts ~ (Population + Tutor_pop + Pred + Water) + (1 | Mom_ID)
##    Data: subset(data_forage, Sex == "M")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1909.0   1933.9   -947.5   1895.0      251 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.7545 -0.6429 -0.3772  0.1818  4.4210 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Mom_ID (Intercept) 0.2143   0.4629  
## Number of obs: 258, groups:  Mom_ID, 24
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.92348    0.25747  11.354   <2e-16 ***
## PopulationAL -0.53618    0.27634  -1.940   0.0523 .  
## Tutor_popAL   0.17609    0.18004   0.978   0.3281    
## Predpred+    -0.08824    0.17711  -0.498   0.6184    
## Waterpred+   -0.11534    0.17571  -0.656   0.5116    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PpltAL Ttr_AL Prdpr+
## PopulatinAL -0.630                     
## Tutor_popAL -0.327  0.028              
## Predpred+   -0.353  0.074  0.029       
## Waterpred+  -0.278  0.006 -0.166  0.099

The mother effect extremely small. The main effect of Population is weakened.

Let’s look at latency.

library(coxme)

model_latency <- coxme(time2peck ~ Population + Water + Pred + Tutor_pop + (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 = 266, 309 (1 observation deleted due to missingness)
##   Iterations= 10 43 
##                     NULL Integrated    Fitted
## Log-likelihood -1344.855  -1320.691 -1314.282
## 
##                   Chisq   df          p   AIC   BIC
## Integrated loglik 48.33 5.00 3.0424e-09 38.33 20.41
##  Penalized loglik 61.15 9.36 1.1628e-09 42.43  8.89
## 
## Model:  time2peck ~ Population + Water + Pred + Tutor_pop + (1 | Mom_ID) 
## Fixed coefficients
##                      coef exp(coef)  se(coef)     z       p
## PopulationAL -0.728498296 0.4826332 0.1556517 -4.68 2.9e-06
## Waterpred+   -0.135832279 0.8729890 0.1244707 -1.09 2.8e-01
## Predpred+     0.003804666 1.0038119 0.1261715  0.03 9.8e-01
## Tutor_popAL  -0.338941299 0.7125243 0.1277746 -2.65 8.0e-03
## 
## Random effects
##  Group  Variable  Std Dev   Variance 
##  Mom_ID Intercept 0.1902467 0.0361938
x
PopulationAL -0.7284983
Waterpred+ -0.1358323
Predpred+ 0.0038047
Tutor_popAL -0.3389413

Now for males:

library(coxme)

model_latency <- coxme(time2peck ~ Population + Water + Pred + Tutor_pop + Population:Pred + (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 = 220, 258 (2 observations deleted due to missingness)
##   Iterations= 5 27 
##                     NULL Integrated    Fitted
## Log-likelihood -1075.391  -1062.626 -1062.612
## 
##                   Chisq   df          p   AIC   BIC
## Integrated loglik 25.53 6.00 0.00027234 13.53 -6.83
##  Penalized loglik 25.56 5.02 0.00011026 15.53 -1.49
## 
## Model:  time2peck ~ Population + Water + Pred + Tutor_pop + Population:Pred +      (1 | Mom_ID) 
## Fixed coefficients
##                               coef exp(coef)  se(coef)     z      p
## PopulationAL           -0.26821863 0.7647406 0.1821960 -1.47 0.1400
## Waterpred+             -0.24415619 0.7833653 0.1364012 -1.79 0.0730
## Predpred+               0.55743518 1.7461881 0.1917012  2.91 0.0036
## Tutor_popAL            -0.04606092 0.9549838 0.1360824 -0.34 0.7400
## PopulationAL:Predpred+ -0.60139981 0.5480439 0.2741814 -2.19 0.0280
## 
## Random effects
##  Group  Variable  Std Dev     Variance   
##  Mom_ID Intercept 0.009023857 0.000081430
x
PopulationAL -0.2682186
Waterpred+ -0.2441562
Predpred+ 0.5574352
Tutor_popAL -0.0460609
PopulationAL:Predpred+ -0.6013998

Weight effects

We here repeat the above analyses but with Weight as a random effect, in order to evaluate potential genetic effects on foraging behavior

library(lme4)

model_pecks <- glmer.nb(Attempts ~ (Population + Tutor_pop + Pred + Water) + (1|Weight), subset(data_forage, Sex == 'F'))

summary(model_pecks)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(0.5412)  ( log )
## Formula: Attempts ~ (Population + Tutor_pop + Pred + Water) + (1 | Weight)
##    Data: subset(data_forage, Sex == "F")
## 
##      AIC      BIC   logLik deviance df.resid 
##   2535.2   2561.3  -1260.6   2521.2      302 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.7323 -0.6496 -0.3837  0.2804  5.2220 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Weight (Intercept) 0.03012  0.1736  
## Number of obs: 309, groups:  Weight, 85
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.5886     0.1731  20.728  < 2e-16 ***
## PopulationAL  -1.0439     0.1752  -5.960 2.53e-09 ***
## Tutor_popAL   -0.2294     0.1616  -1.419  0.15591    
## Predpred+      0.4515     0.1641   2.751  0.00594 ** 
## Waterpred+    -0.1777     0.1610  -1.103  0.26987    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PpltAL Ttr_AL Prdpr+
## PopulatinAL -0.323                     
## Tutor_popAL -0.451 -0.012              
## Predpred+   -0.398 -0.049  0.066       
## Waterpred+  -0.371  0.081 -0.075 -0.134

Let’s redo it for males.

library(lme4)

model_pecks <- glmer.nb(Attempts ~ (Population + Tutor_pop + Pred + Water) + (1|Weight), subset(data_forage, Sex == 'M'))
## Warning in optTheta(g1, interval = interval, tol = tol, verbose = verbose, :
## Model failed to converge with max|grad| = 0.0025586 (tol = 0.002, component 1)
summary(model_pecks)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(0.5999)  ( log )
## Formula: Attempts ~ (Population + Tutor_pop + Pred + Water) + (1 | Weight)
##    Data: subset(data_forage, Sex == "M")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1910.1   1935.0   -948.1   1896.1      251 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.7658 -0.6464 -0.3998  0.2737  4.6450 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Weight (Intercept) 0.2365   0.4864  
## Number of obs: 258, groups:  Weight, 58
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.861574   0.029001  98.672  < 2e-16 ***
## PopulationAL -0.539259   0.029177 -18.482  < 2e-16 ***
## Tutor_popAL   0.063681   0.029325   2.172   0.0299 *  
## Predpred+    -0.001516   0.029262  -0.052   0.9587    
## Waterpred+   -0.114994   0.028232  -4.073 4.64e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PpltAL Ttr_AL Prdpr+
## PopulatinAL -0.024                     
## Tutor_popAL -0.267 -0.005              
## Predpred+   -0.017 -0.252 -0.005       
## Waterpred+  -0.025 -0.012 -0.010 -0.008
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0025586 (tol = 0.002, component 1)

Let’s look at latency.

library(coxme)

model_latency <- coxme(time2peck ~ Population + Water + Pred + Tutor_pop + (1|Weight), 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 = 266, 309 (1 observation deleted due to missingness)
##   Iterations= 6 28 
##                     NULL Integrated    Fitted
## Log-likelihood -1344.855  -1319.825 -1298.743
## 
##                   Chisq   df          p   AIC    BIC
## Integrated loglik 50.06  5.0 1.3469e-09 40.06  22.14
##  Penalized loglik 92.22 22.8 2.5883e-10 46.63 -35.07
## 
## Model:  time2peck ~ Population + Water + Pred + Tutor_pop + (1 | Weight) 
## Fixed coefficients
##                       coef exp(coef)  se(coef)     z       p
## PopulationAL -0.8508169331 0.4270659 0.1414466 -6.02 1.8e-09
## Waterpred+   -0.1364429652 0.8724561 0.1257904 -1.08 2.8e-01
## Predpred+     0.0002519219 1.0002520 0.1344532  0.00 1.0e+00
## Tutor_popAL  -0.3769576247 0.6859451 0.1374750 -2.74 6.1e-03
## 
## Random effects
##  Group  Variable  Std Dev   Variance 
##  Weight Intercept 0.3249550 0.1055957
x
PopulationAL -0.8508169
Waterpred+ -0.1364430
Predpred+ 0.0002519
Tutor_popAL -0.3769576

Now for males:

library(coxme)

model_latency <- coxme(time2peck ~ Population + Water + Pred + Tutor_pop + Population:Pred + (1|Weight), 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 = 220, 258 (2 observations deleted due to missingness)
##   Iterations= 22 113 
##                     NULL Integrated    Fitted
## Log-likelihood -1075.391  -1062.467 -1056.891
## 
##                   Chisq    df          p   AIC    BIC
## Integrated loglik 25.85  6.00 2.3765e-04 13.85  -6.51
##  Penalized loglik 37.00 10.14 6.2658e-05 16.72 -17.68
## 
## Model:  time2peck ~ Population + Water + Pred + Tutor_pop + Population:Pred +      (1 | Weight) 
## Fixed coefficients
##                               coef exp(coef)  se(coef)     z      p
## PopulationAL           -0.29333276 0.7457739 0.1885596 -1.56 0.1200
## Waterpred+             -0.24739385 0.7808331 0.1373494 -1.80 0.0720
## Predpred+               0.56368656 1.7571384 0.1997460  2.82 0.0048
## Tutor_popAL            -0.04119972 0.9596375 0.1397144 -0.29 0.7700
## PopulationAL:Predpred+ -0.58818957 0.5553318 0.2834350 -2.08 0.0380
## 
## Random effects
##  Group  Variable  Std Dev    Variance  
##  Weight Intercept 0.17007777 0.02892645
x
PopulationAL -0.2933328
Waterpred+ -0.2473938
Predpred+ 0.5636866
Tutor_popAL -0.0411997
PopulationAL:Predpred+ -0.5881896