Analysis for male behavior manuscript

I. Relative importance of predictors

Looking at the percentage of contribution of each factor to male mating strategy.

mod <- lm(logit(sigmoid_prop)  ~ pop*(tutor_pop + pred_trt + assay_water),
           data = data)

relative_importance <-   relaimpo::calc.relimp(mod, rela=T)

relative_importance$lmg %>% kable(digits = 3)

	x
tutor_pop	0.638
pop:tutor_pop	0.081
pop	0.108
pred_trt	0.009
assay_water	0.001
pop:pred_trt	0.152
pop:assay_water	0.012

# quick plot - but I need more customization
# plot(relative_importance)

Figure 1

relative_importance_df <- 
  data.frame(
    variable=c("social treatment", 
               "population × social",
               "population origin",
               "predator treatment",
               "assay water",
               "population × predator",
               "population × assay water"),
    importance= relative_importance$lmg
    )




ggplot(relative_importance_df, 
       aes(x = reorder(variable, importance), y = importance)) +

  geom_col(fill = "grey70") + 

  labs(y = "Percent explained",
       x = "Variables") +
  
  scale_y_continuous(labels = scales::percent) +
  
  coord_flip() +  
  theme_classic(base_size = 15)

II. Mating strategy among treatments

HP males

mod <-glmer(cbind(sigmoid, sneak)  ~ tutor_pop+assay_water+pred_trt+ 
              (1|rowID)+(1|ID) + (1|mom_ID), 
            data = data %>% filter(pop=="AH"),
            family = binomial())

Anova(mod, type=3) %>% kable(digits=3)

	Chisq	Df	Pr(>Chisq)
(Intercept)	0.583	1	0.445
tutor_pop	11.266	2	0.004
assay_water	0.181	1	0.671
pred_trt	7.489	1	0.006

summary(mod)$coefficients %>% kable(digits=3)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-0.393	0.515	-0.763	0.445
tutor_popAH	-1.199	0.455	-2.637	0.008
tutor_popAL	0.278	0.454	0.611	0.541
assay_waterplus	-0.062	0.146	-0.425	0.671
pred_trtY	0.999	0.365	2.737	0.006

emmeans(mod, ~tutor_pop) %>% pairs %>% kable(digits=3)

contrast	estimate	SE	df	z.ratio	p.value
solo - AH	1.199	0.455	Inf	2.637	0.023
solo - AL	-0.278	0.454	Inf	-0.611	0.814
AH - AL	-1.477	0.470	Inf	-3.143	0.005

mod_mm <- emmeans(mod, ~tutor_pop|pred_trt:assay_water, type = "response") %>% summary() %>% tibble

AH_prop <- ggplot() +
  geom_line(data = mod_mm, 
            aes(x = tutor_pop, y = prob, group = interaction(pred_trt,assay_water),
                lty = assay_water, color = pred_trt),
            position=position_dodge(width = 0.2))+
  geom_point(data = mod_mm, 
             aes(x = tutor_pop, y = prob, group = interaction(pred_trt,assay_water),
                 color = pred_trt, shape = assay_water),
             size = 3, position=position_dodge(width = 0.2))+
  geom_errorbar(data = mod_mm, 
            aes(x = tutor_pop, ymin = prob - SE, ymax = prob + SE,
            group = interaction(pred_trt,assay_water), color = pred_trt),
            width = 0.2, position=position_dodge(width = 0.2)) +
  labs(title = "HP mating assay",
       x = "Social treatment", y = "Sigmoid proportion")+
  
  scale_shape_manual(name = "Assay water", labels = c("pred-", "pred+"), values = c(16, 1)) +
  scale_color_manual(name = "Predator treatment", labels = c("pred-","pred+"),
                     values = c("lightsalmon1","SteelBlue")) +
  scale_linetype_manual(name = "Assay water", labels = c("pred-", "pred+"), values = c(1, 2)) +

  scale_y_continuous(limits = c(0,1))+
  scale_x_discrete(labels = c("solo","HP tutor","LP tutor"))+
  
  annotate("text", label = c("a","b","a"), y = 1, x = c(1, 2, 3), size = 4)+
  
  
  theme_bw(base_size = 15) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

AH_prop

LP males

mod <-glmer(cbind(sigmoid, sneak)  ~ tutor_pop+assay_water*pred_trt+ 
              (1|rowID)+(1|ID) + (1|mom_ID) , 
            data = data %>% filter(pop=="AL"),
            family = binomial())

Anova(mod, type=3) %>% kable(digits=3)

	Chisq	Df	Pr(>Chisq)
(Intercept)	1.321	1	0.250
tutor_pop	18.454	2	0.000
assay_water	7.290	1	0.007
pred_trt	2.963	1	0.085
assay_water:pred_trt	3.886	1	0.049

summary(mod)$coefficients %>% kable(digits=3)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	0.396	0.344	1.149	0.250
tutor_popAH	-1.489	0.350	-4.259	0.000
tutor_popAL	-0.908	0.351	-2.587	0.010
assay_waterplus	-0.437	0.162	-2.700	0.007
pred_trtY	-0.552	0.321	-1.721	0.085
assay_waterplus:pred_trtY	0.473	0.240	1.971	0.049

emmeans(mod, ~tutor_pop) %>% pairs %>% kable(digits=3)

contrast	estimate	SE	df	z.ratio	p.value
solo - AH	1.489	0.350	Inf	4.259	0.000
solo - AL	0.908	0.351	Inf	2.587	0.026
AH - AL	-0.581	0.353	Inf	-1.647	0.226

emmeans(mod, ~ assay_water*pred_trt) %>% 
  pairs(simple = "each", combine = TRUE) %>% kable(digits=3)

pred_trt	assay_water	contrast	estimate	SE	df	z.ratio	p.value
N	.	minus - plus	0.437	0.162	Inf	2.700	0.028
Y	.	minus - plus	-0.036	0.177	Inf	-0.201	1.000
.	minus	N - Y	0.552	0.321	Inf	1.721	0.341
.	plus	N - Y	0.079	0.323	Inf	0.244	1.000

mod_mm <- emmeans(mod, ~tutor_pop|assay_water:pred_trt, type = "response") %>% summary() %>% tibble

AL_prop <- ggplot() +
  geom_line(data = mod_mm, 
            aes(x = tutor_pop, y = prob, group = interaction(pred_trt,assay_water),
                lty = assay_water, color = pred_trt),
            position=position_dodge(width = 0.2))+
  geom_point(data = mod_mm, 
             aes(x = tutor_pop, y = prob, group = interaction(pred_trt,assay_water),
                 shape = assay_water, color = pred_trt),
             size = 3, position=position_dodge(width = 0.2))+
  geom_errorbar(data = mod_mm, 
            aes(x = tutor_pop, ymin = prob - SE, ymax = prob + SE,
            group = interaction(pred_trt,assay_water), color = pred_trt),
            width = 0.2, position=position_dodge(width = 0.2)) +
  labs(title = "LP mating assay",
       x = "Social treatment", y = "Sigmoid proportion")+
  
  scale_shape_manual(name = "Assay water", labels = c("pred-", "pred+"), values = c(16, 1)) +
  scale_linetype_manual(name = "Assay water", labels = c("pred-", "pred+"), values = c(1, 2)) +
  scale_color_manual(name = "Predator treatment", labels = c("pred-","pred+"),
                     values = c("lightsalmon1","SteelBlue")) +
  scale_y_continuous(limits = c(0,1))+
  scale_x_discrete(labels = c("solo","HP tutor","LP tutor"))+
  
  annotate("text", label = c("a","b","b"), y = 1, x = c(1, 2, 3), size = 4)+
  
  
  theme_bw(base_size = 15) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

AL_prop

Figure 2

egg::ggarrange(AH_prop + theme(legend.position = "none"),
               AL_prop,
               nrow = 1,
               labels = c("A","B"))

III. Tutor behavior correlation

HP males

mod <-glmer(cbind(sigmoid, sneak) ~ sigmoid_prop_tutor+
               (1|ID) , 
            data = data %>% filter(assay_water=="minus", pop=="AH"),
            family = binomial())

summary(mod)$coefficients %>% kable(digits=3)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-1.172	0.410	-2.859	0.004
sigmoid_prop_tutor	1.871	0.906	2.066	0.039

LP males

mod <-glmer(cbind(sigmoid, sneak) ~ sigmoid_prop_tutor+
               (1|ID) , 
            data = data %>% filter(assay_water=="minus", pop=="AL"),
            family = binomial())

summary(mod)$coefficients %>% kable(digits=3)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-1.238	0.344	-3.604	0.000
sigmoid_prop_tutor	0.533	0.637	0.837	0.403

Figure 3

data_bin <- data %>% filter(pop=="AH") %>%
  mutate(bins = cut(sigmoid_prop_tutor, breaks = 8)) %>%
  filter(!is.na(sigmoid_prop), !is.na(sigmoid_prop_tutor)) %>%
  group_by(bins) %>%
  summarise(mean = mean(sigmoid_prop), 
            sd = sd(sigmoid_prop),
            se = sd(sigmoid_prop)/ sqrt(length(sigmoid_prop))) %>%
  mutate(midpoint = seq(0.0625,0.9375,0.125))

HP_sex_cor <- ggplot(data %>% filter(pop=="AH") , aes(x = sigmoid_prop_tutor, y = sigmoid_prop)) +
  geom_point(alpha = 0.1)+
  geom_smooth(method=lm, alpha = 0.2, color = "SteelBlue", fill = "SteelBlue")+
  geom_errorbar(data = data_bin, 
                aes(x = midpoint, y = mean, ymin = mean - se, ymax = mean + se),
                width = 0, color = "black") +

  geom_point(data = data_bin, aes(x = midpoint, y = mean), color = "black", size = 2)+
  
  labs(title = "HP males", 
       x = "tutor sigmoid proportion",
       y = "sigmoid proportion")+

  
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio=1)

data_bin <- data %>% filter(pop=="AL") %>%
  mutate(bins = cut(sigmoid_prop_tutor, breaks = 8)) %>%
  filter(!is.na(sigmoid_prop), !is.na(sigmoid_prop_tutor)) %>%
  group_by(bins) %>%
  summarise(mean = mean(sigmoid_prop), 
            sd = sd(sigmoid_prop),
            se = sd(sigmoid_prop)/ sqrt(length(sigmoid_prop))) %>%
  mutate(midpoint = seq(0.0625,0.9375,0.125))

LP_sex_cor <- ggplot(data %>% filter(pop=="AL") , aes(x = sigmoid_prop_tutor, y = sigmoid_prop)) +
  geom_point(alpha = 0.1)+
  geom_smooth(method=lm, alpha = 0.2, color = "lightsalmon1", fill = "lightsalmon1")+
  geom_errorbar(data = data_bin, 
                aes(x = midpoint, y = mean, ymin = mean - se, ymax = mean + se),
                width = 0, color = "black") +

  geom_point(data = data_bin, aes(x = midpoint, y = mean), color = "black", size = 2)+
  
  labs(title = "LP males", 
       x = "tutor sigmoid proportion",
       y = "sigmoid proportion")+

  
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio=1)

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

HP vs LP male brain size

mod <- lmer(scale(BM) ~ scale(weight) + pop+  (1|mom_ID),
            data =  data %>% filter(assay_water=="minus"))
summary(mod)$coefficients %>% kable(digit = 3)

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	0.081	0.116	16.125	0.695	0.497
scale(weight)	0.612	0.075	136.995	8.148	0.000
popAL	-0.359	0.159	17.653	-2.259	0.037

IV. Path analysis - brain morphology

data_path <- data %>% filter(BM!="NA" ) %>% droplevels() %>%
  mutate(logit_sigmoid_prop = logit(sigmoid_prop))

data_path <- cbind(data_path, psych::dummy.code(data_path$tutor_pop)[,1:2])



mod <- '

#measurement model
brain =~ Tel_resid + Cb_resid + OB_resid + OT_resid + Hyp_resid


#direct effects
logit_sigmoid_prop ~ a*pred_trt
logit_sigmoid_prop ~ b*AH
logit_sigmoid_prop ~ c*AL
logit_sigmoid_prop ~ d*assay_water

#mediator
brain ~ e*pred_trt
brain ~ f*AH
brain ~ g*AL

BM_resid ~ h*pred_trt
BM_resid ~ i*AH
BM_resid ~ j*AL


#mediator to outcome
logit_sigmoid_prop ~ k*brain
logit_sigmoid_prop ~ l*BM_resid

#indirect effects
brain_pred := e*k
brain_AH := f*k 
brain_AL := g*k

BM_pred := h*l
BM_AH := i*l 
BM_AL := j*l

#total effects
tot_pred := e*k+ h*l + a
tot_AH := f*k + i*l + b
tot_AL := g*k + j*l + c


'

fitmod <- sem(mod, data = data_path %>% filter(pop=="AH"))

summary(fitmod, standardized=T)

## lavaan 0.6.14 ended normally after 101 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        24
## 
##                                                   Used       Total
##   Number of observations                           129         138
## 
## Model Test User Model:
##                                                       
##   Test statistic                                96.091
##   Degrees of freedom                                32
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   brain =~                                                              
##     Tel_resid         1.000                               0.083    0.775
##     Cb_resid          1.061    0.312    3.401    0.001    0.088    0.398
##     OB_resid         -0.306    0.258   -1.186    0.236   -0.025   -0.124
##     OT_resid          0.687    0.163    4.214    0.000    0.057    0.633
##     Hyp_resid         0.666    0.224    2.972    0.003    0.055    0.336
## 
## Regressions:
##                        Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   logit_sigmoid_prop ~                                                      
##     pred_trt   (a)        0.587    0.319    1.843    0.065    0.587    0.150
##     AH         (b)       -0.838    0.388   -2.157    0.031   -0.838   -0.203
##     AL         (c)        0.705    0.394    1.788    0.074    0.705    0.167
##     assay_watr (d)        0.164    0.309    0.529    0.597    0.164    0.042
##   brain ~                                                                   
##     pred_trt   (e)       -0.010    0.017   -0.597    0.550   -0.122   -0.061
##     AH         (f)        0.044    0.021    2.092    0.036    0.525    0.250
##     AL         (g)        0.038    0.021    1.795    0.073    0.459    0.214
##   BM_resid ~                                                                
##     pred_trt   (h)        0.057    0.022    2.591    0.010    0.057    0.222
##     AH         (i)        0.012    0.027    0.459    0.646    0.012    0.045
##     AL         (j)        0.020    0.027    0.746    0.456    0.020    0.074
##   logit_sigmoid_prop ~                                                      
##     brain      (k)        2.505    2.319    1.080    0.280    0.209    0.106
##     BM_resid   (l)        3.485    1.238    2.815    0.005    3.485    0.228
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .Tel_resid         0.005    0.002    2.893    0.004    0.005    0.399
##    .Cb_resid          0.041    0.006    7.399    0.000    0.041    0.842
##    .OB_resid          0.041    0.005    7.982    0.000    0.041    0.985
##    .OT_resid          0.005    0.001    5.186    0.000    0.005    0.599
##    .Hyp_resid         0.024    0.003    7.615    0.000    0.024    0.887
##    .logit_sgmd_prp    3.068    0.384    7.986    0.000    3.068    0.797
##    .BM_resid          0.016    0.002    8.031    0.000    0.016    0.948
##    .brain             0.007    0.002    3.436    0.001    0.940    0.940
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     brain_pred       -0.025    0.049   -0.524    0.600   -0.025   -0.006
##     brain_AH          0.109    0.113    0.970    0.332    0.109    0.027
##     brain_AL          0.096    0.102    0.935    0.350    0.096    0.023
##     BM_pred           0.199    0.104    1.906    0.057    0.199    0.051
##     BM_AH             0.043    0.094    0.453    0.651    0.043    0.010
##     BM_AL             0.071    0.098    0.721    0.471    0.071    0.017
##     tot_pred          0.761    0.320    2.374    0.018    0.761    0.194
##     tot_AH           -0.686    0.388   -1.769    0.077   -0.686   -0.167
##     tot_AL            0.872    0.397    2.198    0.028    0.872    0.207

fitmod <- sem(mod, data = data_path %>% filter(pop=="AL"))

summary(fitmod, standardized=T)

## lavaan 0.6.14 ended normally after 93 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        24
## 
##                                                   Used       Total
##   Number of observations                           130         144
## 
## Model Test User Model:
##                                                       
##   Test statistic                               138.781
##   Degrees of freedom                                32
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   brain =~                                                              
##     Tel_resid         1.000                               0.137    0.955
##     Cb_resid          0.602    0.126    4.768    0.000    0.082    0.435
##     OB_resid          0.201    0.159    1.262    0.207    0.027    0.115
##     OT_resid          0.387    0.060    6.471    0.000    0.053    0.595
##     Hyp_resid         0.686    0.114    5.995    0.000    0.094    0.549
## 
## Regressions:
##                        Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   logit_sigmoid_prop ~                                                      
##     pred_trt   (a)       -0.679    0.363   -1.871    0.061   -0.679   -0.193
##     AH         (b)       -1.110    0.381   -2.916    0.004   -1.110   -0.304
##     AL         (c)       -0.688    0.374   -1.841    0.066   -0.688   -0.189
##     assay_watr (d)       -0.275    0.294   -0.937    0.349   -0.275   -0.078
##   brain ~                                                                   
##     pred_trt   (e)       -0.144    0.022   -6.634    0.000   -1.051   -0.524
##     AH         (f)       -0.030    0.027   -1.085    0.278   -0.217   -0.104
##     AL         (g)       -0.016    0.027   -0.592    0.554   -0.118   -0.057
##   BM_resid ~                                                                
##     pred_trt   (h)        0.076    0.025    3.084    0.002    0.076    0.258
##     AH         (i)        0.056    0.031    1.796    0.073    0.056    0.183
##     AL         (j)        0.021    0.031    0.680    0.496    0.021    0.070
##   logit_sigmoid_prop ~                                                      
##     brain      (k)       -2.159    1.368   -1.578    0.115   -0.295   -0.168
##     BM_resid   (l)       -0.081    1.045   -0.077    0.939   -0.081   -0.007
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .Tel_resid         0.002    0.002    0.986    0.324    0.002    0.088
##    .Cb_resid          0.029    0.004    7.830    0.000    0.029    0.811
##    .OB_resid          0.056    0.007    8.053    0.000    0.056    0.987
##    .OT_resid          0.005    0.001    7.293    0.000    0.005    0.646
##    .Hyp_resid         0.020    0.003    7.519    0.000    0.020    0.698
##    .logit_sgmd_prp    2.796    0.348    8.039    0.000    2.796    0.906
##    .BM_resid          0.020    0.002    8.062    0.000    0.020    0.907
##    .brain             0.013    0.003    5.278    0.000    0.716    0.716
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     brain_pred        0.310    0.201    1.540    0.124    0.310    0.088
##     brain_AH          0.064    0.072    0.895    0.371    0.064    0.018
##     brain_AL          0.035    0.063    0.554    0.580    0.035    0.010
##     BM_pred          -0.006    0.080   -0.077    0.939   -0.006   -0.002
##     BM_AH            -0.005    0.059   -0.077    0.939   -0.005   -0.001
##     BM_AL            -0.002    0.022   -0.077    0.939   -0.002   -0.000
##     tot_pred         -0.375    0.298   -1.260    0.208   -0.375   -0.106
##     tot_AH           -1.051    0.378   -2.782    0.005   -1.051   -0.287
##     tot_AL           -0.655    0.376   -1.742    0.081   -0.655   -0.180

lavaanPlot::lavaanPlot(model = fitmod, coefs = TRUE)

V. Brain-Gonopodium-Behavior correlation

gonopodium-behavior

# rowID is added as a observation level random effect to help with overdispersion
mod <-glmer(cbind(sigmoid, sneak)  ~ gono_resid+pop+
              (1|rowID)+(1|ID)+(1|mom_ID), 
            data = data,
            family = binomial())
summary(mod)$coefficients %>% kable(digit = 3)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-0.271	0.294	-0.922	0.356
gono_resid	-5.558	1.884	-2.950	0.003
popAL	-0.639	0.417	-1.534	0.125

mod.fit <- emmip(mod, pop~gono_resid, cov.reduce = range, type = "response", plotit = F)

gono_behav <- ggplot()+
  geom_point(data = data ,
             aes(x = gono_resid, y = sigmoid_prop, color = pop), alpha = 0.3) +

  stat_summary_bin(data=data, 
                   aes(x = gono_resid, y = sigmoid_prop, color = pop),
                   fun.data = mean_se, 
                   bins=8, size =0.5, position = position_dodge(0.02))+
  
  geom_line(size = 1 ,data = mod.fit, aes(x = xvar, y = yvar, color = pop))+
  labs( x = "Residual gonopodium length", y = "Sigmoid proportion") +
  
  scale_color_manual("Population origin",
                     labels = c("HP","LP"),
                     values=c("SteelBlue", "lightsalmon1")) +
  
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom",
        aspect.ratio=1) 

gono_behav

brain size-behavior

# ID is added as a observation level random effect to help with overdispersion
mod <-glmer(cbind(sigmoid, sneak)  ~ pop*BM_resid+
              (1|rowID)+(1|ID) + (1|mom_ID), 
            data = data,
            family = binomial(),
            control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000)))

## boundary (singular) fit: see help('isSingular')

summary(mod)$coefficients %>% kable(digit = 3)

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-1.118	0.421	-2.658	0.008
popAL	0.374	0.520	0.719	0.472
BM_resid	4.098	1.763	2.324	0.020
popAL:BM_resid	-3.721	2.321	-1.603	0.109

emtrends(mod, ~pop, "BM_resid") %>% test

##  pop BM_resid.trend   SE  df z.ratio p.value
##  AH           4.098 1.76 Inf   2.324  0.0201
##  AL           0.376 1.51 Inf   0.249  0.8031

mod.fit <- emmip(mod, pop~BM_resid, cov.reduce = range, type = "response", plotit = F)

BM_behav <- ggplot()+
  geom_point(data = data,
             aes(x = BM_resid, y = sigmoid_prop, color = pop), alpha = 0.3) +
  geom_line(size = 1 ,data = mod.fit, aes(x = xvar, y = yvar, color = pop))+
  
  stat_summary_bin(data=data, 
                   aes(x = BM_resid, y = sigmoid_prop, color = pop),
                   fun.data = mean_se, 
                   bins=8, size =0.5, position = position_dodge(0.025))+
  
  labs( x = "Residual brain mass", y = "Sigmoid proportion") +
  
  scale_color_manual("Population origin",
                     labels = c("HP","LP"),
                     values=c("SteelBlue", "lightsalmon1")) +
  scale_x_continuous(limits =c(-0.2506015, 0.4115183)) +
  
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom",
        aspect.ratio=1) 

BM_behav

brain size-gonopodia

mod <-lmer(scale(BM)  ~ gono_resid + pop + weight +
              (1|mom_ID), 
            data = data %>% filter(assay_water=="minus") )

summary(mod)$coefficients %>% kable(digit = 3)

	Estimate	Std. Error	df	t value	Pr(>|t|)
(Intercept)	-2.404	0.350	131.603	-6.874	0.000
gono_resid	-1.836	0.994	136.471	-1.847	0.067
popAL	-0.369	0.168	17.287	-2.195	0.042
weight	20.339	2.600	136.588	7.823	0.000

mod <-lmer(BM_resid  ~ gono_resid+pop +  (1|mom_ID), 
            data = data %>% filter(assay_water=="minus") )
mod.fit <- emmip(mod, pop~gono_resid, cov.reduce = range, type = "response", plotit = F)

data_minus <- data %>% filter(assay_water=="minus", BM_resid !="NA")

BM_gono <- ggplot()+
  geom_point(data = data_minus,
             aes(x = gono_resid, y = BM_resid, color = pop), alpha = 0.3) +
  geom_line(size = 1 ,data = mod.fit, aes(x = xvar, y = yvar, color = pop))+
  
  stat_summary_bin(data=data_minus, 
                   aes(x = gono_resid, y = BM_resid, color = pop),
                   fun.data = mean_se, 
                   bins=8, size =0.5, position = position_dodge(0.025))+
  
  labs( x = "Residual gonopodium length", y = "Residual brain mass") +
  
  scale_color_manual("Population origin",
                     labels = c("HP","LP"),
                     values=c("SteelBlue", "lightsalmon1")) +
  
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom",
        aspect.ratio=1) 

BM_gono

Figure 5

egg::ggarrange(gono_behav + theme(legend.position = "none"),
               BM_behav, 
               BM_gono + theme(legend.position = "none"), nrow = 1,labels = c("A","B","C"))