## This script conducts simple individual-based simulations to produce Figure A1 # #### Things you'll need #### library(foreach) library(doParallel)## neeed in order to run simulation runs in parallel library(ggplot2) library(RColorBrewer) # setup local environment for parallel computing # registerDoParallel(detectCores() - 1) getDoParWorkers() # check #### Setup #### reps = 100 # how many replicate simulations at each density? Seas.Len = 100 # Breeding season length, T in the text area.cov.m.t = 1 # Area covered per male per time step t. Assuming ideal gas dynamics, equal to 8Dv/pi, where # D = detection distance, v = speed (distance per time step) (Hutchinson & Waser 2007) N.total = 500 # total number of individuals in the population r = 1 # sex ratio (females/males) # calculate number of males and females in each simulation N.females = (r/(1+r))*N.total N.males = (1/(1+r))*N.total # Specify vector of areas to consider A = c(10000, 5000, 2500, 1250, 833, 625, 500, 416, 357, 313, 278, 250, 167, 125, 100) # Calculate corresponding densities d = N.total/A # calculate expected number of males encountered over the breeding season at each density exp.no.enc = (1/(1 + r))*N.total*area.cov.m.t*Seas.Len/A #### Best-of-all sampling (as in Eshel 1979 and Kokko & Rankin 2006) #### ## Parallel loop over simulated densities, and simulate multiple replicates within each ## Best.All.Sim.Dat <<- foreach(this.d = 1:length(d), .combine = rbind) %:% # turn multiple foreach loops into one loop foreach(this.r = 1:reps, .combine = rbind) %dopar% { # what is the current value of density? density = d[this.d] # what's the expected number of males encountered at this density? this.exp.no.enc = exp.no.enc[this.d] # make dataframes to hold all the information from the individuals in each simulation # Females need to update their mating status (Y/N) and the number of encountered males Best.Sim.F = data.frame("Sex" = c(rep("female", times = N.females)), "mated" = factor("N", levels = c("N", "Y")), "n.enc" = 0) # Males need to store their sexual trait value, and record the number of matings they achieve # each male gets an ID that's used when looking up traits of each male a female encounters Best.Sim.M = data.frame("Sex" = "male", "Sex.trait" = rnorm(n = N.males, mean = 0, sd = 1), "no.matings" = 0, "ID" = 1:N.males) ## Begin Simulation ## # for each female for (f in 1:N.females){ # draw actual number of males encountered from Poisson distribution with mean = this.exp.no.enc Best.Sim.F$n.enc[f] = rpois(1, lambda = this.exp.no.enc) # if she's encountered at least one male, she picks the most preferred from among them (largest trait) if(Best.Sim.F$n.enc[f] > 0){ Best.Sim.F$mated[f] = "Y" # update her mating status her.pool = sample(Best.Sim.M$ID, size = Best.Sim.F$n.enc[f], replace = T) # create a sample of that many males # find most preferred (largest) male in sample, store his ID this.male = her.pool[which(Best.Sim.M$Sex.trait[her.pool] == max(Best.Sim.M$Sex.trait[her.pool]))] Best.Sim.M$no.matings[this.male] = Best.Sim.M$no.matings[this.male] + 1 # update that male's no. matings } # end procedure for choosing from encountered males } # end loop over females ## to calculate Beta, need to standardize male traits to the sd (already done bc sd ~ 1) if (sum(Best.Sim.M$no.matings) > 0){ # can only do this is at least one male mated ## standardize fitness to the mean Best.Sim.M$rel.fit = Best.Sim.M$no.matings / mean(Best.Sim.M$no.matings) # calculate Beta from regression coeff Beta = lm(rel.fit ~ Sex.trait, data = Best.Sim.M)$coef[2] # Opportunity for selection, I.male.matings I.male.matings = sd(Best.Sim.M$no.matings)/mean(Best.Sim.M$no.matings) # mean matings per mated male mean.matings.per.mated.male = mean(Best.Sim.M$no.matings[which(Best.Sim.M$no.matings >= 1)]) } ## if no males had mated, then... if (sum(Best.Sim.M$no.matings) == 0){ # Beta is NA Beta = NA # Opportunity for selection, I.male.matings, is NA I.male.matings = NA # mean matings per mated male is NA (no such males exist) mean.matings.per.mated.male = NA } # now fill in other measures that don't require at least one male to have mated # mean matings per male mean.male.matings = mean(Best.Sim.M$no.matings) # sd matings per male sd.male.matings = sd(Best.Sim.M$no.matings) # proportion of males mated p.males.mated = length(which(Best.Sim.M$no.matings >= 1)) / length(Best.Sim.M$no.matings) # proportion of females mated p.females.mated = length(which(Best.Sim.F$mated == "Y")) / length(Best.Sim.F$mated) ### bind all results variables into an object that will be combined by rows across ### replicate simulations at different densities as the foreach loop runs as.data.frame(cbind(density, this.exp.no.enc, Beta, I.male.matings, mean.matings.per.mated.male, mean.male.matings, sd.male.matings, p.males.mated, p.females.mated)) } # end foreach loop over simulated densities and replicate simulations at each density # view structure of resulting object, just as a check str(Best.All.Sim.Dat) #### Best-of-n sampling (Janetos 1980, Real 1990,...) #### # specify some different values of n to try out ns = c(2, 10, 20) # note that values > 10 probably unrealistic except in cases where costs are negligible ## Parallel loop over simulated densities, and simulate multiple replicates within each ## Best.n.Sim.Dat <<- foreach(this.n = 1:length(ns), .combine = rbind) %:% # turn multiple foreach loops into one loop foreach(this.d = 1:length(d), .combine = rbind) %:% foreach(this.r = 1:reps, .combine = rbind) %dopar% { # what is the current value of n? n = ns[this.n] # what is the current value of density? density = d[this.d] # what's the expected number of males encountered at this density? this.exp.no.enc = exp.no.enc[this.d] # make dataframes to hold all the information from the individuals in each simulation # Females need to update their mating status (Y/N) and the number of encountered males Best.Sim.F = data.frame("Sex" = c(rep("female", times = N.females)), "mated" = factor("N", levels = c("N", "Y")), "n.enc" = 0) # Males need to store their sexual trait value, and record the number of matings they achieve # each male gets an ID that's used when looking up traits of each male a female encounters Best.Sim.M = data.frame("Sex" = "male", "Sex.trait" = rnorm(n = N.males, mean = 0, sd = 1), "no.matings" = 0, "ID" = 1:N.males) ## Begin Simulation ## # for each female for (f in 1:N.females){ # draw actual number of males encountered from Poisson distribution with mean = this.exp.no.enc Best.Sim.F$n.enc[f] = rpois(1, lambda = this.exp.no.enc) # if she's encountered at least n males, she picks the most preferred from among them (largest trait) if(Best.Sim.F$n.enc[f] >= n){ Best.Sim.F$mated[f] = "Y" # update her mating status her.pool = sample(Best.Sim.M$ID, size = n, replace = T) # create a sample of n males # find most preferred (largest) male in sample, store his ID this.male = her.pool[which(Best.Sim.M$Sex.trait[her.pool] == max(Best.Sim.M$Sex.trait[her.pool]))] Best.Sim.M$no.matings[this.male] = Best.Sim.M$no.matings[this.male] + 1 # update that male's no. matings } # end procedure for choosing from encountered males } # end loop over females ## to calculate Beta, need to standardize male traits to the sd (already done bc sd ~ 1) if (sum(Best.Sim.M$no.matings) > 0){ # can only do this is at least one male mated ## standardize fitness to the mean Best.Sim.M$rel.fit = Best.Sim.M$no.matings / mean(Best.Sim.M$no.matings) # calculate Beta from regression coeff Beta = lm(rel.fit ~ Sex.trait, data = Best.Sim.M)$coef[2] # Opportunity for selection, I.male.matings I.male.matings = sd(Best.Sim.M$no.matings)/mean(Best.Sim.M$no.matings) # mean matings per mated male mean.matings.per.mated.male = mean(Best.Sim.M$no.matings[which(Best.Sim.M$no.matings >= 1)]) } ## if no males had mated, then... if (sum(Best.Sim.M$no.matings) == 0){ # Beta is NA Beta = NA # Opportunity for selection, I.male.matings, is NA I.male.matings = NA # mean matings per mated male is NA (no such males exist) mean.matings.per.mated.male = NA } # now fill in other measures that don't require at least one male to have mated # mean matings per male mean.male.matings = mean(Best.Sim.M$no.matings) # sd matings per male sd.male.matings = sd(Best.Sim.M$no.matings) # proportion of males mated p.males.mated = length(which(Best.Sim.M$no.matings >= 1)) / length(Best.Sim.M$no.matings) # proportion of females mated p.females.mated = length(which(Best.Sim.F$mated == "Y")) / length(Best.Sim.F$mated) ### bind all results variables into an object that will be combined by rows across ### replicate simulations at different densities as the foreach loop runs as.data.frame(cbind(n, density, this.exp.no.enc, Beta, I.male.matings, mean.matings.per.mated.male, mean.male.matings, sd.male.matings, p.males.mated, p.females.mated)) } # end foreach loop over simulated densities and replicate simulations at each density # check object structure... str(Best.n.Sim.Dat) #### Summary statistics #### # start with "best-of-all" sampling means.Best.All = as.data.frame(aggregate(cbind(I.male.matings, Beta, p.males.mated) ~ density, Best.All.Sim.Dat, mean, na.rm = T, na.action=na.pass)) sds = as.vector(aggregate(I.male.matings~density, Best.All.Sim.Dat, sd, na.rm = T, na.action=na.pass)[2]) colnames(sds) <- c("I.mates.sd") means.Best.All$I.mates.sd = sds$I.mates.sd sds = as.vector(aggregate(Beta~density, Best.All.Sim.Dat, sd, na.rm = T, na.action=na.pass)[2]) colnames(sds) <- c("Beta.sd") means.Best.All$Beta.sd = sds$Beta.sd sds = as.vector(aggregate(p.males.mated~density, Best.All.Sim.Dat, sd, na.rm = T, na.action=na.pass)[2]) colnames(sds) <- c("p.males.mated.sd") means.Best.All$p.males.mated.sd = sds$p.males.mated.sd # now for best-of-n -- build up a dataframe with the summary stats for each n # for (this.n in 1:length(ns)){ if (this.n == 1){ means.Best.n = as.data.frame(aggregate(cbind(I.male.matings, Beta, p.males.mated) ~ density, Best.n.Sim.Dat[which(Best.n.Sim.Dat$n == ns[this.n]),], mean, na.rm = T, na.action=na.pass)) sds = as.vector(aggregate(I.male.matings~density, Best.n.Sim.Dat[which(Best.n.Sim.Dat$n == ns[this.n]),], sd, na.rm = T, na.action=na.pass)[2]) colnames(sds) <- c("I.mates.sd") means.Best.n$I.mates.sd = sds$I.mates.sd sds = as.vector(aggregate(Beta~density, Best.n.Sim.Dat[which(Best.n.Sim.Dat$n == ns[this.n]),], sd, na.rm = T, na.action=na.pass)[2]) colnames(sds) <- c("Beta.sd") means.Best.n$Beta.sd = sds$Beta.sd sds = as.vector(aggregate(p.males.mated~density, Best.n.Sim.Dat[which(Best.n.Sim.Dat$n == ns[this.n]),], sd, na.rm = T, na.action=na.pass)[2]) colnames(sds) <- c("p.males.mated.sd") means.Best.n$p.males.mated.sd = sds$p.males.mated.sd # add value of n to df means.Best.n$n = ns[this.n] } # if not on the first value of n, make a separate df and rbind it to means.Best.n if (this.n > 1){ new = as.data.frame(aggregate(cbind(I.male.matings, Beta, p.males.mated) ~ density, Best.n.Sim.Dat[which(Best.n.Sim.Dat$n == ns[this.n]),], mean, na.rm = T, na.action=na.pass)) sds = as.vector(aggregate(I.male.matings~density, Best.n.Sim.Dat[which(Best.n.Sim.Dat$n == ns[this.n]),], sd, na.rm = T, na.action=na.pass)[2]) colnames(sds) <- c("I.mates.sd") new$I.mates.sd = sds$I.mates.sd sds = as.vector(aggregate(Beta~density, Best.n.Sim.Dat[which(Best.n.Sim.Dat$n == ns[this.n]),], sd, na.rm = T, na.action=na.pass)[2]) colnames(sds) <- c("Beta.sd") new$Beta.sd = sds$Beta.sd sds = as.vector(aggregate(p.males.mated~density, Best.n.Sim.Dat[which(Best.n.Sim.Dat$n == ns[this.n]),], sd, na.rm = T, na.action=na.pass)[2]) colnames(sds) <- c("p.males.mated.sd") new$p.males.mated.sd = sds$p.males.mated.sd # add n to df new$n = ns[this.n] # now rbind to add to dataframe means.Best.n = rbind(means.Best.n, new) } } # end loop over values of n #### CREATE FIGURE A1A #### # pick colors # #display.brewer.all(n = NULL, type = "qual", select = NULL, # colorblindFriendly = TRUE) myCols = brewer.pal(n = 10, name = "Paired")[c(2,4,8,10)] myCols = c("#000000", myCols) # plot # quartz(height = 4, width = 4.5) par(mar = c(5,5,4,1)) ggplot(means.Best.All, aes(x = density, y = Beta), color = myCols[1]) + ylab(expression(paste("Standardized Selection Gradient, ", beta))) + xlab("Total Density") + geom_point(color = myCols[1]) + ylim(-0.1, 3) + geom_smooth(color = myCols[1], se = F, span = 0.7) + geom_errorbar(data = means.Best.All, aes(ymin = Beta - Beta.sd, ymax = Beta + Beta.sd), width = 0.2, color = myCols[1]) + geom_point(data = means.Best.n[which(means.Best.n$n == ns[1]),], aes(x = density, y = Beta), color = myCols[2]) + geom_smooth(data = means.Best.n[which(means.Best.n$n == ns[1]),], aes(x = density, y = Beta), color = myCols[2], se = F, span = 0.7) + geom_errorbar(data = means.Best.n[which(means.Best.n$n == ns[1]),], aes(ymin = Beta - Beta.sd, ymax = Beta + Beta.sd), width = 0.2, color = myCols[2]) + geom_point(data = means.Best.n[which(means.Best.n$n == ns[2]),], aes(x = density, y = Beta), color = myCols[3]) + geom_smooth(data = means.Best.n[which(means.Best.n$n == ns[2]),], aes(x = density, y = Beta), color = myCols[3], se = F, span = 0.7) + geom_errorbar(data = means.Best.n[which(means.Best.n$n == ns[2]),], aes(ymin = Beta - Beta.sd, ymax = Beta + Beta.sd), width = 0.2, color = myCols[3]) + geom_point(data = means.Best.n[which(means.Best.n$n == ns[3]),], aes(x = density, y = Beta), color = myCols[4]) + geom_smooth(data = means.Best.n[which(means.Best.n$n == ns[3]),], aes(x = density, y = Beta), color = myCols[4], se = F, span = 0.7) + geom_errorbar(data = means.Best.n[which(means.Best.n$n == ns[3]),], aes(ymin = Beta - Beta.sd, ymax = Beta + Beta.sd), width = 0.2, color = myCols[4]) + theme_bw() + theme(plot.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=12)) # a legend, should you need it. quartz() plot.new() legend("bottomright", legend = c("Best-of-all", "Best-of-20", "Best-of-10", "Best-of-2"), col= c(myCols[1], myCols[4:2]), # reorder those so they match lines on the plot, easier to see pch = 16) ## what is the CV of Beta over density for each strategy? sd(means.Best.All$Beta, na.rm = T) / mean(means.Best.All$Beta, na.rm = T) sd(means.Best.n$Beta[which(means.Best.n$n == ns[1])], na.rm = T) / mean(means.Best.n$Beta[which(means.Best.n$n == ns[1])], na.rm = T) sd(means.Best.n$Beta[which(means.Best.n$n == ns[2])], na.rm = T) / mean(means.Best.n$Beta[which(means.Best.n$n == ns[2])], na.rm = T) sd(means.Best.n$Beta[which(means.Best.n$n == ns[3])], na.rm = T) / mean(means.Best.n$Beta[which(means.Best.n$n == ns[3])], na.rm = T) #### Figure A1B #### # same plots, but for Opportunity for Sexual Selection quartz(height = 4, width = 4.5) par(mar = c(5,5,4,1)) ggplot(means.Best.All, aes(x = density, y = I.male.matings), color = myCols[1]) + ylab("Opportunity for Sexual Selection, I") + xlab("Total Density") + geom_point(color = myCols[1]) + ylim(-0.1, 20) + geom_smooth(color = myCols[1], se = F, span = 0.7) + geom_errorbar(data = means.Best.All, aes(ymin = I.male.matings - I.mates.sd, ymax = I.male.matings + I.mates.sd), width = 0.2, color = myCols[1]) + geom_point(data = means.Best.n[which(means.Best.n$n == ns[1]),], aes(x = density, y = I.male.matings), color = myCols[2]) + geom_smooth(data = means.Best.n[which(means.Best.n$n == ns[1]),], aes(x = density, y = I.male.matings), color = myCols[2], se = F, span = 0.7) + geom_errorbar(data = means.Best.n[which(means.Best.n$n == ns[1]),], aes(ymin = I.male.matings - I.mates.sd, ymax = I.male.matings + I.mates.sd), width = 0.2, color = myCols[2]) + geom_point(data = means.Best.n[which(means.Best.n$n == ns[2]),], aes(x = density, y = I.male.matings), color = myCols[3]) + geom_smooth(data = means.Best.n[which(means.Best.n$n == ns[2]),], aes(x = density, y = I.male.matings), color = myCols[3], se = F, span = 0.7) + geom_errorbar(data = means.Best.n[which(means.Best.n$n == ns[2]),], aes(ymin = I.male.matings - I.mates.sd, ymax = I.male.matings + I.mates.sd), width = 0.2, color = myCols[3]) + geom_point(data = means.Best.n[which(means.Best.n$n == ns[3]),], aes(x = density, y = I.male.matings), color = myCols[4]) + geom_smooth(data = means.Best.n[which(means.Best.n$n == ns[3]),], aes(x = density, y = I.male.matings), color = myCols[4], se = F, span = 0.7) + geom_errorbar(data = means.Best.n[which(means.Best.n$n == ns[3]),], aes(ymin = I.male.matings - I.mates.sd, ymax = I.male.matings + I.mates.sd), width = 0.2, color = myCols[4]) + theme_bw() + theme(plot.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=12)) # a legend, should you need it. quartz() plot.new() legend("bottomright", legend = c("Best-of-all", "Best-of-20", "Best-of-10", "Best-of-2"), col= c(myCols[1], myCols[4:2]), # reorder those so they match lines on the plot, easier to see pch = 16) #### Figure A1C #### #to help interpret changes in opportunity for sexual selection at low denisty, plot p.males.mated # quartz(height = 4, width = 4.5) par(mar = c(5,5,4,1)) ggplot(means.Best.All, aes(x = density, y = p.males.mated), color = myCols[1]) + ylab("Proportion of Males Mated") + xlab("Total Density") + geom_point(color = myCols[1]) + ylim(0, 1) + geom_smooth(color = myCols[1], se = F, span = 0.7) + geom_errorbar(data = means.Best.All, aes(ymin = p.males.mated - p.males.mated.sd, ymax = p.males.mated + p.males.mated.sd), width = 0.2, color = myCols[1]) + geom_point(data = means.Best.n[which(means.Best.n$n == ns[1]),], aes(x = density, y = p.males.mated), color = myCols[2]) + geom_smooth(data = means.Best.n[which(means.Best.n$n == ns[1]),], aes(x = density, y = p.males.mated), color = myCols[2], se = F, span = 0.7) + geom_errorbar(data = means.Best.n[which(means.Best.n$n == ns[1]),], aes(ymin = p.males.mated - p.males.mated.sd, ymax = p.males.mated + p.males.mated.sd), width = 0.2, color = myCols[2]) + geom_point(data = means.Best.n[which(means.Best.n$n == ns[2]),], aes(x = density, y = p.males.mated), color = myCols[3]) + geom_smooth(data = means.Best.n[which(means.Best.n$n == ns[2]),], aes(x = density, y = p.males.mated), color = myCols[3], se = F, span = 0.7) + geom_errorbar(data = means.Best.n[which(means.Best.n$n == ns[2]),], aes(ymin = p.males.mated - p.males.mated.sd, ymax = p.males.mated + p.males.mated.sd), width = 0.2, color = myCols[3]) + geom_point(data = means.Best.n[which(means.Best.n$n == ns[3]),], aes(x = density, y = p.males.mated), color = myCols[4]) + geom_smooth(data = means.Best.n[which(means.Best.n$n == ns[3]),], aes(x = density, y = p.males.mated), color = myCols[4], se = F, span = 0.7) + geom_errorbar(data = means.Best.n[which(means.Best.n$n == ns[3]),], aes(ymin = p.males.mated - p.males.mated.sd, ymax = p.males.mated + p.males.mated.sd), width = 0.2, color = myCols[4]) + theme_bw() + theme(plot.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=12)) # a legend, should you need it. quartz() plot.new() legend("bottomright", legend = c("Best-of-all", "Best-of-20", "Best-of-10", "Best-of-2"), col= c(myCols[1], myCols[4:2]), # reorder those so they match lines on the plot, easier to see pch = 16)