# R script for the manuscript # "Social network position predicts male mating success in a small passerine" # Kristina B. Beck, Damien R. Farine, Bart Kempenaers # Email: kbbeck.mail@gmail.com # Load required packages require(c("lme4", "dplyr", "data.table", "ggplot2", "arm", "effects")) # Load data load("data.RData") # 1) Probability to acquire a breeding partner ###----------------------------------------------------------------------------- # Example shown for males # For female data select: data[[2]] and run: # m <- glm(breeding ~ N_male_asso + betweenness + avg_asso_male + sex_ratio + date_diff + age, data=females, family="binomial") males <- data[[1]] m <- glm(breeding ~ N_female_asso + betweenness + avg_asso_female + sex_ratio + age + date_diff, data=males, family="binomial") m2 <- standardize(m) summary(m2) # Extract coefficients for observed data degree <- m2$coefficients[2] betw <- m2$coefficients[3] avg_asso <- m2$coefficients[4] sex_ratio <- m2$coefficients[5] # Permutations males$ii <- paste(males$N_female_asso, males$betweenness, males$avg_asso_female, males$sex_ratio, sep = "_") rand_degree_epp <- rep(0,1000) rand_betw_epp <- rep(0,1000) rand_asso_epp <- rep(0,1000) sex_epp <- rep(0,1000) for (i in c(1:1000)) { # Run either restricted or unrestricted permutation males[, ii_rand:=sample(ii),] # unrestricted #males[, ii_rand:=sample(ii),by=Location] # control for spatial location rand <- data.frame(do.call('rbind', strsplit(as.character(males$ii_rand),'_',fixed=TRUE))) setnames(rand, c("X1","X2","X3","X4"), c("N_fem_rand", "betweenness_rand","avg_asso_female_rand","sex_ratio_rand")) males_r <- cbind(males, rand) males_r$N_fem_rand <- as.numeric(as.character(males_r$N_fem_rand)) males_r$betweenness_rand <- as.numeric(as.character(males_r$betweenness_rand)) males_r$avg_asso_female_rand <- as.numeric(as.character(males_r$avg_asso_female_rand)) males_r$sex_ratio_rand <- as.numeric(as.character(males_r$sex_ratio_rand)) # GLM m_r <- glm(breeding ~ N_fem_rand + betweenness_rand + avg_asso_female_rand + sex_ratio_rand + age + date_diff, data=males_r, family="binomial") m_r2 <- standardize(m_r) rand_degree_epp[i] <- m_r2$coefficients[2] rand_betw_epp[i] <- m_r2$coefficients[3] rand_asso_epp[i] <- m_r2$coefficients[4] sex_epp[i] <- m_r2$coefficients[5] } # Get P values sum(abs(degree) < abs(rand_degree_epp))/1000 sum(abs(betw) < abs(rand_betw_epp))/1000 sum(abs(avg_asso) < abs(rand_asso_epp))/1000 sum(abs(sex_ratio) < abs(sex_epp))/1000 # Figure 1A m <- glm(breeding ~ N_female_asso + age + betweenness + avg_asso_female + sex_ratio + date_diff, data=males, family="binomial") e <- allEffects(m, xlevels=50)$`betweenness` %>% data.frame %>% data.table ggplot() + geom_count(data=males, aes(y = breeding, x = betweenness), alpha=.5,shape = 1, stroke = 1.0) + scale_size(range = c(2, 10)) + geom_line(data = e, aes(y = fit, x = betweenness), size=1) + geom_ribbon(data = e, aes(x = betweenness , ymin=lower, ymax=upper, linetype=NA), alpha=.3, fill="gray46") + theme_classic() + theme(axis.text.x = element_text( color="black", size=15), axis.text.y = element_text( color="black", size=15)) + xlab("Betweenness") + ylab("Probability to breed") + ggtitle("A") + theme(axis.title.x = element_text(size=16), axis.title.y = element_text(size=16), plot.title = element_text(size=19,hjust = 0), legend.position = "none") + scale_y_continuous(limits = c(0,1)) # 2) Probability to acquire EPY ###----------------------------------------------------------------------------- # MALES # Example shown for males # For female data select: data[[4]] and run: # m <- glm(epp ~ N_male_asso + betweenness + avg_asso_male + sex_ratio, data=females, family="binomial") adult_males <- data[[3]] m <- glm(epp ~ N_female_asso + betweenness + avg_asso_female + sex_ratio, data=adult_males, family="binomial") m2 <- standardize(m) summary(m2) degree <- m2$coefficients[2] betw <- m2$coefficients[3] avg_asso <- m2$coefficients[4] sex_ratio <- m2$coefficients[5] # Permutation adult_males$ii <- paste(adult_males$N_female_asso, adult_males$betweenness, adult_males$avg_asso_female, adult_males$sex_ratio, sep = "_") rand_degree_epp <- rep(0,1000) rand_betw_epp <- rep(0,1000) rand_asso_epp <- rep(0,1000) sex_epp <- rep(0,1000) for (i in c(1:1000)) { #adult_males[, ii_rand:=sample(ii),] # unrestricted adult_males[, ii_rand:=sample(ii),by=Location] # control for spatial locationl rand <- data.frame(do.call('rbind', strsplit(as.character(adult_males$ii_rand),'_',fixed=TRUE))) setnames(rand, c("X1","X2","X3","X4"), c("N_fem_rand", "betweenness_rand","avg_asso_female_rand","sex_ratio_rand")) males_r <- cbind(adult_males, rand) males_r$N_fem_rand <- as.numeric(as.character(males_r$N_fem_rand)) males_r$betweenness_rand <- as.numeric(as.character(males_r$betweenness_rand)) males_r$avg_asso_female_rand <- as.numeric(as.character(males_r$avg_asso_female_rand)) males_r$sex_ratio_rand <- as.numeric(as.character(males_r$sex_ratio_rand)) # GLM m_r <- glm(epp ~ N_fem_rand + betweenness_rand + avg_asso_female_rand + sex_ratio_rand, data=males_r, family="binomial") m_r2 <- standardize(m_r) rand_degree_epp[i] <- m_r2$coefficients[2] rand_betw_epp[i] <- m_r2$coefficients[3] rand_asso_epp[i] <- m_r2$coefficients[4] sex_epp[i] <- m_r2$coefficients[5] } # Get P values sum(abs(degree) < abs(rand_degree_epp))/1000 sum(abs(betw) < abs(rand_betw_epp))/1000 sum(abs(avg_asso) < abs(rand_asso_epp))/1000 sum(abs(sex_ratio) < abs(sex_epp))/1000 # Figure 1B m <- glm(epp ~ N_female_asso + betweenness + avg_asso_female + sex_ratio, data=adult_males, family="binomial") e <- allEffects(m, xlevels=50)$`N_female_asso` %>% data.frame %>% data.table ggplot() + geom_count(data=adult_males, aes(y = epp, x = N_female_asso), alpha=.5,shape = 1, stroke = 1.0) + scale_size(range = c(2, 10)) + geom_line(data = e, aes(y = fit, x = N_female_asso), size=1) + geom_ribbon(data = e, aes(x = N_female_asso , ymin=lower, ymax=upper, linetype=NA), alpha=.3, fill="gray46") + theme_classic() + theme(axis.text.x = element_text( color="black", size=15), axis.text.y = element_text( color="black", size=15)) + xlab("Number of female associates") + ylab("Probability to sire EPY") + ggtitle("B") + theme(axis.title.x = element_text(size=16), axis.title.y = element_text(size=16), plot.title = element_text(size=19,hjust = 0),legend.position = "none") + scale_y_continuous(limits = c(0,1)) # 3) Does the number of female associates predict number of neighbours/proportion of familiar females? ###----------------------------------------------------------------------------- males_breeding_neigh <- data[[5]] # Number of neighbours (N) mm <- lm(N ~ rescale(N_female_asso), data=males_breeding_neigh) summary(mm) obs_m <- mm[[1:2]] # Proportion of familiar neighbours mm <- glm(cbind(N_fam_neig, N-N_fam_neig) ~ rescale(N_female_asso), data=males_breeding_neigh, family="binomial") summary(mm) obs_m2 <- mm[[1:2]] # Permutation rand_m <- rep(0,1000) rand_m2 <- rep(0,1000) for (i in c(1:1000)) { males_rand <- males_breeding_neigh males_rand[, N_female_asso_rand:= sample(N_female_asso),] # unrestricted #males_rand[, N_female_asso_rand:= sample(N_female_asso),by=Location] # control for spatial location males_rand$N_female_asso_rand <- as.numeric(males_rand$N_female_asso_rand) # model mm_r <- lm(N ~ rescale(N_female_asso_rand), data=males_rand) rand_m[i] <- mm_r[[1:2]] mm_r <- glm(cbind(N_fam_neig, N-N_fam_neig) ~ rescale(N_female_asso_rand), data=males_rand, family="binomial") rand_m2[i] <- mm_r[[1:2]] } # Get the p-values sum(abs(obs_m) < abs(rand_m))/1000 sum(abs(obs_m2) < abs(rand_m2))/1000