Markdown _WBCs.Rmd

---
title: "WBC & Sociality MS"
output: html_document
date: "2026-01-12"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(tidyr)
library(magrittr)
library(dplyr)
library(tidyverse)
library(bisonR)
library(igraph)
library(brms)
library(beepr)
library(ggplot2)
library(bayestestR)
library(reshape2)
#Prior predictive checks
library(marginaleffects)
library(modelsummary)
#extract tables
library(sjPlot)
#Marginal effects
library(emmeans)
#Correlation matrix
library(corrplot)
library(scales)

```

#Load data
```{r}
#set directory:
setwd()

#Edgelists for groups of interest
cayo_edgelist = read.csv("cayo_edgelists.csv")
#WBC data and individual attributes
wbc_behav = read.csv("wbc_behav.csv")

```

#Generate Bison networks
```{r}

#Create list of groups analysed to loop through
group.year.list = c("F2015","V2015","R2015","KK2015","S2019")

#ADD SAMPLING EFFORT DYAD IN NUMBER OF FOCALS
#Only for S2019 (short-focal) compute effort by dividing dyad obs effort in minutes by 5, for all other group_years divided by 10
cayo_edgelist %<>% mutate(total_focals = ifelse(group.year == "S2019", ((focal.sampling.effort.duration.id.1*60)/5)+((focal.sampling.effort.duration.id.2*60)/5),
                                              ((focal.sampling.effort.duration.id.1*60)/10)+((focal.sampling.effort.duration.id.2*60)/10)))


#ADD EDGE WEIGHTS
cayo_edgelist %<>% mutate(grooming.proximity_events = aggregated.proximity.events.scan + grooming.bouts.focal,
                                dsi = grooming.proximity_events/total_focals)

####BISON NETWORKS FOR EACH GROUP YEAR####

#Convert individuals ids into factors
cayo_edgelist %<>% mutate(id.1 = factor(id.1),
                          id.2 = factor(id.2))

#PROXIMITY + GROOMING NETROWKS###
#Priors should be general for all group_years
priors <- get_default_priors("count_conjugate")#define priors
priors$edge <- "gamma(0.1,1)" # set less flat prior for proximity networks

for(i in 1:length(group.year.list)){
        #subset to group_year
        edgelist = cayo_edgelist[cayo_edgelist$group.year == group.year.list[i],]
        #run models for prox
        fit_bison <- bison_model(
                (grooming.proximity_events | total_focals) ~ dyad(id.1, id.2), 
                 data= edgelist, 
                 model_type = "count_conjugate",
                priors=priors
                )
        assign(paste0(group.year.list[i]),fit_bison)
      }

#Check if bison models make sense##
#bisonR::plot_predictions(S2019, num_draws=20)#replace group year

#remove unnecessary data
rm(edgelist,fit_bison,priors)
```

Define thresholds as minimum empirical edge weight (different from zero)
```{r}

#Loop to get thresholds
for(j in 1:length(group.year.list)){
    GROUPYEAR = group.year.list[j]#define group year
    #subset to group year of interest
    only_group = cayo_edgelist %>% filter(group.year == GROUPYEAR)
    #Define threshold
    thres = min(only_group$dsi[only_group$dsi != 0])
    #store values
    assign(paste0(GROUPYEAR,"_thres"), thres)
  }

rm(GROUPYEAR,only_group,thres)
```

#Compute network metrics
```{r}
#List with bison networks 
nets = group.year.list
thres = c("F2015_thres","V2015_thres","R2015_thres","KK2015_thres","S2019_thres")

for(k in 1:length(nets)){
    fit_edge = get(nets[k])
    num_draws = 100 #
    nodes = names(fit_edge$node_to_idx)
    #Generate df with nodes where data will be imputed
    df_nodal = data.frame(node_id = as.factor(nodes))
    #Set number of nodes and number of posterior samples
    num_nodes = fit_edge$num_nodes
    #Create a list of igraph networks from edgemodel to represent network posterior
    network <- bison_to_igraph(fit_edge, num_draws)   
    
    # Create an empty list for imputed versions of the dataframe
    imputed_data <- list()
    
    #Get threshold
    threshold = get(thres[k])
    #Loop through each posterior sample and calculate strength of top mixture.
  for (i in 1:num_draws) {
    #Add node label to dataframe
    imputed_data[[i]] <- df_nodal
    #make zero dyad non-seen interacting
    edge_exist =  1* (E(network[[i]])$weight >= threshold)
    E(network[[i]])$weight = E(network[[i]])$weight * edge_exist
    
    #Compute degree using strength (setting edge weight as 1 for all edges > 0)
    edges_deg = 1 * (E(network[[i]])$weight > 0)
    network_deg = network[[i]]
    E(network_deg)$weight <- edges_deg
    #Degree 
    imputed_data[[i]]$degree <- strength(network_deg)[match(df_nodal$node_id, names(strength(network_deg)))]
    
    #Compute weak/strong partners (upper quantile)
    #define threshold for strong/weak partners as upper quantile
    thres_str = round(quantile(E(network[[i]])$weight[E(network[[i]])$weight >= threshold], 0.75), digits = 4)
    edges_weak = 1 * (E(network[[i]])$weight < thres_str & E(network[[i]])$weight >= threshold)
    network_weak = network[[i]]
    E(network_weak)$weight <- edges_weak
    #Degree for weak partners
    imputed_data[[i]]$weak_Q <- strength(network_weak)[match(df_nodal$node_id, names(strength(network_weak)))]
    #Network does not need filtering
    edges_strong = 1 * (E(network[[i]])$weight >= thres_str)
    network_str = network[[i]]
    E(network_str)$weight <- E(network_str)$weight * edges_strong
    imputed_data[[i]]$top_strQ <- strength(network_str)[match(df_nodal$node_id, names(strength(network_str)))]
    
    #Compute weak/strong partners (10%/90%)
    #define threshold for strong/weak partners as upper quantile
    thres_str2 = round(quantile(E(network[[i]])$weight[E(network[[i]])$weight >= threshold], 0.9), digits = 4)
    edges_weak2 = 1 * (E(network[[i]])$weight < thres_str2 & E(network[[i]])$weight >= threshold)
    network_weak2 = network[[i]]
    E(network_weak2)$weight <- edges_weak2
    #Degree for weak partners
    imputed_data[[i]]$weak_10 <- strength(network_weak2)[match(df_nodal$node_id, names(strength(network_weak2)))]
    #Network does not need filtering
    edges_strong2 = 1 * (E(network[[i]])$weight >= thres_str2)
    network_str2 = network[[i]]
    E(network_str2)$weight <- E(network_str2)$weight * edges_strong2
    imputed_data[[i]]$top_str90 <- strength(network_str2)[match(df_nodal$node_id, names(strength(network_str2)))]
  
    #Save data
    assign(paste0(nets[k],"_bison"), imputed_data)
  }
}

rm(df_nodal,fit_edge,imputed_data, network,network_deg,network_weak,network_weak2,network_str,network_str2)

```

Merge group years and subset to keep only individuals with blood smear data
```{r}

#MERGE DATAFRAMES##
networks <- list()
ids = data.frame(id = as.character(c(F2015_bison[[1]]$node_id,KK2015_bison[[1]]$node_id,
                                     R2015_bison[[1]]$node_id,V2015_bison[[1]]$node_id, S2019_bison[[1]]$node_id)))


num_draws = 100
  for(i in 1:num_draws){
    networks[[i]]<- ids
    networks[[i]]$groupyear = c(rep("F2015", length(F2015_bison[[1]]$node_id)),
                                 rep("KK2015", length(KK2015_bison[[1]]$node_id)),
                                 rep("R2015", length(R2015_bison[[1]]$node_id)),
                                 rep("V2015", length(V2015_bison[[1]]$node_id)),
                                 rep("S2019", length(S2019_bison[[1]]$node_id)))

    networks[[i]]$deg = c(F2015_bison[[i]]$degree,KK2015_bison[[i]]$degree,
                          R2015_bison[[i]]$degree,V2015_bison[[i]]$degree, S2019_bison[[i]]$degree)
    
    networks[[i]]$top_strQ = c(F2015_bison[[i]]$top_strQ,KK2015_bison[[i]]$top_strQ,
                          R2015_bison[[i]]$top_strQ,V2015_bison[[i]]$top_strQ, S2019_bison[[i]]$top_strQ)
    
    networks[[i]]$top_str90 = c(F2015_bison[[i]]$top_str90,KK2015_bison[[i]]$top_str90,
                          R2015_bison[[i]]$top_str90,V2015_bison[[i]]$top_str90, S2019_bison[[i]]$top_str90)
    
    networks[[i]]$weak_Q = c(F2015_bison[[i]]$weak_Q,KK2015_bison[[i]]$weak_Q,
                          R2015_bison[[i]]$weak_Q,V2015_bison[[i]]$weak_Q, S2019_bison[[i]]$weak_Q)
    
    networks[[i]]$weak_10 = c(F2015_bison[[i]]$weak_10,KK2015_bison[[i]]$weak_10,
                          R2015_bison[[i]]$weak_10,V2015_bison[[i]]$weak_10, S2019_bison[[i]]$weak_10)
    
  
      #Add data from blood smears, making sure that info is added to id for respective year (some ids appear in more than one groupyear)
    
      networks[[i]] = merge(networks[[i]],wbc_behav, by = c("id","groupyear"))
    
      #Convert characters to numeric/factors
      networks[[i]]$perc_rank = as.numeric(networks[[i]]$perc_rank)
      networks[[i]]$sex = factor(networks[[i]]$sex)
      networks[[i]]$age = as.numeric(networks[[i]]$age)
      networks[[i]]$groupyear = factor(networks[[i]]$groupyear)
      networks[[i]]$id = factor(networks[[i]]$id)
      #Add a variable to separate group that suffered hurricane impact
      networks[[i]]$is_hurricane = ifelse(networks[[i]]$groupyear == "S2019", 1, 0)
      networks[[i]]$is_hurricane = factor(networks[[i]]$is_hurricane)
    }

```

###ANALYSES FOR PAPER###

```{r}
#Subset tensor to a single dataframe for rank analyses
data = networks[[1]]
data$id = as.factor(data$id)
```


#Descriptive statistics
```{r}
#total ind
length(unique(data$id))#65,

#Final sample size = 65 animals
table(data$sex)#39 F and 26 M
table(data$groupyear)#A bit spread between groups, but most of samples from F and S
range(as.numeric(data$age))#7-22


data %>% group_by(groupyear,sex) %>% summarise(n())
data %>% group_by(groupyear) %>% summarise(avg = mean(as.numeric(age)), sd = sd(as.numeric(age)))


```

#Given the small sample size, each sociality metric will be tested in separate models. 
I can not control for the effect of group membership because there are too few animals for KK, R and V 2015. 
Including a factor for hurricane might be interesting but still data is from only one group***

#First question: Does sociality influence overall health measured as NLR?

```{r}

#SOCIAL STATUS
#Log-transformed NLR: number of weak partners (< 4rd quartile)
NLR_rank = brm(data = data,
                      family = gaussian(),
                      log(NLR) ~ scale(perc_rank) + sex + scale(age) + is_hurricane,
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)
#Save model
#saveRDS(NLR_rank, file = paste0(directory,"brms_models/NLR_rank.rds"))

NLR_rank <- readRDS(paste0(directory,"brms_models/NLR_rank.rds"))#load model
pp_check(NLR_rank)#looks okay

#WEAK PARTNERS

#Log-transformed NLR: number of weak partners (< 4rd quartile)
NLR_weak = brm_multiple(data = networks,
                      family = gaussian,
                      log(NLR) ~ scale(weak_Q) + sex + scale(age) + is_hurricane,
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)

#Save model
#saveRDS(NLR_weak, file = paste0(directory,"brms_models/NLR_weak.rds"))

NLR_weak <- readRDS(paste0(directory,"brms_models/NLR_weak.rds"))#load model
pp_check(NLR_weak)#looks okay
max(NLR_weak$rhats)#Check rhat ~1:okay

#STRONG PARTNERS
NLR_str = brm_multiple(data = networks,
                      family = gaussian,
                      log(NLR) ~ scale(top_strQ) + sex + scale(age) + is_hurricane,
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)


#Save model
#saveRDS(NLR_str, file = paste0(directory,"brms_models/NLR_str.rds"))

NLR_str <- readRDS(paste0(directory,"brms_models/NLR_str.rds"))#load model
pp_check(NLR_str)#looks okay
max(NLR_str$rhats)#Check rhat ~1:okay

#Get tables with model outputs
#Social status
tab_model(NLR_rank,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "NLR",
          show.ci = 0.89)

#Weak partners
tab_model(NLR_weak,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "NLR",
          show.ci = 0.89)

#Strong partners
tab_model(NLR_str,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "NLR",
          show.ci = 0.89)

```


#Second question: Does sociality influence inflammation levels?
```{r}

#SOCIAL STATUS
neu_rank = brm(data = data,
                      family = poisson,
                      neutrophils ~ scale(perc_rank) + scale(age) + sex + is_hurricane +
                      offset(log(total+1)) + (1|id),#add observation-level effect to account for overdispersion
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 10000)

#Save model
#saveRDS(neu_rank, file = paste0(directory,"brms_models/neu_rank.rds"))

neu_rank <- readRDS(paste0(directory,"brms_models/neu_rank.rds"))#load model
pp_check(neu_rank)#looks okay


#WEAK PARTNERS
neu_weak = brm_multiple(data = networks,
                      family = poisson,
                      neutrophils ~  scale(weak_Q) + scale(age) + sex + is_hurricane +
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)


#Save model
#saveRDS(neu_weak, file = paste0(directory,"brms_models/neu_weak.rds"))

neu_weak <- readRDS(paste0(directory,"brms_models/neu_weak.rds"))#load model
pp_check(neu_weak)#looks okay
max(neu_weak$rhats)#Check rhat ~1:okay

#STRONG PARTNERS
neu_str = brm_multiple(data = networks,
                      family = poisson,
                      neutrophils ~ scale(top_strQ) + sex + scale(age) + is_hurricane +
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)

#Save model
#saveRDS(neu_str, file = paste0(directory,"brms_models/neu_str.rds"))

neu_str <- readRDS(paste0(directory,"brms_models/neu_str.rds"))#load model
pp_check(neu_str)#looks okay
max(neu_str$rhats)#Check rhat ~1:okay

#Get tables with model outputs
#Social status
tab_model(neu_rank,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Neutrophils",
          show.ci = 0.89)

#Weak partners
tab_model(neu_weak,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Neutrophils",
          show.ci = 0.89)

#Strong partners
tab_model(neu_str,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Neutrophils",
          show.ci = 0.89)

```

#Third question: Does sociality influence eosinophil levels?
```{r}

#SOCIAL STATUS
eosin_rank = brm(data = data,
                        family = poisson,
                       eosinophils ~ scale(perc_rank) + scale(age) + sex + is_hurricane +
                        offset(log(total+1)) + (1|id),
                        cores = 4,
                        prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                        iter = 10000)
  
#Save model
#saveRDS(eosin_rank, file = paste0(directory,"brms_models/eosin_rank.rds"))

eosin_rank <- readRDS(paste0(directory,"brms_models/eosin_rank.rds"))#load model
pp_check(eosin_rank)#looks okay


#WEAK PARTNERS
eosin_weak = brm_multiple(data = networks,
                      family = poisson,
                      eosinophils ~ scale(weak_Q)*sex + scale(age) + is_hurricane +
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 10000)#Had to increase iterations to reach at least 1000 for ESS

#Save model
#saveRDS(eosin_weak, file = paste0(directory,"brms_models/eosin_weak.rds"))

eosin_weak <- readRDS(paste0(directory,"brms_models/eosin_weak.rds"))#load model
pp_check(eosin_weak)#looks okay
max(eosin_weak$rhats)#chack rhats value

#Marginal effect of sex
emm = emtrends(eosin_weak, pairwise ~ sex, var = "weak_Q")
confint(emm, level = .89)
#Females: -0.045, CI = -0.087, -0.003
#Males: 0.002, CI = -0.049, 0.054
#Females vs males at average = -0.048, CI = -0.095, -0.00064

#Main effect model
eosin_weak_main = brm_multiple(data = networks,
                      family = poisson,
                      eosinophils ~ scale(weak_Q) + sex + scale(age) + is_hurricane +
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)#Had to increase iterations to reach at least 1000 for ESS

#Save model
saveRDS(eosin_weak_main, file = paste0(directory,"brms_models/eosin_weak_main.rds"))


#STRONG PARTNERS
eosin_str = brm_multiple(data = networks,
                      family = poisson,
                      eosinophils ~ scale(top_strQ)*sex + scale(age) + is_hurricane +
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)

#Save model
#saveRDS(eosin_str, file = paste0(directory,"brms_models/eosin_str.rds"))

eosin_str <- readRDS(paste0(directory,"brms_models/eosin_str.rds"))#load model
pp_check(eosin_str)#looks okay
max(eosin_str$rhats)#chack rhats value

#Marginal effect of sex
emm = emtrends(eosin_str, pairwise ~ sex, var = "top_strQ")
confint(emm, level = .89)
#Females: -0.321, CI = -1.008, 0.344
#Males: 1.04, CI = -0.0059, 2.073
#Females vs males at average = -1.37, CI = -2.44, -0.291

#Main effect model
eosin_str_main = brm_multiple(data = networks,
                      family = poisson,
                      eosinophils ~ scale(top_strQ) + sex + scale(age) + is_hurricane +
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)

#Save model
#saveRDS(eosin_str_main, file = paste0(directory,"brms_models/eosin_str_main.rds"))

eosin_str_main <- readRDS(paste0(directory,"brms_models/eosin_str_main.rds"))#load model
pp_check(eosin_str_main)#looks okay
max(eosin_str_main$rhats)#chack rhats value

#Get tables with model outputs
#Social status
tab_model(eosin_rank,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Eosinophils",
          show.ci = 0.89)

#Weak partners
tab_model(eosin_weak,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Eosinophils",
          show.ci = 0.89)
#Weak partners
tab_model(eosin_weak_main,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Eosinophils",
          show.ci = 0.89)

#Strong partners
tab_model(eosin_str,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Eosinophils",
          show.ci = 0.89)
#Strong partners
tab_model(eosin_str_main,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Eosinophils",
          show.ci = 0.89)

```

#Does sociality mediate hurricane effects on NLR?
```{r}
#Social status
NLR_rank_hurr = brm(data = data,
                      family = gaussian(),
                      log(NLR) ~ scale(perc_rank)*is_hurricane + sex + scale(age),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)
#Save model
saveRDS(NLR_rank_hurr, file = paste0(directory,"brms_models/NLR_rank_hurr.rds"))

#NLR_rank_hurr <- readRDS(paste0(directory,"brms_models/NLR_rank_hurr.rds"))#load model
#pp_check(NLR_rank_hurr)

#weak partners
NLR_weak_hurr = brm_multiple(data = networks,
                      family = gaussian(),
                      log(NLR) ~ scale(weak_Q)*is_hurricane + scale(age) + sex,
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 10000)#

#Save model
saveRDS(NLR_weak_hurr, file = paste0(directory,"brms_models/NLR_weak_hurr.rds"))

#NLR_weak_hurr <- readRDS(paste0(directory,"brms_models/NLR_weak_hurr.rds"))#load model
#pp_check(NLR_weak_hurr)#


NLR_str_hurr = brm_multiple(data = networks,
                      family = gaussian(),
                      log(NLR) ~ scale(top_strQ)*is_hurricane + sex + scale(age),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)
#Save model
saveRDS(NLR_str_hurr, file = paste0(directory,"brms_models/NLR_str_hurr.rds"))

#NLR_str_hurr <- readRDS(paste0(directory,"brms_models/NLR_str_hurr.rds"))#load model
pp_check(NLR_str_hurr)#

#Get tables with model outputs
#Social status
tab_model(NLR_rank_hurr,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "NLR",
          show.ci = 0.89)

#Weak partners
tab_model(NLR_weak_hurr,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "NLR",
          show.ci = 0.89)

#Strong partners
tab_model(NLR_str_hurr,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "NLR",
          show.ci = 0.89)

```

#Does sociality mediate hurricane effects on neutrophils?
```{r}
#Social status
neu_rank_hurr = brm(data = data,
                      family = poisson,
                      neutrophils ~ scale(perc_rank)*is_hurricane + sex + scale(age)+
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)
#Save model
saveRDS(neu_rank_hurr, file = paste0(directory,"brms_models/neu_rank_hurr.rds"))

neu_rank_hurr <- readRDS(paste0(directory,"brms_models/neu_rank_hurr.rds"))#load model
#pp_check(neu_rank_hurr)#

#weak partners
neu_weak_hurr = brm_multiple(data = networks,
                      family = poisson,
                      neutrophils ~ scale(weak_Q)*is_hurricane + scale(age) + sex +
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 10000)#

#Save model
saveRDS(neu_weak_hurr, file = paste0(directory,"brms_models/neu_weak_hurr.rds"))

neu_weak_hurr <- readRDS(paste0(directory,"brms_models/neu_weak_hurr.rds"))#load model
#pp_check(neu_weak_hurr)#


neu_str_hurr = brm_multiple(data = networks,
                      family = poisson,
                      neutrophils ~ scale(top_strQ)*is_hurricane + sex + scale(age) +
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)
#Save model
saveRDS(neu_str_hurr, file = paste0(directory,"brms_models/neu_str_hurr.rds"))

neu_str_hurr <- readRDS(paste0(directory,"brms_models/neu_str_hurr.rds"))#load model
#pp_check(neu_str_hurr)#

#Get tables with model outputs
#Social status
tab_model(neu_rank_hurr,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Neutrophils",
          show.ci = 0.89)

#Weak partners
tab_model(neu_weak_hurr,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Neutrophils",
          show.ci = 0.89)

#Strong partners
tab_model(neu_str_hurr,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Neutrophils",
          show.ci = 0.89)


```

#Does sociality mediate hurricane effects on eosinophils?
```{r}
#Social status
eosin_rank_hurr = brm(data = data,
                      family = poisson,
                      eosinophils ~ scale(perc_rank)*is_hurricane + sex + scale(age)+
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)
#Save model
saveRDS(eosin_rank_hurr, file = paste0(directory,"brms_models/eosin_rank_hurr.rds"))

eosin_rank_hurr <- readRDS(paste0(directory,"brms_models/eosin_rank_hurr.rds"))#load model
#pp_check(eosin_rank_hurr)#

#weak partners
eosin_weak_hurr = brm_multiple(data = networks,
                      family = poisson,
                      eosinophils ~ scale(weak_Q)*is_hurricane + scale(age) + sex +
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 10000)#

#Save model
saveRDS(eosin_weak_hurr, file = paste0(directory,"brms_models/eosin_weak_hurr.rds"))

eosin_weak_hurr <- readRDS(paste0(directory,"brms_models/eosin_weak_hurr.rds"))#load model
#pp_check(eosin_weak_hurr)#


eosin_str_hurr = brm_multiple(data = networks,
                      family = poisson,
                      eosinophils ~ scale(top_strQ)*is_hurricane + sex + scale(age) +
                      offset(log(total+1)) + (1|id),
                      cores = 4,
                      prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                      iter = 1000)
#Save model
saveRDS(eosin_str_hurr, file = paste0(directory,"brms_models/eosin_str_hurr.rds"))

eosin_str_hurr <- readRDS(paste0(directory,"brms_models/eosin_str_hurr.rds"))#load model
#pp_check(eosin_str_hurr)#

#Post-hoc test 3-way interactions
#weak partners
eosin_weak_int = brm_multiple(data = networks,
                          family = poisson,
                          eosinophils ~ scale(weak_Q)*sex*is_hurricane + scale(age) + 
                            offset(log(total+1)) + (1|id),
                          cores = 4,
                          prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                          iter = 10000)

#Save model
#saveRDS(eosin_weak_int, file = paste0(directory,"brms_models/eosin_weak_int.rds"))
eosin_weak_int <- readRDS(paste0(directory,"brms_models/eosin_weak_int.rds"))#load model

#Strong partners
eosin_str_int = brm_multiple(data = networks,
                         family = poisson,
                         eosinophils ~ scale(top_strQ)*sex*is_hurricane + scale(age) + 
                           offset(log(total+1)) + (1|id),
                         cores = 4,
                         prior = c(prior(student_t(5, 0, 2.5), class = "b")),
                         iter = 1000)

#Save model
#saveRDS(eosin_str_int, file = paste0(directory,"brms_models/eosin_str_int.rds"))
eosin_str_int <- readRDS(paste0(directory,"brms_models/eosin_str_int.rds"))#load model

#Get tables with model outputs
#Social status
tab_model(eosin_rank_hurr,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Eosinophils",
          show.ci = 0.89)

#Weak partners
tab_model(eosin_weak_hurr,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Eosinophils",
          show.ci = 0.89)

#Strong partners
tab_model(eosin_str_hurr,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Eosinophils",
          show.ci = 0.89)

#Triple interactions

#Weak partners
tab_model(eosin_weak_int,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Eosinophils", show.ci = 0.89)

#Strong partners
tab_model(eosin_str_int,transform =  NULL, show.r2 = F, show.p = T, show.se = T, show.stat = T, show.icc = F, show.re.var = F, dv.labels = "Eosinophils", show.ci = 0.89)


```

FIGURES

Figure 1: Bison networks for all group years
```{r}

group_year = "S2019"#Specify network
fit_bison = get(group_year)
threshold = get(paste0(group_year,"_thres"))

#Plot for bison networks
ci = 0.9;edge_scaler = 12
edgelist_bison <- get_edgelist(fit_bison, ci = ci, transform = T)
plot_edgelist <- edgelist_bison %>% 
  mutate(from = node_1, to = node_2, 
         uncertainty = round(.$'95%'-.$'5%', 5), weight = round(.$median, 5),
         scaled_weight = rescale(weight, to = c(0, 1)),
         scaled_range = rescale(uncertainty, to = c(0, 2)),
         over_threshold = if_else(median >= threshold, "yes", "no"))
plot_edge <- plot_edgelist %>%
  select(from, to, scaled_weight)
graph_edge <- graph_from_data_frame(plot_edge, directed = FALSE)
E(graph_edge)$over_threshold <- E(graph_edge)$scaled_weight >= threshold

#Set attribute for individuals sampled
names_graph =  vertex_attr(graph_edge,"name")
data2 = data
data2$id = as.character(data2$id)
data2$groupyear = as.character(data2$groupyear)
names_sampled = data2$id[data2$groupyear == group_year]
names_graph = names_graph %in% names_sampled
graph_edge <- set_vertex_attr(graph_edge, "sampled", value= names_graph)
#assign colour to the vertex based on sampled
V(graph_edge)$color = c("#2C9392","khaki")[as.numeric(as.factor(V(graph_edge)$sampled))]

plot_range <- plot_edgelist %>%
  select(from, to, scaled_range, over_threshold)
graph_range <- graph_from_data_frame(plot_range, directed = FALSE)

coords <- layout_nicely(graph_edge)
plot(graph_edge, edge.width = edge_scaler * E(graph_range)$scaled_range, 
     edge.color =  ifelse(E(graph_edge)$over_threshold, adjustcolor("grey72", 0.6), "#00000000"),  
     layout=coords, vertex.label=NA, vertex.size=5, vertex.color = adjustcolor("#2C9392", 1))
plot(graph_edge, edge.width = edge_scaler * E(graph_edge)$scaled_weight, 
     edge.color = ifelse(E(graph_edge)$over_threshold, "grey35", "#00000000"),
     vertex.label=NA, vertex.size=7,
     vertex.label.dist=0.5, vertex.label.color="#2C9392", layout=coords, add=TRUE)


```


Figure 2: Effect of the number of weak partners on eosinophil count by sex
```{r}

#Load model
eosin_weak <- readRDS(paste0(directory,"brms_models/eosin_weak.rds"))
#Extract draws from posterior
post <- eosin_weak %>%
  tidybayes::epred_draws(newdata = expand_grid(weak_Q = seq(0,37, by=5), 
                                    age = mean(data$age),
                                    is_hurricane = "0",
                                    sex = c("F","M"),
                                    total = seq(88,410, by=1)), #
              ndraws = 100,summary = FALSE,#by using false, preserve values for proportions
              re_formula = NA)#

#Compute proportions
post$prop = (post$.epred/post$total)*100

post %>%
  ggplot(aes(x = weak_Q, y = prop)) +
  ggdist::stat_lineribbon(aes(),alpha=0.2, .width = c(0.5, 0.8), linewidth = 0.7,
                   show.legend = FALSE, colour = "gray30") +
   scale_colour_manual(values = c("gray50","khaki")) +
  scale_fill_manual(values = c("khaki","grey")) +
         scale_x_continuous(breaks = seq(0,40,10)) +
         #scale_y_continuous(breaks=seq(0,1,0.2)) +
         theme_bw() +
         labs(x = "Number of weak partners", y = "Eosinophils (%)") +
        theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
        text = element_text(size = 18))+
  facet_grid(~sex)

#Posterior for males and females
#Extract posterior
post <- as_draws_df(eosin_weak) 
#extract estimates per sex
post %<>%
  mutate(b_weak_males = b_scaleweak_Q + `b_scaleweak_Q:sexM` * 1,
         b_weak_females = b_scaleweak_Q + `b_scaleweak_Q:sexM`* 0)

post %<>% select(b_weak_males,b_weak_females)

#Posterior for females
post %>% 
  ggplot(aes(x = b_weak_females, fill = after_stat(x > 0))) + ggdist::stat_halfeye(.width = c(0.66,0.89))+
  geom_vline(xintercept = 0, linetype =  "dashed", color = "gray28", size = 1)+
  scale_fill_manual(values = c("gray80", "goldenrod2")) + theme_classic(base_size = 18)+ theme(axis.line = element_blank(),
      axis.text.x=element_blank(),
      axis.text.y=element_blank(),
      axis.ticks=element_blank(),
      axis.title.x=element_blank(),
      axis.title.y=element_blank(),
      legend.position="none")

#Posterior for males
post %>% 
  ggplot(aes(x = b_weak_males, fill = after_stat(x > 0))) + ggdist::stat_halfeye(.width = c(0.66,0.89))+
  geom_vline(xintercept = 0, linetype =  "dashed", color = "gray28", size = 1)+
  scale_fill_manual(values = c("gray80", "goldenrod2")) + theme_classic(base_size = 18)+ theme(axis.line=element_blank(),
      axis.text.x=element_blank(),
      axis.text.y=element_blank(),
      axis.ticks=element_blank(),
      axis.title.x=element_blank(),
      axis.title.y=element_blank(),
      legend.position="none") 

```

Figure 3: Effect of the strength to strong partners on eosinophil count by sex
```{r}

#Load model
eosin_str <- readRDS(paste0(directory,"brms_models/eosin_str.rds"))
#Extract draws from posterior
post <- eosin_str %>%
  tidybayes::epred_draws(newdata = expand_grid(top_strQ = seq(0,1.5, by=0.1), 
                                    age = mean(data$age),
                                    is_hurricane = "0",
                                    sex = c("F","M"),
                                    total = seq(88,410, by=1)), #
              ndraws = 100,summary = FALSE,
              re_formula = NA)#
#Compute proportions
post$prop = (post$.epred/post$total)*100

post %>%
  ggplot(aes(x = top_strQ, y = prop)) +
  ggdist::stat_lineribbon(aes(),alpha=0.2, .width = c(0.5, 0.8), linewidth = 0.7,
                   show.legend = FALSE, colour = "gray30") +
   scale_colour_manual(values = c("gray50","khaki")) +
  scale_fill_manual(values = c("khaki","grey")) +
         scale_x_continuous(breaks = seq(0,1.5,0.5)) +
         scale_y_continuous(breaks=seq(0,40,10)) +
         theme_bw() +
         labs(x = "Strength to strong partners", y = "Eosinophils (%)") +
        theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
        text = element_text(size = 18))+
  facet_grid(~sex)

#Posterior for males and females
#Extract posterior
post <- as_draws_df(eosin_str) 
#extract estimates per sex
post %<>%
  mutate(b_str_males = b_scaletop_strQ + `b_scaletop_strQ:sexM` * 1,
         b_str_females = b_scaletop_strQ + `b_scaletop_strQ:sexM`* 0)

post %<>% select(b_str_males,b_str_females)

#Posterior for females
post %>% 
  ggplot(aes(x = b_str_females, fill = after_stat(x > 0))) + ggdist::stat_halfeye(.width = c(0.66,0.89))+
  geom_vline(xintercept = 0, linetype =  "dashed", color = "gray28", size = 1)+
  scale_fill_manual(values = c("gray80", "goldenrod2")) + theme_classic(base_size = 18)+ theme(axis.line=element_blank(),
      axis.text.x=element_blank(),
      axis.text.y=element_blank(),
      axis.ticks=element_blank(),
      axis.title.x=element_blank(),
      axis.title.y=element_blank(),
      legend.position="none") 

#Posterior for males
post %>% 
  ggplot(aes(x = b_str_males, fill = after_stat(x < 0))) + ggdist::stat_halfeye(.width = c(0.66,0.89))+
  geom_vline(xintercept = 0, linetype =  "dashed", color = "gray28", size = 1)+
  scale_fill_manual(values = c("gray80", "goldenrod2")) + theme_classic(base_size = 18)+ theme(axis.line=element_blank(),
      axis.text.x=element_blank(),
      axis.text.y=element_blank(),
      axis.ticks=element_blank(),
      axis.title.x=element_blank(),
      axis.title.y=element_blank(),
      legend.position="none") 
  


```

Figure 4: Results for hurricane effect for all three measures
```{r}

#Subset to data of interest
data_plot = data %>% mutate(neu_perc = neutrophils/total,
                            eosin_perc = eosinophils/total,
                            is_hurricane = ifelse(is_hurricane == 0, "Not exposed", "Exposed"),
                            is_hurricane = factor(is_hurricane, levels = c("Not exposed", "Exposed"))) %>%
                    select(is_hurricane, NLR,Neutrophils = neu_perc, Eosinophils = eosin_perc)
#transpose dataframe for multiple plots
for_plot = melt(data_plot, id = c("is_hurricane"))

ggplot(for_plot, aes(is_hurricane, value, color = is_hurricane)) + geom_boxplot(linewidth = 0.7,show.legend = FALSE) + facet_wrap(~variable, scale = "free") + theme_classic(base_size = 14) + scale_color_manual(values = c("black","darkgoldenrod"))+ labs(x = "Hurricane exposure", y = "Value")


```



SUPPLEMENTARY ANALYSES

Checking for overdispersion in count models
```{r}

#Neutrophils (without including sociality measures)
neu_poiss = brm(data = data,
                family = poisson,
                neutrophils ~ offset(log(total+1)),
                iter = 1000)


neu_poiss_OD = brm(data = data,
                   family = poisson,
                   neutrophils ~ offset(log(total+1)) + (1|id),#add observation-level effect to account for overdispersion
                   cores = 4,
                   iter = 1000)


neu_NB = brm(data = data,
             family = negbinomial(),
             neutrophils ~ offset(log(total+1)),
             cores = 4,
             iter = 1000)

#Posterior predictive check
pp_check(neu_poiss)
pp_check(neu_poiss_OD)
pp_check(neu_NB)

#Compare models using loo criteria
neu_poiss = add_criterion(neu_poiss, "loo")
neu_poiss_OD = add_criterion(neu_poiss_OD, "loo", moment_match = TRUE)
neu_NB = add_criterion(neu_NB, "loo")


loo_compare(neu_poiss,neu_poiss_OD,neu_NB, criterion = "loo")
#Best model is poisson with random effect

#Eosinophils (without social metrics)
eosin_poiss =    brm(data = data,
                      family = poisson,
                     eosinophils ~ offset(log(total+1)),
                      cores = 4,
                      iter = 1000)

eosin_poiss_OD =    brm(data = data,
                      family = poisson,
                     eosinophils ~ offset(log(total+1)) + (1|id),
                      cores = 4,
                      iter = 1000)

eosin_NB =    brm(data = data,
                      family = negbinomial(),
                     eosinophils ~ offset(log(total+1)),
                      cores = 4,
                      iter = 1000)


#Posterior predictive check
pp_check(eosin_poiss)
pp_check(eosin_poiss_OD)
pp_check(eosin_NB)

#Compare models using loo criteria
eosin_poiss = add_criterion(eosin_poiss, "loo")
eosin_poiss_OD = add_criterion(eosin_poiss_OD, "loo")
eosin_NB = add_criterion(eosin_NB, "loo")

loo_compare(eosin_poiss,eosin_poiss_OD,eosin_NB, criterion = "loo")
#Best model is the one including observation-level random effect

```


SUPPLEMENTARY FIGURES

Figure S1: Plot thresholds and edge weights from Bison 
```{r}

#Plot edge weights and threshold separately per groupyear
fit_edge = KK2015
#Define threshold for minimum empirical edge weight (different from zero)
threshold_net = KK2015_thres

#Extract 100 random draws from bison edge model
draws = 100
weights = data.frame(value = as.numeric())

for(i in 1:draws){
temp = data.frame(value = exp(fit_edge$edge_samples[i,]))
weights = rbind(weights,temp)
  }

#weights goes up to 60 decimals so need to be rounded
weights$value = round(weights$value,digits = 4)

#Define threshold for strong partners
#Quantiles
thres_Q = round(quantile(weights[weights >= threshold_net], 0.75), digits = 4)
#90%
thres_90 = round(quantile(weights[weights >= threshold_net], 0.9), digits = 4)

###Plot with ggplot
max_x = max(weights$value)

#Plot
ggplot(weights, aes(x = value)) +
  geom_density(colour = "gray48",fill= "orange2", alpha = .7)+
  xlab("Edge weights") + ylab("Density") + 
  scale_x_continuous(breaks = seq(0,max_x, 0.05)) +
  scale_y_continuous(expand = c(0.01,0)) + 
  geom_vline(xintercept = c(threshold_net, thres_Q,thres_90),
             linewidth = c(0.4),
             lty = c(2,1,1), color = c("black","red","blue"))+ 
  theme_bw()+
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=12,face="bold"), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + ggtitle("KK2015")

rm(fit_edge,temp,weights,F2015,V2015,R2015,KK2015,S2019,GROUPYEAR,only_group,thres,cayo.data.edgelist)

```



Figure S2: Correlation between sociality measures used in the study
```{r}

#Get values from 20 draws from each posterior
net_metrics = data.frame(degree = "",
                         weak_Q = "",
                         weak_10 = "",
                         top_str_Q = "",
                         top_str_90 = "",
                         rank = "")

n_draws = 50

for(k in 1:n_draws){
   post_n = networks[[k]]
   temp  = data.frame(degree = rep("", nrow(post_n)),
                      weak_Q = rep("", nrow(post_n)),
                      weak_10 = rep("", nrow(post_n)),
                      top_str_Q = rep("", nrow(post_n)),
                      top_str_90 = rep("", nrow(post_n)),
                      rank = rep("", nrow(post_n)))
   
   temp$degree = post_n$deg
   temp$weak_Q = post_n$weak_Q
   temp$weak_10 = post_n$weak_10
   temp$top_str_Q = post_n$top_strQ
   temp$top_str_90 = post_n$top_str90
   temp$rank = post_n$perc_rank

   net_metrics = rbind(net_metrics,temp)
}

net_metrics = net_metrics[-1,]
net_metrics %<>% mutate_if(is.character, as.numeric)


#Plot correlogram
cor.mat1 <- net_metrics %>% select(degree, "# weak < 4th quartile" = weak_Q, "# weak < 90%" = weak_10, "sum > 3rd quartile" = top_str_Q,
                                   "sum > 89%" = top_str_90,rank)
M <- cor(cor.mat1, method = "pearson", use = "complete.obs")
corrplot(M, type = "upper", order = "original",tl.col = "black", tl.srt = 45)

```


Figure S3: Plot proportion of strong and weak partner for same and different sex dyads
```{r}
#Get edge weights from Bison networks 
#Create dataframe where to store data
str_weights = data.frame(value = as.numeric(),
                         dyad_sex = as.character(),
                         sampled = as.character(),
                         group_year = as.character())

for(k in 1:length(group.year.list)){
    groupyear = group.year.list[k]
    #Get Bison network(change according to group year needed)
    fit_edge = get(groupyear)
    threshold_net = get(paste0(groupyear,"_thres"))
    #Get edgelist for that group year
    subset_groupyear = cayo.data.edgelist[cayo.data.edgelist$group.year == groupyear,]
    #Add a column to specify dyad sex and whether includes a sampled individual (considers at least one in the dyad sampled)
    subset_groupyear$dyad_sex = paste0(subset_groupyear$sex.id1,subset_groupyear$sex.id2)
    subset_groupyear$dyad_sex = tolower(subset_groupyear$dyad_sex)#make all lower case
    subset_groupyear$sampled = ifelse(subset_groupyear$id.1 %in% data$id, "yes", "no")
    subset_groupyear$sampled = ifelse(subset_groupyear$id.2 %in% data$id, "yes", subset_groupyear$sampled)
    
    #Extract 100 random draws from bison edge model
    draws = 100
    weights = data.frame(value = as.numeric(),
                         dyad_sex = as.character(),
                         sampled = as.character())
   
     for(i in 1:draws){
      temp = data.frame(value = exp(fit_edge$edge_samples[i,]),
                        dyad_sex = subset_groupyear$dyad_sex,
                        sampled = subset_groupyear$sampled)
      weights = rbind(weights,temp)
     }
    
    #weights goes up to 60 decimals so need to be rounded
    weights$value = round(weights$value,digits = 4)
    
    #Exclude all edges below minimum empirical edge weight
    weights = weights[weights$value >= threshold_net,] 
    #compute threshold for strong partner
    thres_str = round(quantile(weights$value, 0.75), digits = 4)
    #Create column to indicate whether tie is weak or strong
    weights$type = ifelse(weights$value >= thres_str, "Strong","Weak")
    #Assign a column for group year
    weights$group_year = groupyear
    #Append to main dataframe
    str_weights = rbind(str_weights,weights)
}
#Change labeling of dyad type
str_weights$dyad_type = ifelse(str_weights$dyad_sex == "ff", "Female-female", "none")
str_weights$dyad_type = ifelse(str_weights$dyad_sex == "fm", "Female-male", str_weights$dyad_type)
str_weights$dyad_type = ifelse(str_weights$dyad_sex == "mm", "Male-male", str_weights$dyad_type)

#Report proportions
#All dyads
table(str_weights$dyad_type,str_weights$type)
#                 Strong  Weak
#  Female-female  36336 97849
#  Female-male    27167 81399
#  Male-male       4170 23360

#Sampled dyads
table(str_weights$dyad_type[str_weights$sampled == "yes"],str_weights$type[str_weights$sampled == "yes"])
#                Strong  Weak
#  Female-female  13564 37689c
#  Female-male     9429 31417
#  Male-male       1337  7925 


#Fig S3A: Proportion of intra and intersexual dyads per relationship strength among all groups sampled
ggplot(str_weights,aes(x = type,fill = dyad_type)) + 
  geom_bar(position = "fill", width = 0.6, alpha=0.4, colour = "white") + scale_fill_manual(values = c("gray28","darkgoldenrod","gray"))+
  xlab("Strength") + ylab("Proportion among all dyads")+labs(fill = "Dyad") +
  theme_bw()+
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=12,face="bold"), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#Fig S3B: Proportion of intra and intersexual dyads per relationship strength among sampled individuals
ggplot(str_weights[str_weights$sampled == "yes",],aes(x = type,fill = dyad_type)) + 
  geom_bar(position = "fill", width = 0.6, alpha=0.4, colour = "white") + scale_fill_manual(values = c("gray28","darkgoldenrod","gray"))+
  xlab("Strength") + ylab("Proportion among sampled dyads")+labs(fill = "Dyad") +
  theme_bw()+
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=12,face="bold"), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())


```