rm(list=ls())

myPath <- '/Path/to/my/files'

# point colors

control_col <- '#ffffbf'
b_downCol <- '#2c7bb6'
b_upCol <- '#d7191c'

#################
# Figure 2 code #
#################

repro_over <- read.table(paste0(myPath, 'repro_mat_over'))
repro_under <- read.table(paste0(myPath, 'repro_mat_under'))

quartz(width=8, height=4)

par(mfrow=c(1, 2))

par(mar=c(5, 4, 2, 1))

# Panel A: uncorrected values

plot(c(-0.5, 5), c(-0.5, 5), type='n', xlab='real reproductive difference', ylab='measured', axes=F, main='uncorrected')

axis(1, 0:5, 0:5)
axis(2, 0:5, 0:5)


realDiff_over <- repro_over$full_uninf - repro_over$full_inf

biasedDiff_over <- repro_over$cap_uninf - repro_over$cap_inf

unbiasedDiff_over <- repro_over$cont_uninf - repro_over$cont_inf

realDiff_under <- repro_under$full_uninf - repro_under$full_inf

biasedDiff_under <- repro_under$cap_uninf - repro_under$cap_inf

points(realDiff_over[1:50], biasedDiff_over[1:50], pch=21, bg= b_upCol)
points(realDiff_over[1:50], unbiasedDiff_over[1:50], pch= 21, bg= control_col)
points(realDiff_under[1:50], biasedDiff_under[1:50], pch= 21, bg= b_downCol)

legend('topleft', fill=c(control_col, b_downCol, b_upCol), legend=c('control', 'low rate', 'high rate'), bty='n')

# Panel B: corrected values

# divide the difference in reproductive success by the mean number of offspring captured from uninfected parents

realDiff_corr <- realDiff_over/repro_over$full_uninf

biasedDiff_corr <- biasedDiff_over/repro_over$cap_uninf

unbiasedDiff_corr <- unbiasedDiff_over/repro_over$cont_uninf

realDiff_down_corr <- realDiff_under/repro_under$full_uninf

biasedDiff_down_corr <- biasedDiff_under/repro_under$cap_uninf

plot(c(-0.2, 0.8), c(-0.2, 0.8), type='n', xlab='real reproductive difference', ylab='measured', axes=F, main='corrected')

axis(1, seq(-0.2, 0.8, by=0.2), seq(-0.2, 0.8, by=0.2))
axis(2, seq(-0.2, 0.8, by=0.2), seq(-0.2, 0.8, by=0.2))


points(realDiff_corr[1:50], biasedDiff_corr[1:50], pch=21, bg=b_upCol)
points(realDiff_corr[1:50], unbiasedDiff_corr[1:50], pch=21, bg= control_col)
points(realDiff_down_corr[1:50], biasedDiff_down_corr[1:50], pch=21, bg= b_downCol)

#################
# Figure 3 code #
#################

het_geno_down <- read.table(paste0(myPath, 'RR_under_allele'))

het_geno_up <-  read.table(paste0(myPath, 'RR_over_allele'))

allele_geno_down <- read.table(paste0(myPath, 'RR_under_het'))

allele_geno_up <- read.table(paste0(myPath, 'RR_over_het'))


quartz(width=8, height=4)

par(mfrow=c(1, 2))

plot(c(1, 7.5), c(1, 7.5), type='n', main='heterozygote advantage', axes=F, xlab='real relative risk', ylab='sampled relative risk')

axis(1, seq(1, 7, by=1), seq(1, 7, by=1))

axis(2, seq(1, 7, by=1), seq(1, 7, by=1))

points(het_geno_down$full_RR[1:50], het_geno_down$cap_RR[1:50], bg= b_downCol, pch=21, cex=1)
points(het_geno_down$full_RR[1:50], het_geno_down$cont_RR[1:50], bg= control_col, pch=21, cex=1)
points(het_geno_up$full_RR[1:50], het_geno_up$cap_RR[1:50], bg= b_upCol, pch=21, cex=1)

legend('topleft', fill=c(control_col, b_downCol, b_upCol), legend=c('control', 'low rate', 'high rate'), bty='n')	

plot(c(1, 5), c(1, 5), type='n', main='resistance allele', axes=F, xlab='real relative risk', ylab='sampled relative risk')

axis(1, seq(1, 5, by=1), seq(1, 5, by=1))

axis(2, seq(1, 5, by=1), seq(1, 5, by=1))

points(allele_geno_down$full_RR[1:50], allele_geno_down$cap_RR[1:50], bg= b_downCol, pch=21, cex=1)
points(allele_geno_down$full_RR[1:50], allele_geno_down$cont_RR[1:50], bg= control_col, pch=21, cex=1)
points(allele_geno_up$full_RR[1:50], allele_geno_up$cap_RR[1:50], bg= b_upCol, pch=21, cex=1)


#################
# Figure 4 code #
#################

allele_change_down <- read.table(paste0(myPath, 'allele_change_under_allele'))

allele_change_up <-  read.table(paste0(myPath, 'allele_change_over_allele'))

het_change_down <- read.table(paste0(myPath, 'allele_change_under_het'))

het_change_up <- read.table(paste0(myPath, 'allele_change_over_het'))

quartz(width=8, height=4)

par(mfrow=c(1, 2))

plot(c(-0.02, 0.02), c(-0.06, 0.06), main='heterozygote advantage', type='n', xlab='true change in frequency', ylab='measured change in frequency', axes=F)

axis(1, seq(-0.02, 0.02, by = 0.01), seq(-0.02, 0.02, by = 0.01))

axis(2, seq(-0.06, 0.06, by = 0.03), seq(-0.06, 0.06, by = 0.03), las=2)

points(het_change_down$full_change[1:50], het_change_down$cap_change[1:50],  bg= b_downCol, pch=21)

points(het_change_up$full_change[1:50], het_change_up$cont_change[1:50], bg= control_col, pch=21)

points(het_change_up$full_change[1:50], het_change_up$cap_change[1:50],  bg= b_upCol, pch=21)

legend('topleft', fill=c(control_col, b_downCol, b_upCol), legend=c('control', 'low rate', 'high rate'), bty='n')


plot(c(-0.04, 0.04), c(-0.07, 0.05), main='resistance allele', type='n', xlab='true change in frequency', ylab='measured change in frequency', axes=F)

axis(1, seq(-0.06, 0.04, by = 0.02), seq(-0.06, 0.04, by = 0.02))

axis(2, seq(-0.04, 0.04, by = 0.02), seq(-0.04, 0.04, by = 0.02), las=2)

points(allele_change_down$full_change[1:50], allele_change_down$cap_change[1:50], bg= b_downCol, cex=1, pch=21)

points(allele_change_up$full_change[1:50], allele_change_up$cont_change[1:50], bg= control_col, cex=1, pch=21)

points(allele_change_up$full_change[1:50], allele_change_up$cap_change[1:50],  bg= b_upCol, pch=21)



#################
# Figure 5 code #
#################

repro_over <- read.table(paste0(myPath, 'repro_mat_over'))
repro_under <- read.table(paste0(myPath, 'repro_mat_under'))


repro_over_resid <- resid(lm((repro_over$cap_uninf - repro_over$cap_inf) ~ (repro_over$full_uninf - repro_over$full_inf)))
repro_cont_resid <- resid(lm((repro_over$cont_uninf - repro_over$cont_inf) ~ (repro_over$full_uninf - repro_over$full_inf)))
repro_under_resid <- resid(lm((repro_under $cap_uninf - repro_under $cap_inf) ~ (repro_under $full_uninf - repro_under $full_inf)))


p_diff_over <- repro_over$cap_p_uninf - repro_over$cap_p_inf
p_diff_under <- repro_under$cap_p_uninf - repro_under$cap_p_inf
p_diff_cont <- repro_under$cont_p_uninf - repro_under$cont_p_inf

############


allele_geno_down <- read.table(paste0(myPath, 'RR_under_het'))

allele_geno_up <- read.table(paste0(myPath, 'RR_over_het'))


RR_over_resid <- resid(lm(allele_geno_up $cap_RR  ~ allele_geno_up $full_RR ))
RR_cont_resid <- resid(lm(allele_geno_up $cont_RR  ~ allele_geno_up $full_RR ))
RR_under_resid <- resid(lm(allele_geno_down $cap_RR  ~ allele_geno_down $full_RR ))

p_diff_over_RR <- allele_geno_up $cap_p_uninf - allele_geno_up $cap_p_inf
p_diff_under_RR <- allele_geno_down $cap_p_uninf - allele_geno_down $cap_p_inf
p_diff_cont_RR <- allele_geno_down $cont_p_uninf - allele_geno_down $cont_p_inf


############

quartz(width=8, height=4)

par(mfrow=c(1, 2))

plot(c(-0.3, 0.3), c(-0.04, 0.06), main='reproductive success', type='n', xlab='residuals', ylab='capture probability difference', axes=F)

axis(1, seq(-0.3, 0.3, by = 0.2), seq(-0.3, 0.3, by = 0.2))

axis(2, seq(-0.04, 0.06, by = 0.02), seq(-0.04, 0.06, by = 0.02), las=2)

points(repro_under_resid[1:50], p_diff_under[1:50], bg= b_downCol, cex=1, pch=21)

points(repro_cont_resid[1:50], p_diff_cont[1:50], bg= control_col, cex=1, pch=21)

points(repro_over_resid[1:50], p_diff_over[1:50],  bg= b_upCol, pch=21)


plot(c(-0.4, 0.4), c(-0.04, 0.06), main='relative risk', type='n', xlab='residuals', ylab='capture probability difference', axes=F)

axis(1, seq(-0.4, 0.4, by = 0.2), seq(-0.4, 0.4, by = 0.2))

axis(2, seq(-0.04, 0.06, by = 0.02), seq(-0.04, 0.06, by = 0.02), las=2)

points(RR_under_resid[1:50], p_diff_under_RR[1:50], bg= b_downCol, cex=1, pch=21)

points(RR_cont_resid[1:50], p_diff_cont_RR[1:50], bg= control_col, cex=1, pch=21)

points(RR_over_resid[1:50], p_diff_over_RR[1:50],  bg= b_upCol, pch=21)