We will here (re)analyze the foraging data using planned comparisons and zero inflated models. Planned comparisons allow us to test specific hypotheses that involve different comparisons (or linear combinations of comparisons) among treatments.

load(here:::here('data/data_behav_SS.Rda'))

But first, let’s look at some bigger picture data. HP-LP comparisons in both drainages for sneak, sigmoid, and proportion of sneaks. This means we will NOT use the cross treatments.

data_hplp <- data_behav_SS %>%
  mutate(treatment = factor(paste(F_pred, M_pred, sep = ''))) %>%
  filter(treatment %in% c("HighHigh", "LowLow"))

data_treat <- data_behav_SS %>%
  mutate(treatment = factor(paste(F_pred, M_pred, sep = '')))

Let’s do a quick check for zero inflation on sigmoid and skeak.

p_sneak <- ggplot(data = data_hplp, aes(x = Sneak)) +
  geom_histogram(aes(fill = treatment), alpha = 0.5, position = 'identity', bins = 30) +
  theme_bw()

p_sigm <- ggplot(data = data_hplp, aes(x = Sigmoid)) +
  geom_histogram(aes(fill = treatment), alpha = 0.5, position = 'identity', bins = 30) +
  theme_bw()

grid.arrange(p_sneak, p_sigm, nrow = 1)

ARIPO

We will first look at population differences in sneaking population by comparing males from LP vs HP origin regardless of the treatment they are in.

model_aripo_sigmo <- glmer.nb(Sigmoid ~ Block + F_pred + M_pred +
                                (1|Mesocosm_unique),
                        data = subset(data_behav_SS,
                                      Drainage == 'Aripo'))

kable(summary(model_aripo_sigmo)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.9565082 0.5474799 1.7471113 0.0806180
Block3 -0.8111451 0.6769069 -1.1983112 0.2307959
Block4 -0.2009999 0.5904307 -0.3404292 0.7335333
F_predLow 0.1882997 0.4943265 0.3809218 0.7032613
M_predLow 0.0168325 0.4963089 0.0339153 0.9729447 
model_aripo_sneak <- glmer.nb(Sneak ~ Block + F_pred + M_pred +
                                (1|Mesocosm_unique),
                        data = subset(data_behav_SS,
                                      Drainage == 'Aripo'))

kable(summary(model_aripo_sneak)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.8633368 0.3922769 4.7500547 0.0000020
Block3 -0.2008741 0.4042547 -0.4968998 0.6192597
Block4 -0.4112512 0.4324432 -0.9509949 0.3416070
F_predLow 0.0885014 0.3605116 0.2454885 0.8060782
M_predLow -0.5881940 0.3319361 -1.7720094 0.0763930 
model_aripo_p <- glmer(cbind(Sneak, Sigmoid) ~ Block + F_pred + M_pred + 
                       (1|Mesocosm_unique),
                       family = 'binomial',
                       data = subset(data_behav_SS, Drainage == 'Aripo'))

kable(summary(model_aripo_p)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.9550840 0.3734605 2.5573898 0.0105461
Block3 0.0479255 0.4463776 0.1073654 0.9144991
Block4 -0.3374621 0.4031549 -0.8370531 0.4025627
F_predLow -0.1985380 0.3597659 -0.5518533 0.5810488
M_predLow -0.6778109 0.3520597 -1.9252728 0.0541952

Aripo fish show no population (treatment) differences in the amount of sneaking nor sigmoid display, however, out of those, there is an almost significantly lower proportion of sneaking in LP males.

We plot it to see the effect size.

sneak_summary <- data_treat %>%
  filter(Drainage == 'Aripo') %>%
  mutate(p_sneak = Sneak / (Sneak + Sigmoid)) %>%
  drop_na(p_sneak, M_Pop) %>%
  group_by(M_Pop) %>%
  summarize(mean_sneak = mean(p_sneak, na.rm = T),
            se = sqrt(mean_sneak * (1-mean_sneak)/length(p_sneak)))

sneak_plot <- sneak_summary %>%
  ggplot(aes(x =M_Pop, y = mean_sneak, col = M_Pop)) +
    geom_point(size = 3) +
    geom_errorbar(aes(ymin = mean_sneak - se, ymax = mean_sneak + se),
                  width = 0.2) + labs(y = 'Proportion Sneaking',
                                      x = 'Population') +
   scale_x_discrete(labels=c("HP", "LP")) +
   theme_bw() +
   theme(legend.position = "none")
   
sneak_plot

GUANAPO

We now repeat the same for the Guanapo drainage

model_guanapo_sigmo <- glmer.nb(Sigmoid ~ Block + F_pred + M_pred +
                                (1|Mesocosm_unique),
                        data = subset(data_behav_SS,
                                      Drainage == 'Guanapo'))

kable(summary(model_guanapo_sigmo)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.1819834 0.3228401 0.5636950 0.5729617
Block2 0.9759657 0.3762047 2.5942412 0.0094800
Block4 1.5389460 0.3658301 4.2067232 0.0000259
F_predLow 0.6392626 0.2984273 2.1421052 0.0321850
M_predLow -0.0974818 0.3001196 -0.3248099 0.7453250 
model_guanapo_sneak <- glmer.nb(Sneak ~ Block + F_pred * M_pred +
                                (1|Mesocosm_unique),
                        data = subset(data_behav_SS,
                                      Drainage == 'Guanapo'))

kable(summary(model_guanapo_sneak)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.5484831 0.3072854 1.784930 0.0742726
Block2 0.2852382 0.2729132 1.045161 0.2959486
Block4 0.4572638 0.2953369 1.548279 0.1215552
F_predLow 1.3571763 0.3489016 3.889854 0.0001003
M_predLow 1.4539009 0.3563183 4.080343 0.0000450
F_predLow:M_predLow -1.0583053 0.4765099 -2.220952 0.0263542 
model_guanapo_p <- glmer(cbind(Sneak, Sigmoid) ~ Block + F_pred + M_pred + 
                       (1|Mesocosm_unique),
                       family = 'binomial',
                       data = subset(data_behav_SS, Drainage == 'Guanapo'))

kable(summary(model_guanapo_p)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.7243305 0.3397073 2.1322190 0.0329888
Block2 -0.5588947 0.3609522 -1.5483901 0.1215284
Block4 -1.0452381 0.3552916 -2.9419159 0.0032619
F_predLow 0.1088546 0.2886788 0.3770788 0.7061151
M_predLow 0.8521738 0.2894191 2.9444287 0.0032355

Guanapo fish show a significantly higher counts and proportion of sneaking in LP males, The opposite to Aripo!

We plot it to see the effect size.

sneak_summary <- data_treat %>%
  filter(Drainage == 'Guanapo') %>%
  mutate(p_sneak = Sneak / (Sneak + Sigmoid)) %>%
  drop_na(p_sneak, M_Pop) %>%
  group_by(M_Pop) %>%
  summarize(mean_sneak = mean(p_sneak, na.rm = T),
            se = sqrt(mean_sneak * (1-mean_sneak)/length(p_sneak)))

sneak_plot <- sneak_summary %>%
  ggplot(aes(x =M_Pop, y = mean_sneak, col = M_Pop)) +
    geom_point(size = 3) +
    geom_errorbar(aes(ymin = mean_sneak - se, ymax = mean_sneak + se),
                  width = 0.2) + labs(y = 'Proportion Sneaking',
                                      x = 'Population') +
   scale_x_discrete(labels=c("HP", "LP")) +
   theme_bw() +
   theme(legend.position = "none")
   
sneak_plot

FORAGING BEHAVIOR

ARIPO

model_peck_aripo <- mixed_model(Peck ~ Block + F_pred + M_pred, 
                        zi_fixed = ~ Block + F_pred * M_pred,
                        random = ~ 1|Mesocosm_unique,
                        family = zi.negative.binomial(),
                        penalized = FALSE,
                        control = list(iter_EM = 0, optimizer = 'nlminb'),
                        data = subset(data_behav_SS,
                                      Drainage == 'Aripo'))

kable(summary(model_peck_aripo)$coef_table)
Estimate Std.Err z-value p-value
(Intercept) 2.935830 1.0336006 2.840391 0.0045058
Block3 1.889040 0.8703676 2.170394 0.0299770
Block4 -2.391675 0.9602092 -2.490786 0.0127461
F_predLow -2.849707 0.7899740 -3.607343 0.0003093
M_predLow 2.035173 0.7584094 2.683476 0.0072861 
kable(summary(model_peck_aripo)$coef_table_zi)
Estimate Std.Err z-value p-value
(Intercept) -0.3781031 1.221774 -0.3094705 0.7569636
Block3 8.0695998 4.729208 1.7063321 0.0879463
Block4 -10.9051210 161.663934 -0.0674555 0.9462191
F_predLow -8.2678644 5.031186 -1.6433230 0.1003161
M_predLow -3.2623717 3.383276 -0.9642641 0.3349135
F_predLow:M_predLow 3.0953504 3.770199 0.8210044 0.4116438

Aripo females with HP males peck significantly less times (when they peck) than those with LP males. Overall, Females from Aripo LP environments forage less than HP.

GUANAPO

model_peck_guanapo <- mixed_model(Peck ~ Block + F_pred + M_pred, 
                        zi_fixed = ~ Block + F_pred + M_pred,
                        random = ~ 1|Mesocosm_unique,
                        family = zi.negative.binomial(),
                        penalized = FALSE,
                        control = list(iter_EM = 0, optimizer = 'nlminb'),
                        data = subset(data_behav_SS,
                                      Drainage == 'Guanapo'))

kable(summary(model_peck_guanapo)$coef_table)
Estimate Std.Err z-value p-value
(Intercept) 2.5949602 1.2595402 2.0602441 0.0393752
Block2 -0.8209976 1.5473926 -0.5305684 0.5957179
Block4 -1.4508394 1.3614286 -1.0656742 0.2865709
F_predLow 1.3469878 0.6371119 2.1142090 0.0344974
M_predLow -0.4726668 0.7225364 -0.6541771 0.5129977 
kable(summary(model_peck_guanapo)$coef_table_zi)
Estimate Std.Err z-value p-value
(Intercept) 3.4161017 1.1839309 2.8853894 0.0039093
Block2 -2.3480409 1.1617029 -2.0212060 0.0432584
Block4 -4.5628108 1.5998724 -2.8519842 0.0043447
F_predLow -0.7468000 0.8513961 -0.8771475 0.3804065
M_predLow 0.4374854 0.8413805 0.5199615 0.6030904

Guanapo LP females eat more (when they eat) than HP. Nothing else is significant.

Finally, we want to analyze the effect of sexual behaviors on foraging.

data_behav_SS <- data_behav_SS %>%
  mutate(p_sneak = Sneak/(Sigmoid+Sneak))

model_peck <- mixed_model(Peck ~ Block + Drainage * Sneak + F_pred +Sigmoid, 
                        zi_fixed = ~ Block + Drainage * Sneak  + F_pred +Sigmoid,
                        random = ~ 1|Mesocosm_unique,
                        family = zi.negative.binomial(),
                        penalized = FALSE,
                        control = list(iter_EM = 0, optimizer = 'nlminb'),
                        data = data_behav_SS)

kable(summary(model_peck)$coef_table)
Estimate Std.Err z-value p-value
(Intercept) 3.3753594 0.9566007 3.5284935 0.0004179
Block2 -1.2797812 1.4441471 -0.8861848 0.3755180
Block3 2.1041217 1.3573097 1.5502148 0.1210900
Block4 -1.4529371 0.9801914 -1.4822993 0.1382607
DrainageGuanapo -0.9599715 1.1262008 -0.8523982 0.3939931
Sneak -0.1198408 0.0621761 -1.9274424 0.0539245
F_predLow -1.1320942 0.9283714 -1.2194411 0.2226768
Sigmoid 0.0104113 0.0543235 0.1916541 0.8480132
DrainageGuanapo:Sneak 0.2311706 0.0874141 2.6445451 0.0081801 
kable(summary(model_peck)$coef_table_zi)
Estimate Std.Err z-value p-value
(Intercept) 3.9782013 1.4262996 2.7891764 0.0052842
Block2 -0.7410771 1.1955793 -0.6198477 0.5353581
Block3 3.3033753 2.1240886 1.5551966 0.1198993
Block4 -1.9465670 1.1736976 -1.6584911 0.0972184
DrainageGuanapo -0.9521430 1.4985760 -0.6353652 0.5251903
Sneak -1.1152046 0.4600155 -2.4242760 0.0153389
F_predLow -0.6314906 1.1216988 -0.5629770 0.5734506
Sigmoid -0.5201838 0.2172164 -2.3947728 0.0166307
DrainageGuanapo:Sneak 1.2211883 0.5315127 2.2975712 0.0215862

There is a negative effect of Sneaking on Foraging, but only in the Aripo, not the Guanapo. This is true of the probability of foraging, but not the counts. The counts (given foraging) show an almost significant negative effect of sneaking (p = 0.054) in Aripo but not in Guanapo. There is also a negative effect of Sigmoid on Foraging (prob, not counts).

ARIPO

model_peck_aripo <- mixed_model(Peck ~ Block  + F_pred + Sneak + Sigmoid, 
                        zi_fixed = ~ Block + F_pred + Sneak  + Sigmoid,
                        random = ~ 1|Mesocosm_unique,
                        family = zi.negative.binomial(),
                        penalized = FALSE,
                        control = list(iter_EM = 0, optimizer = 'nlminb'),
                        data = subset(data_behav_SS, Drainage == 'Aripo'))

kable(summary(model_peck_aripo)$coef_table)
Estimate Std.Err z-value p-value
(Intercept) 2.8219709 1.2651265 2.2305839 0.0257087
Block3 0.7379403 1.9182888 0.3846868 0.7004695
Block4 -0.7031144 1.5141241 -0.4643704 0.6423824
F_predLow -0.0684277 1.9798240 -0.0345625 0.9724286
Sneak -0.1028260 0.0652064 -1.5769313 0.1148113
Sigmoid -0.0046730 0.0904817 -0.0516462 0.9588106 
kable(summary(model_peck_aripo)$coef_table_zi)
Estimate Std.Err z-value p-value
(Intercept) 1.2935160 1.2999449 0.9950545 0.3197098
Block3 2.2714497 2.3347552 0.9728856 0.3306102
Block4 1.3683606 1.5834287 0.8641757 0.3874914
F_predLow 2.6347179 2.3777793 1.1080582 0.2678367
Sneak -0.5860189 0.2757735 -2.1250005 0.0335866
Sigmoid -0.6217988 0.2432249 -2.5564764 0.0105738

With respect to the prob of foraging. Both Sneak and Sigmoid affect Foraging negatively in Aripo. This effect is not apparent on the counts.

GUANAPO

model_peck_guanapo <- mixed_model(Peck ~ Block + F_pred + Sneak + Sigmoid, 
                        zi_fixed = ~ Block + F_pred + Sneak + Sigmoid,
                        random = ~ 1|Mesocosm_unique,
                        family = zi.negative.binomial(),
                        penalized = FALSE,
                        control = list(iter_EM = 0, optimizer = 'nlminb'),
                        data = subset(data_behav_SS, Drainage == 'Guanapo'))

kable(summary(model_peck_guanapo)$coef_table)
Estimate Std.Err z-value p-value
(Intercept) 2.5092047 1.4280570 1.7570760 0.0789049
Block2 -1.6942229 1.6189962 -1.0464650 0.2953464
Block4 -1.7620220 1.5372497 -1.1462172 0.2517053
F_predLow 1.5579419 1.1342316 1.3735660 0.1695765
Sneak 0.0035990 0.0753018 0.0477941 0.9618803
Sigmoid 0.0094285 0.0538091 0.1752211 0.8609059 
kable(summary(model_peck_guanapo)$coef_table_zi)
Estimate Std.Err z-value p-value
(Intercept) 4.5208712 1.7335194 2.6079149 0.0091096
Block2 -2.7768245 2.0223469 -1.3730703 0.1697305
Block4 -5.2614623 3.0009285 -1.7532782 0.0795543
F_predLow 1.9673057 2.3758253 0.8280515 0.4076413
Sneak 0.1022289 0.1385019 0.7381051 0.4604506
Sigmoid -0.7557163 0.4543522 -1.6632829 0.0962558

No effect of sexual behaviors on foraging for Guanapo.