Simulation code for 'QUANTIFYING INBREEDING AVOIDANCE THROUGH EXTRA-PAIR REPRODUCTION'
Reid, Arcese, Keller, Germain, Duthie, Losdat, Wolak & Nietlisbach Evolution, 2015

This simulation code quantifies the consequences of observation failure (ie pre-banding mortality) for
estimates of the pattern of occurrence of extra-pair reproduction in relation to the kinship
between a female and her socially-paired male, and the degree of inbreeding avoidance through extra-pair reproduction.

The code is written for R, and requires packages kinship2 and lme4 and associated dependencies.

It requires five datafiles, as specified below. These files are provided as textfiles (header text will need to be deleted or bracketed out).

Some basic code to summarise the simulation output is also provided.

The data come from the long-term song sparrow field study on Mandarte Island, BC, Canada.
The data provided are sufficient to replicate the analyses presented in the above paper, and are therefore a highly restricted subset of the full Mandarte dataset.
If you are interested in running further analyses that require further data then please get in touch with
Prof Peter Arcese (University of British Columbia), Prof Lukas Keller (University of Zurich) or Prof Jane Reid (University of Aberdeen).
We are always happy to develop collaborations with researchers who have good ideas for new analyses.
We'd also appreciate it if you could let us know if you are intending to make use of the dataset in order to facilitate coordination of different ongoing research efforts
and allow us to keep track of all outputs from the long-term field study.


####Load required libraries

library(kinship2)
library(lme4)


####Read in datafile that lists all breeding attempts made during 2007-2012
breed.data <-read.table("C:\\Documents and Settings\\nhy406\\My Documents\\Song Sparrow Data\\EPP,Relatedness&Inbreeding\\EPP,Relatedness,Inbreeding2013\\Final data\\breed.data.dryad.use.txt",header=T)
nrow(breed.data)
#head(breed.data)

####Read in pedigree data
ped.data <-read.table("C:\\Documents and Settings\\nhy406\\My Documents\\Song Sparrow Data\\EPP,Relatedness&Inbreeding\\EPP,Relatedness,Inbreeding2013\\Final data\\ped.data.dryad.use.txt",header=T)
nrow(ped.data)
#head(ped.data)

####Read in list of first neighbours
neighbour1.data <-read.table("C:\\Documents and Settings\\nhy406\\My Documents\\Song Sparrow Data\\EPP,Relatedness&Inbreeding\\EPP,Relatedness,Inbreeding2013\\Final data\\neighbours1.dryad.use.txt",header=T)
nrow(neighbour1.data)
#head(neighbour1.data)

####Read in list of second neighbours
neighbour2.data <-read.table("C:\\Documents and Settings\\nhy406\\My Documents\\Song Sparrow Data\\EPP,Relatedness&Inbreeding\\EPP,Relatedness,Inbreeding2013\\Final data\\neighbours2.dryad.use.txt",header=T)
nrow(neighbour2.data)
#head(neighbour2.data)

####Read in list of non-neighbours
neighbour.na.data <-read.table("C:\\Documents and Settings\\nhy406\\My Documents\\Song Sparrow Data\\EPP,Relatedness&Inbreeding\\EPP,Relatedness,Inbreeding2013\\Final data\\neighbours.na.dryad.use.txt",header=T)
nrow(neighbour.na.data)
#head(neighbour.na.data)



###################################################################################################################################
##Calculate pairwise coefficients of kinship among adults

kin.coefs <- kinship(ped.data$indID, ped.data$sireID, ped.data$damID)
kin.coefs.df <- as.data.frame(kin.coefs)

###Function to extract kinship coefficients from the kin matrix dataframe

kinship.func <- function(DamSire,kin.frame){
		DamSire <- as.character(DamSire)
		kpair <- kin.frame[DamSire[1],DamSire[2]]
		kpair
		}

##################################################################################################################################
##Loop that assigns random sires to offspring, lets some offspring die due to inbreeding depression,
##then calculate key metrics across all offspring and across only those offspring that survived to hypothetical observation.


#####Specify number of iterations
##NB, the simulation takes 1-2 seconds per iteration
iter <- 100

#####Set up results vectors
results.mean.kdiff <- numeric(length=iter)
results.mean.kdiff.surv <- numeric(length=iter)
results.mean.kdiff.complete <- numeric(length=iter)

results.mean.kdiff.diff <- numeric(length=iter)
results.mean.kdiff.diff.complete <- numeric(length=iter)

results.coef.all <- numeric(length=iter)
results.coef.surv <- numeric(length=iter)
results.coef.diff <- numeric(length=iter)
results.coef.complete.clutch <- numeric(length=iter)
results.coef.complete.diff <- numeric(length=iter)

results.prop.band <- numeric(length=iter)

results.leth.equiv.reg <- numeric(length=iter)
results.leth.equiv.reg.int <- numeric(length=iter)
results.leth.equiv.reg.corr <- numeric(length=iter)

results.ID <- numeric(length=iter)
results.max <- numeric(length=iter)

results.mean.f.off.die <- numeric(length=iter)
results.mean.f.off.surv <- numeric(length=iter)
results.mean.f.off.diff <- numeric(length=iter)

results.leth.equiv.reg.soc <- numeric(length=iter)
results.leth.equiv.reg.int.soc <- numeric(length=iter)


######Start main loop
for (i in 1:iter) {


########Randomly select one male from each candidate neighbour set per nest record (nestrec)
##This gives the potential extra-pair sires for each egg

####First neighbours of each female
neighbour.sel1 <- by(neighbour1.data, neighbour1.data$nestrec, function(x){
    samp<- sample(as.character(unique(x$neighbour1ID)), 1)
    x[as.character(x$neighbour1ID) == samp, ]})

fems.neighbour1.single <- as.data.frame(do.call(rbind, neighbour.sel1))
fems.neighbour1.single <- unique(subset(fems.neighbour1.single, select = c("nestrec", "year", "socsirenID", "gendamnID", "neighbour1ID")))
nrow(fems.neighbour1.single)


####Second neighbours of each female
neighbour.sel2 <- by(neighbour2.data, neighbour2.data$nestrec, function(x){
    samp<- sample(as.character(unique(x$neighbour2ID)), 1)
    x[as.character(x$neighbour2ID) == samp, ]})

fems.neighbour2.single <- as.data.frame(do.call(rbind, neighbour.sel2))
fems.neighbour2.single <- unique(subset(fems.neighbour2.single, select = c("nestrec", "year", "socsirenID", "gendamnID", "neighbour2ID")))
nrow(fems.neighbour2.single)


####Non-neighbours of each female
neighbour.sel.na <- by(neighbour.na.data, neighbour.na.data$nestrec, function(x){
    samp<- sample(as.character(unique(x$neighbour.naID)), 1)
    x[as.character(x$neighbour.naID) == samp, ]})

fems.neighbour.na.single <- as.data.frame(do.call(rbind, neighbour.sel.na))
fems.neighbour.na.single <- unique(subset(fems.neighbour.na.single, select = c("nestrec", "year", "socsirenID", "gendamnID", "neighbour.naID")))
nrow(fems.neighbour.na.single)


##########Select a single random extra-pair sire for each nestrec
##according to some specified probability of an extra-pair sire being a first neighbour,
##a second neighbour or a non-neighbour
##Set these probabilities equal to the observed frequencies in the real song sparrow data

##Create a dataset that lists the three randomly selected males for each nestrec (ie, a first neighbour, a second neighbour and a non-neighbour)
fems.neighbour.single.all.data <- merge(fems.neighbour1.single, fems.neighbour2.single)
fems.neighbour.single.all.data <- merge(fems.neighbour.single.all.data, fems.neighbour.na.single)

##Create random number between 0 and 1 specifiying the neighbour category to select a single random EP male from
fems.neighbour.single.all.data$rand <- runif(nrow(fems.neighbour.single.all.data), 0, 1)

##NB, default thresholds are 0.11 and 0.025
##Set to 0.02 and 0.01 or similar for tighter focus on first-neighbour males
##Set to 0.98 and 0.95 or similar for tighter focus on non-neighbour males

fems.neighbour.single.all.data$rand.male <- ifelse(fems.neighbour.single.all.data$rand > 0.11, fems.neighbour.single.all.data$neighbour1ID,
					ifelse(fems.neighbour.single.all.data$rand < 0.025, fems.neighbour.single.all.data$neighbour.naID,
					fems.neighbour.single.all.data$neighbour2ID))

rand.male.data <- data.frame(fems.neighbour.single.all.data$nestrec, fems.neighbour.single.all.data$rand.male)
colnames(rand.male.data) <- c("nestrec", "rand.male")
#head(rand.male.data)
nrow(rand.male.data)


##Merge randomly selected male to each observed nestrec
##NB, this assigns a single random extra-pair male to all extra-pair offspring from the same brood (ie nestrec)

breed.data.cut1.rand.sire <- merge(breed.data, rand.male.data)
#head(breed.data.cut1.rand.sire)
nrow(breed.data.cut1.rand.sire)


##############Expand dataset to create one row per egg laid

breed.data.exp <- data.frame(rep(breed.data.cut1.rand.sire$nestrec, times = breed.data.cut1.rand.sire$clutch), rep(breed.data.cut1.rand.sire$year, times = breed.data.cut1.rand.sire$clutch),
rep(breed.data.cut1.rand.sire$socsireID, times = breed.data.cut1.rand.sire$clutch), rep(breed.data.cut1.rand.sire$gendamID, times = breed.data.cut1.rand.sire$clutch),
rep(breed.data.cut1.rand.sire$clutch, times = breed.data.cut1.rand.sire$clutch), rep(breed.data.cut1.rand.sire$banded, times = breed.data.cut1.rand.sire$clutch),
rep(breed.data.cut1.rand.sire$rand.male, times = breed.data.cut1.rand.sire$clutch))
colnames(breed.data.exp) <- colnames(breed.data.cut1.rand.sire)
#head(breed.data.exp)
nrow(breed.data.exp)


##############Assign a genetic sire to each egg
##Create random variable to choose within-pair or extra-pair sire
breed.data.exp$rand <- runif(nrow(breed.data.exp), 0, 1)

##Assign sire
##Assign socsire if rand < 0.76 and rand.male otherwise
##Change this value to change the assumed extra-pair paternity rate, which is 0.24 across the real song sparrow dataset

breed.data.exp$sire <- ifelse(breed.data.exp$rand < 0.76, breed.data.exp$socsireID, breed.data.exp$rand.male)
#head(breed.data.exp)
nrow(breed.data.exp)

##Define each offspring as extra-pair or within-pair
breed.data.exp$epy <- ifelse(breed.data.exp$rand < 0.76, 0, 1)
#head(breed.data.exp)


##############Calculate kinship between each female and her social and genetic mate

####Create DamSire object for female and her paired social mate
DamSocSire.all <- cbind(as.numeric(as.character(breed.data.exp$gendamID)),as.numeric(as.character(breed.data.exp$socsireID)))
nrow(DamSocSire.all)

k.dam.socsire.all <- apply(DamSocSire.all, 1, kinship.func, kin.frame = kin.coefs.df)
length(k.dam.socsire.all)

##Add to main dataframe
breed.data.exp$k.dam.socsire <- k.dam.socsire.all


####Create DamSire object for female and the sire of each offspring
DamGenSire.all <- cbind(as.numeric(as.character(breed.data.exp$gendamID)),as.numeric(as.character(breed.data.exp$sire)))
nrow(DamGenSire.all)

k.dam.gensire.all <- apply(DamGenSire.all, 1, kinship.func, kin.frame = kin.coefs.df)
length(k.dam.gensire.all)

##Add to main dataframe
breed.data.exp$k.dam.gensire <- k.dam.gensire.all


####Calculate kdiff - the difference in kinship through extra-pair reproduction for each egg
breed.data.exp$kdiff <- breed.data.exp$k.dam.gensire- breed.data.exp$k.dam.socsire
#head(breed.data.exp)

##Count the total numbers of extra-pair and within-pair offspring assigned
table(breed.data.exp$epy)


############Calculate kdiff across all extra-pair offspring (ie, all extra-pair eggs)
####Reduce dataset to extra-pair offspring
breed.data.epy <- subset(breed.data.exp, breed.data.exp$epy == 1)
nrow(breed.data.epy)

##Summarise kdiff across extra-pair offspring
summary(breed.data.epy$kdiff)


###########Count the numbers of extra-pair and within-pair offspring produced within each clutch
#NB, for each clutch, this also provides original clutch size and the female's kinship with her socially paired male

clutch.off.all <- aggregate(breed.data.exp$epy, list(breed.data.exp$nestrec, breed.data.exp$year,
breed.data.exp$socsireID, breed.data.exp$gendamID, breed.data.exp$clutch, breed.data.exp$k.dam.socsire), sum)
colnames(clutch.off.all) <- c("nestrec", "year", "socsireID", "gendamID", "clutch", "k.dam.socsire", "clutch.epy.all")

clutch.off.all$clutch.wpy.all <- clutch.off.all$clutch - clutch.off.all$clutch.epy.all
#head(clutch.off.all)
nrow(clutch.off.all)


##Test whether probability of producing an extra-pair offspring varies with pair relatedness across all pairs and offspring.
##It should not, because paternity has been randomly assigned
##NB, can also fit a mixed model to account for repeat observations of females and broods.

model.epr.all <- glm(cbind(clutch.epy.all, clutch.wpy.all) ~ k.dam.socsire,
family = binomial, data = clutch.off.all)
#summary(model.epr.all)


##############Make all eggs die or survive to banding according to how inbred they are

###########log-linear decrease 
##NB The inbreeding load that's created depends on the parameters max and ID (where ID is the log-linear slope)
##Hence a function is written below to extract the realised degree of inbreeding depression
##measured as (haploid) lethal equivalents
##Keep min = 0
min <- 0
max <- runif(1, 0.25, 2)
ID <- runif(1, 0, 3)

##Calculate individual probability of egg survival as a function of coefficient of inbreeding (f)
breed.data.exp$surv.rand <- runif(nrow(breed.data.exp), min, max)

breed.data.exp$surv.rand.f <- exp(-ID*breed.data.exp$k.dam.gensire - breed.data.exp$surv.rand)

plot(breed.data.exp$k.dam.gensire, breed.data.exp$surv.rand.f,
xlab = "Egg's coefficient of inbreeding", ylab = "Egg's probability of surviving to banding")


#######Assign each offspring as survived or dead
##Specify proportion of offspring that die before banding
##Real observed values are: all years 0.17
##2007 - 0.09, 2008 - 0.24, 2009 - 0.15, 2010 - 0.17, 2011 - 0.18, 2012 - 0.18

##This code allows for a year-specific mortality rate
##where prop.dead applies to 2009-2012 (which were similar in the real dataset)
#Set all three values to 0.17 for a constant rate
prop.dead.2007 <- 0.17
prop.dead.2008 <- 0.17
prop.dead <- 0.17

q12.2007 <- quantile(breed.data.exp$surv.rand.f, prop.dead.2007)

q12.2008 <- quantile(breed.data.exp$surv.rand.f, prop.dead.2008)

q12 <- quantile(breed.data.exp$surv.rand.f, prop.dead)

breed.data.exp$surv <- ifelse(breed.data.exp$year == 2007 & breed.data.exp$surv.rand.f > q12.2007, 1, 
		ifelse(breed.data.exp$year == 2008 & breed.data.exp$surv.rand.f > q12.2008, 1,
		ifelse(breed.data.exp$year > 2008 & breed.data.exp$surv.rand.f > q12, 1, 0)))

#table(breed.data.exp$surv)
#table(breed.data.exp$surv, breed.data.exp$year)


####Set up dataset to calculate realised haploid lethal equivalents for survival based on regression across pooled categories
##This is across all offspring with no observation failure

##Group offspring into categories in f that are equal in terms of numbers of individuals

quants <- quantile(breed.data.exp$k.dam.gensire, probs=seq(0, 1, 0.1))
breed.data.exp$fcat <- ifelse(breed.data.exp$k.dam.gensire < quants[2], 1,
			ifelse(breed.data.exp$k.dam.gensire >= quants[2] & breed.data.exp$k.dam.gensire < quants[3], 2,
			ifelse(breed.data.exp$k.dam.gensire >= quants[3] & breed.data.exp$k.dam.gensire < quants[4], 3,
			ifelse(breed.data.exp$k.dam.gensire >= quants[4] & breed.data.exp$k.dam.gensire < quants[5], 4,
			ifelse(breed.data.exp$k.dam.gensire >= quants[5] & breed.data.exp$k.dam.gensire < quants[6], 5,
			ifelse(breed.data.exp$k.dam.gensire >= quants[6] & breed.data.exp$k.dam.gensire < quants[7], 6,
			ifelse(breed.data.exp$k.dam.gensire >= quants[7] & breed.data.exp$k.dam.gensire < quants[8], 7,
			ifelse(breed.data.exp$k.dam.gensire >= quants[8] & breed.data.exp$k.dam.gensire < quants[9], 8,
			ifelse(breed.data.exp$k.dam.gensire >= quants[9] & breed.data.exp$k.dam.gensire < quants[10], 9, 10)))))))))


#Calculate numbers of survivors and non-survivors and mean f within each category of f

fcat.data <- data.frame(aggregate(breed.data.exp$surv, by = list(breed.data.exp$fcat), length))
colnames(fcat.data) <- c("fcat", "total.n")

fcat.data.surv <- data.frame(aggregate(breed.data.exp$surv, by = list(breed.data.exp$fcat), sum))
colnames(fcat.data.surv) <- c("fcat", "surv.n")

fcat.data.f <- data.frame(aggregate(breed.data.exp$k.dam.gensire, by = list(breed.data.exp$fcat), mean))
colnames(fcat.data.f) <- c("fcat", "mean.f")

#Merge to create dataset for analysis
fcat.data.all <- merge(fcat.data, fcat.data.surv)
fcat.data.all <- merge(fcat.data.all, fcat.data.f)

##Calculate raw proportion survivors
fcat.data.all$prop.surv <- fcat.data.all$surv.n/fcat.data.all$total.n

##Apply small noise approximation to deal with the rare cases where the proportion of survivors is zero
fcat.data.all$prop.surv.corr <- (1 + fcat.data.all$surv.n)/(2 + fcat.data.all$total.n)

fcat.data.all$prop.surv.use <- ifelse(fcat.data.all$prop.surv > 0, fcat.data.all$prop.surv, fcat.data.all$prop.surv.corr)
#fcat.data.all


##Calculate log proportion survivors
fcat.data.all$prop.surv.log <- log(fcat.data.all$prop.surv.use)
fcat.data.all$prop.surv.corr.log <- log(fcat.data.all$prop.surv.corr)
fcat.data.all

##Regression of log proportion survivors on mean f across categories to estimate lethal equivalents
model.ID.reg <- lm(fcat.data.all$prop.surv.log ~ fcat.data.all$mean.f)

##Slope
leth.equiv.reg <- -coef(model.ID.reg)[2]
#leth.equiv.reg

##Intercept
leth.equiv.reg.int <- -coef(model.ID.reg)[1]
#leth.equiv.reg.int


##Regression of log proportion survivors on mean f across categories to estimate lethal equivalents
##using full small noise approximation
model.ID.reg.corr <- lm(fcat.data.all$prop.surv.corr.log ~ fcat.data.all$mean.f)

leth.equiv.reg.corr <- -coef(model.ID.reg.corr)[2]
#leth.equiv.reg.corr


###############Select extra-pair offspring that survived to banding and look at distribution of kdiff
breed.data.epy.surv <- subset(breed.data.exp, breed.data.exp$surv == 1 & breed.data.exp$epy == 1)
nrow(breed.data.epy.surv)
#summary(breed.data.epy.surv$kdiff)


###############Count the numbers of extra-pair and within-pair offspring that survived from each clutch
##Reduce the dataset to only individuals that survived to banding
##NB, tend to lose some clutches here, where all simulated offspring died

breed.data.exp.surv <- subset(breed.data.exp, breed.data.exp$surv == 1)
nrow(breed.data.exp.surv)
#head(breed.data.exp.surv)

##Total offspring that survived per clutch
clutch.off.surv <- aggregate(breed.data.exp.surv$nestrec, list(breed.data.exp.surv$nestrec, breed.data.exp.surv$year,
breed.data.exp.surv$socsireID, breed.data.exp.surv$gendamID, breed.data.exp.surv$clutch, breed.data.exp.surv$k.dam.socsire), length)
colnames(clutch.off.surv) <- c("nestrec", "year", "socsireID", "gendamID", "clutch", "k.dam.socsire", "clutch.off.surv")
#head(clutch.off.surv)
nrow(clutch.off.surv)

##Total extra-pair offspring that survived per clutch
clutch.epy.surv <- aggregate(breed.data.exp.surv$epy, list(breed.data.exp.surv$nestrec), sum)
colnames(clutch.epy.surv) <- c("nestrec", "clutch.epy.surv")
#head(clutch.epy.surv)
nrow(clutch.epy.surv)

##Create full dataset and count within-pair offspring that survived per clutch
clutch.off.surv.use <- merge(clutch.off.surv, clutch.epy.surv)
clutch.off.surv.use$clutch.wpy.surv <- clutch.off.surv.use$clutch.off.surv - clutch.off.surv.use$clutch.epy.surv
#head(clutch.off.surv.use)
nrow(clutch.off.surv.use)


##Test whether probability of producing an observed extra-pair offspring varies with pair relatedness across observed offspring from all broods
model.epr.surv <- glm(cbind(clutch.epy.surv, clutch.wpy.surv) ~ k.dam.socsire,
family = binomial, data = clutch.off.surv.use)


##Test whether probability of producing an observed extra-pair offspring varies with pair relatedness across observed offspring
##from broods where all offspring were observed
#head(clutch.off.surv.use)

clutch.off.surv.complete <- subset(clutch.off.surv.use, clutch.off.surv.use$clutch - clutch.off.surv.use$clutch.off.surv == 0)
nrow(clutch.off.surv.complete)

model.epr.surv.complete <- glm(cbind(clutch.epy.surv, clutch.wpy.surv) ~ k.dam.socsire,
family = binomial, data = clutch.off.surv.complete)
summary(model.epr.surv.complete)


##Calculate kdiff across extra-pair offspring from broods where all offspring were observed

##Get list of nest records where all offspring were observed
complete.nestrecs <- data.frame(clutch.off.surv.complete$nestrec, rep(1, times = length(clutch.off.surv.complete$nestrec)))
colnames(complete.nestrecs) <- c("nestrec", "complete")
#head(complete.nestrecs)
nrow(complete.nestrecs)

##Merge back to main dataset for observed extra-pair offspring
breed.data.epy.surv.complete <- merge(breed.data.epy.surv, complete.nestrecs)
#head(breed.data.epy.surv.complete)
nrow(breed.data.epy.surv.complete)

#summary(breed.data.epy.surv.complete$kdiff)

#####
##Create dataset that allows calculation of the probability that an egg would survive to hypothetical observation in relation to
##the kinship between a female and her socially paired male done per clutch
##This is to allow estimation of the degree of bias in inbreeding depression in early survival estimated from
##observed hatching success across social pairings (see Supporting Information)

clutch.off.surv.die <- aggregate(breed.data.exp$surv, list(breed.data.exp$nestrec, breed.data.exp$year,
breed.data.exp$socsireID, breed.data.exp$gendamID, breed.data.exp$clutch, breed.data.exp$k.dam.socsire), sum)
colnames(clutch.off.surv.die) <- c("nestrec", "year", "socsireID", "gendamID", "clutch", "k.dam.socsire", "clutch.off.surv")
clutch.off.surv.die$clutch.off.die <- clutch.off.surv.die$clutch - clutch.off.surv.die$clutch.off.surv
#head(clutch.off.surv.die)
#nrow(clutch.off.surv.die)

##Fit model that predicts the probability that an offspring will hatch in relation to k between the female and its socially paired male
##(as done for the real song sparrow data)

model.dead.mixed.bin.est <- glmer(cbind(clutch.off.surv, clutch.off.die) ~ k.dam.socsire + as.factor(year) + (1 | gendamID),
family = binomial, data = clutch.off.surv.die)
#summary(model.dead.mixed.bin.est)


##Use regression method to estimate lethal equivalents in relation to k.dam.socsire across all offspring
####Set up dataset to calculate realised haploid lethal equivalents for survival based on regression across pooled categories
##based on k.dam.socsire between a female and her observed socially paired male

quants.soc <- quantile(breed.data.exp$k.dam.socsire, probs=seq(0, 1, 0.1))
breed.data.exp$fcat.soc <- ifelse(breed.data.exp$k.dam.socsire < quants.soc[2], 1,
			ifelse(breed.data.exp$k.dam.socsire >= quants.soc[2] & breed.data.exp$k.dam.socsire < quants.soc[3], 2,
			ifelse(breed.data.exp$k.dam.socsire >= quants.soc[3] & breed.data.exp$k.dam.socsire < quants.soc[4], 3,
			ifelse(breed.data.exp$k.dam.socsire >= quants.soc[4] & breed.data.exp$k.dam.socsire < quants.soc[5], 4,
			ifelse(breed.data.exp$k.dam.socsire >= quants.soc[5] & breed.data.exp$k.dam.socsire < quants.soc[6], 5,
			ifelse(breed.data.exp$k.dam.socsire >= quants.soc[6] & breed.data.exp$k.dam.socsire < quants.soc[7], 6,
			ifelse(breed.data.exp$k.dam.socsire >= quants.soc[7] & breed.data.exp$k.dam.socsire < quants.soc[8], 7,
			ifelse(breed.data.exp$k.dam.socsire >= quants.soc[8] & breed.data.exp$k.dam.socsire < quants.soc[9], 8,
			ifelse(breed.data.exp$k.dam.socsire >= quants.soc[9] & breed.data.exp$k.dam.socsire < quants.soc[10], 9, 10)))))))))


#Calculate numbers of survivors and non-survivors and mean social f within each category of social f

fcat.soc.data <- data.frame(aggregate(breed.data.exp$surv, by = list(breed.data.exp$fcat.soc), length))
colnames(fcat.soc.data) <- c("fcat.soc", "total.n.soc")

fcat.soc.data.surv <- data.frame(aggregate(breed.data.exp$surv, by = list(breed.data.exp$fcat.soc), sum))
colnames(fcat.soc.data.surv) <- c("fcat.soc", "surv.n.soc")

fcat.soc.data.f <- data.frame(aggregate(breed.data.exp$k.dam.socsire, by = list(breed.data.exp$fcat.soc), mean))
colnames(fcat.soc.data.f) <- c("fcat.soc", "mean.f.soc")

#Merge to create dataset for analysis
fcat.soc.data.all <- merge(fcat.soc.data, fcat.soc.data.surv)
fcat.soc.data.all <- merge(fcat.soc.data.all, fcat.soc.data.f)

##Calculate raw proportion survivors
fcat.soc.data.all$prop.surv.soc <- fcat.soc.data.all$surv.n.soc/fcat.soc.data.all$total.n.soc

##Apply small noise approximation to deal with the rare cases where the proportion of survivors is zero
fcat.soc.data.all$prop.surv.soc.corr <- (1 + fcat.soc.data.all$surv.n.soc)/(2 + fcat.soc.data.all$total.n.soc)

fcat.soc.data.all$prop.surv.soc.use <- ifelse(fcat.soc.data.all$prop.surv.soc > 0, fcat.soc.data.all$prop.surv.soc, fcat.soc.data.all$prop.surv.soc.corr)


##Calculate log proportion survivors
fcat.soc.data.all$prop.surv.soc.log <- log(fcat.soc.data.all$prop.surv.soc.use)
fcat.soc.data.all$prop.surv.soc.corr.log <- log(fcat.soc.data.all$prop.surv.soc.corr)
fcat.soc.data.all

##Regression of log proportion survivors on mean f across categories to estimate lethal equivalents
model.ID.reg.soc <- lm(fcat.soc.data.all$prop.surv.soc.log ~ fcat.soc.data.all$mean.f.soc)

leth.equiv.reg.soc <- -coef(model.ID.reg.soc)[2]
leth.equiv.reg.soc

leth.equiv.reg.int.soc <- -coef(model.ID.reg.soc)[1]
leth.equiv.reg.int.soc



#####Collect required statistics
####Difference in kinship between a female and her socially-paired versus extra-pair male
##All eggs
mean.kdiff <- mean(breed.data.epy$kdiff)
##Survived eggs
mean.kdiff.surv <- mean(breed.data.epy.surv$kdiff)
##Eggs from completely observed clutches
mean.kdiff.complete <- mean(breed.data.epy.surv.complete$kdiff)


####Regression slopes of extra-pair reproduction on the kinship between a female and her socially paired male
##All eggs
coef.all <- coef(summary(model.epr.all))[2]
##Survived eggs
coef.surv <- coef(summary(model.epr.surv))[2]
##Eggs from completely observed clutches
coef.complete.clutch <- coef(summary(model.epr.surv.complete))[2]

#Proportion of offspring that died before banding
prop.band <- table(breed.data.exp$surv)[1]/sum(table(breed.data.exp$surv))

##Mean f of offspring that died and survived
mean.f.off.die <- tapply(breed.data.exp$k.dam.gensire, breed.data.exp$epy, mean)[1]
mean.f.off.surv <- tapply(breed.data.exp$k.dam.gensire, breed.data.exp$epy, mean)[2]


#######Fill in results vectors
#Difference in kinship between a female and her socially-paired versus extra-pair male
#All eggs
results.mean.kdiff[i] <- mean.kdiff
#Survived eggs
results.mean.kdiff.surv[i] <- mean.kdiff.surv
#Eggs from completely observed clutches
results.mean.kdiff.complete[i] <- mean.kdiff.complete

#Bias in kdiff due to failure to observe eggs that died before banding
results.mean.kdiff.diff[i] <- mean.kdiff.surv - mean.kdiff
#Bias when analysis is restricted to fully observed clutches
results.mean.kdiff.diff.complete[i] <- mean.kdiff.complete - mean.kdiff

#Regression of the probability of extra-pair reproduction on a female's kinship with her socially paired male
#All eggs
results.coef.all[i] <- coef.all
#Survived eggs
results.coef.surv[i] <- coef.surv
#Eggs from completely observed clutches
results.coef.complete.clutch[i] <- coef.complete.clutch

#Bias in regression due to failure to observe eggs that died before banding
results.coef.diff[i] <- coef.surv - coef.all
#Bias when analysis is restricted to fully observed clutches
results.coef.complete.diff[i] <- coef.complete.clutch - coef.all

#Proportion of eggs that survived to hypothetical banding
results.prop.band[i] <- prop.band

#Estimate of lethal equivalents in egg survival to banding
results.leth.equiv.reg[i] <- leth.equiv.reg
results.leth.equiv.reg.int[i] <- leth.equiv.reg.int
results.leth.equiv.reg.corr[i] <- leth.equiv.reg.corr

##Values of the parameters max and ID that control the magnitude of inbreeding depression
results.ID[i] <- ID
results.max[i] <- max

##Inbreeding coefficients of eggs that died and survived
results.mean.f.off.die[i] <- mean.f.off.die
results.mean.f.off.surv[i] <- mean.f.off.surv
results.mean.f.off.diff[i] <- mean.f.off.die - mean.f.off.surv

results.leth.equiv.reg.soc[i] <- leth.equiv.reg.soc
results.leth.equiv.reg.int.soc[i] <- leth.equiv.reg.int.soc

}



#####Some basic data description and summarising

####Mean kdiff
##Mean kdiff across all offspring
#results.mean.kdiff
summary(results.mean.kdiff)
sd(results.mean.kdiff)
hist(results.mean.kdiff)

##Mean kdiff across offspring that survived to observation
#results.mean.kdiff.surv
summary(results.mean.kdiff.surv)
hist(results.mean.kdiff.surv)

##Bias in mean kdiff caused by observation failure
#results.mean.kdiff.diff
summary(results.mean.kdiff.diff)
hist(results.mean.kdiff.diff)

##Mean kdiff estimated across extra-pair offspring from completely observed broods
summary(results.mean.kdiff.complete)
hist(results.mean.kdiff.complete)

##Bias in mean kdiff estimated across EPO from completely observed broods
summary(results.mean.kdiff.diff.complete)
hist(results.mean.kdiff.diff.complete)


####Regression of probability of extra-pair reproduction on a female's kinship with her socially paired male
##Across all simulated eggs
#results.coef.all
summary(results.coef.all)
sd(results.coef.all)
hist(results.coef.all)

##Across eggs that survived to observation
#results.coef.surv
summary(results.coef.surv)
hist(results.coef.surv)

##Bias due to observation failure
#results.coef.diff
hist(results.coef.diff)
summary(results.coef.diff)


####Regression of probability of extra-pair reproduction on kinship between a female and her socially paired male
##across completely observed broods only
#results.coef.complete.clutch
summary(results.coef.complete.clutch)
sd(results.coef.complete.clutch)
hist(results.coef.complete.clutch)

##Bias due to selecting only completely observed broods
summary(results.coef.complete.diff)
hist(results.coef.complete.diff)


####Proportion of offspring that died before banding
#results.prop.band
summary(results.prop.band)

####Lethal equivalents
#results.leth.equiv.reg
summary(results.leth.equiv.reg)
hist(results.leth.equiv.reg)


##Compare true lethal equivalents to those estimated from social pairings
summary(results.leth.equiv.reg)
summary(results.leth.equiv.reg.soc)

