## Code for carrying out quantitative genetic analysis of fitness data

## MALES ##

destfile.m <- "raw_male_fitness_data.csv"
url.m <- "https://zenodo.org/record/571168/files/raw_male_fitness_data.csv"
download.file(url.m, destfile.m)

# load in non-standardised male fitness data
male <- read.table("raw_male_fitness_data.csv", header = T, sep = ",")
male <- na.omit(male) # remove NA rows - vials with no offspring
# calculate proportion red-eyes of total
male$prop <- male$total_red / (male$total_red + male$total_brown)
hist(male$prop)
# calculate proportion per male
male$prop.per.m <- male$prop / male$male_density
hist(male$prop.per.m)
# calculate relative fecundity by dividing proportion by median proportion
m.median <- median(male$prop.per.m)
male$rel.fec <- male$prop.per.m / m.median
hist(male$rel.fec)
# calculate standardized relative fecundity
male$rel.st <- (male$rel.fec-mean(male$rel.fec))/sd(male$rel.fec)
hist(male$rel.st)


## FEMALES ##

destfile.f <- "raw_female_fitness_data.csv"
url.f <- "https://zenodo.org/record/571168/files/raw_female_fitness_data.csv"
download.file(url.f, destfile.f)

# load in non-standardised female fitness data
female <- read.table("raw_female_fitness_data.csv", header=T, sep = ",")
# function to calculate mean number of eggs per female, omitting NA or 0/1 eggs
female$mean.eggs <- apply(female[ , c(4:8) > 1], 1, mean, na.rm = T)
hist(female$mean.eggs)
# calculate relative fecundity by dividing eggs by median eggs
f.median <- median(female$mean.eggs)
female$rel.fec <- female$mean.eggs / f.median
hist(female$rel.fec)
# calculate standardized relative fecundity
female$rel.st <- (female$rel.fec-mean(female$rel.fec))/sd(female$rel.fec)
hist(female$rel.st)


## ANALYSIS ##
library(MCMCglmm)

# combine male and female data - use relative standardized fecundity values
data <- data.frame(sex = factor(c(rep("m", length(male$sex)), rep("f", length(female$sex)))), 
                   line = factor(c(male$line, female$line)), 
                   rep = factor(c(male$rep, female$rep)), 
                   rel.st = c(male$rel.st, female$rel.st))
hist(data$rel.st)

# MCMC prior and model:
prior <- list(R = list(V = diag(2)/2, nu = 0.02),
              G = list(G1 = list(V = diag(2)/2, nu = 2, alpha.mu = c(0,0), alpha.V = diag(2)/2)))
model <- MCMCglmm(rel.st ~ sex * rep - 1, 
                  random = ~us(sex):line, rcov = ~idh(sex):units, 
                  family = "gaussian", 
                  nitt = 100000, burnin = 25000, thin=50, 
                  data = data, prior = prior)
# various model checks:
autocorr.plot(model$Sol)
autocorr.plot(model$VCV)
plot(model$Sol)
plot(model$VCV)
summary(model$VCV) # this is a summary of the variance/covariance estimates from the model
# use these estimates to calculate heritability and genetic covariance:
# heritability according to Lynch and Walsh (1998) equation - i.e. the additive genetic variance in female fecundity divided by the total phenotypic variance in fecundity measured

# mean female heritability and 95% credible interval estimate, calculated from MCMC posterior distribution of estimates
femaleh2 <-  ((model$VCV[,1]) / (model$VCV[,4]+model$VCV[,1]+model$VCV[,5]+model$VCV[,6]))
mean(femaleh2) # 0.4211726
HPDinterval(femaleh2) # CI: 0.372102 - 0.4700871

# mean male heritability and 95% credible interval estimate, calculated from MCMC posterior distribution of estimates
maleh2 <-  ((model$VCV[,4]) / (model$VCV[,4]+model$VCV[,1]+model$VCV[,5]+model$VCV[,6]))
mean(maleh2) # 0.04152749
HPDinterval(maleh2) # CI: 0.01833213 - 0.06863062


# genetic covariance taken directly from the model:
mean(model$VCV[,2]) # 0.03007155
HPDinterval(model$VCV[,2]) # CI: -0.02920491 -  0.08505323


# genetic correlation between sexes, rmf, also calculated sensu Lynch and Walsh (covariance divided by square root of product of male and female variance)
rmf <- model$VCV[,2] / sqrt((model$VCV[,1]) * (model$VCV[,4]))
mean(rmf) # 0.1263571
HPDinterval(rmf) # CI: -0.1022529 - 0.3813369


## REFERENCES ##
##Lynch, M. & Walsh, B. 1998. Genetics and analysis of quantitative traits. Sinauer, Sunderland, Mass.