#Set working directory setwd("~/Desktop/Analysis") #Save corresponding excel sheets as csv files in working directory and add to global environment Q <- read.csv("Quantity.csv") library(arm) library(fmsb) # Determining multicolinearity for predictor variables VIF(lm(Q$Overall.MHC.Hs ~ Q$Year + Q$Age)) VIF(lm(Q$Year ~ Q$Overall.MHC.Hs + Q$Age)) VIF(lm(Q$Age ~ Q$Year + Q$Overall.MHC.Hs)) VIF(lm(Q$DBB.Heterozygosity ~ Q$DCB.Heterozygosity + Q$DAB.Heterozygosity + Q$Year + Q$Age)) VIF(lm(Q$DCB.Heterozygosity ~ Q$DBB.Heterozygosity + Q$DAB.Heterozygosity + Q$Year + Q$Age)) VIF(lm(Q$DAB.Heterozygosity ~ Q$DCB.Heterozygosity + Q$DBB.Heterozygosity + Q$Year + Q$Age)) VIF(lm(Q$Year ~ Q$DCB.Heterozygosity + Q$DBB.Heterozygosity + Q$DAB.Heterozygosity + Q$Age)) VIF(lm(Q$Age ~ Q$DCB.Heterozygosity + Q$DBB.Heterozygosity + Q$DAB.Heterozygosity + Q$Year)) # Determining if MHC heterozygosity influences copulation success overall CopSuc <- cbind(Q$Copulations, Q$Pairings - Q$Copulations) Model1 <- glm(CopSuc ~ Year + Age + Overall.MHC.Hs, family = "binomial", data = Q) summary(Model1) s.Model1 <- standardize(Model1) summary(s.Model1) M1 <- data.frame(summary(s.Model1)$coefficients) m.Overall.MHC.Hs <- mean(Q$Overall.MHC.Hs) sd.Overall.MHC.Hs <- sd(Q$Overall.MHC.Hs) intercept1 <- coef(s.Model1) [1] slope1 <- coef(s.Model1) [4] x1 <- seq(0, 1.5, 0.25) x1.m <- (x1 - m.Overall.MHC.Hs)/(2*sd.Overall.MHC.Hs) y1 <- invlogit(slope1*x1.m+intercept1)*100 intercept1.e <- summary(s.Model1)$coef[1,2] slope1.e <- summary(s.Model1)$coef[4,2] boots.i1 <- rnorm(1000, mean=intercept1, sd=intercept1.e) boots.s1 <- rnorm(1000, mean=slope1, sd=slope1.e) hist(boots.i1 + boots.s1, main= "Bootstrapped Coef") abline(v = (intercept1 + slope1), col="red", lwd=2) errors1 <- array(NA, dim = c(length(boots.i1), length(x1))) for (i in 1:length(x1)) errors1[,i] <- boots.i1 + boots.s1*x1.m[i] errors1.t <- invlogit(errors1) uppers1 <- rep(NA, length(x1)) for (i in 1:length(x1)) uppers1[i] <- quantile(errors1.t[,i], probs = 0.975) lowers1 <- rep(NA, length(x1)) for (i in 1:length(x1)) lowers1[i] <- quantile(errors1.t[,i], probs = 0.025) # Determining if MHC heterozygosity influences breeding success overall BreedSuc <- cbind(Q$Offspring, Q$Copulations - Q$Offspring) Model2 <- glm(BreedSuc ~ Year + Overall.MHC.Hs, family = "binomial", data = Q) summary(Model2) s.Model2 <- standardize(Model2) summary(s.Model2) M2 <- data.frame(summary(s.Model2)$coefficients) m.Overall.MHC.Hs <- mean(Q$Overall.MHC.Hs) sd.Overall.MHC.Hs <- sd(Q$Overall.MHC.Hs) intercept2 <- coef(s.Model2) [1] slope2 <- coef(s.Model2) [3] x2 <- seq(0, 1.5, 0.25) x2.m <- (x2 - m.Overall.MHC.Hs)/(2*sd.Overall.MHC.Hs) y2 <- (invlogit(slope2*x2.m+intercept2))*100 intercept2.e <- summary(s.Model2)$coef[1,2] slope2.e <- summary(s.Model2)$coef[3,2] boots.i2 <- rnorm(1000, mean=intercept2, sd=intercept2.e) boots.s2 <- rnorm(1000, mean=slope2, sd=slope2.e) hist(boots.i2 + boots.s2, main= "Bootstrapped Coef") abline(v = (intercept2 + slope2), col="red", lwd=2) errors2 <- array(NA, dim = c(length(boots.i2), length(x2))) for (i in 1:length(x2)) errors2[,i] <- boots.i2 + boots.s2*x2.m[i] errors2.t <- invlogit(errors2) uppers2 <- rep(NA, length(x2)) for (i in 1:length(x2)) uppers2[i] <- quantile(errors2.t[,i], probs = 0.975) lowers2 <- rep(NA, length(x2)) for (i in 1:length(x2)) lowers2[i] <- quantile(errors2.t[,i], probs = 0.025) # Determining if MHC heterozygosity influences offspring success overall OffSuc <- cbind(Q$Surviving.Offspring, Q$Offspring - Q$Surviving.Offspring) Model3 <- glm(OffSuc ~ Year + Overall.MHC.Hs, family = "binomial", data = Q) summary(Model3) s.Model3 <- standardize(Model3) summary(s.Model3) M3 <- data.frame(summary(s.Model3)$coefficients) # Determining if heterozygosity at individual MHC Loci influences copulation success DBBCopSuc <- cbind(Q$Copulations, Q$Pairings - Q$Copulations) Model4 <- glm(DBBCopSuc ~ Year + Age + DBB.Heterozygosity + DCB.Heterozygosity + DAB.Heterozygosity, family = "binomial", data = Q) summary(Model4) s.Model4 <- standardize(Model4) summary(s.Model4) M4 <- data.frame(summary(s.Model4)$coefficients) m.DAB.Heterozygosity <- mean(Q$DAB.Heterozygosity) sd.DAB.Heterozygosity <- sd(Q$DAB.Heterozygosity) intercept4 <- coef(s.Model4) [1] slope4 <- coef(s.Model4) [6] x4 <- seq(0, 1, 0.1) x4.m <- (x4 - m.DAB.Heterozygosity)/(2*sd.DAB.Heterozygosity) y4 <- invlogit(slope4*x4.m+intercept4)*100 intercept4.e <- summary(s.Model4)$coef[1,2] slope4.e <- summary(s.Model4)$coef[6,2] boots.i4 <- rnorm(1000, mean=intercept4, sd=intercept4.e) boots.s4 <- rnorm(1000, mean=slope4, sd=slope4.e) hist(boots.i4 + boots.s4, main= "Bootstrapped Coef") abline(v = (intercept4 + slope4), col="red", lwd=2) errors4 <- array(NA, dim = c(length(boots.i4), length(x4))) for (i in 1:length(x4)) errors4[,i] <- boots.i4 + boots.s4*x4.m[i] errors4.t <- invlogit(errors4) uppers4 <- rep(NA, length(x4)) for (i in 1:length(x4)) uppers4[i] <- quantile(errors4.t[,i], probs = 0.975) lowers4 <- rep(NA, length(x4)) for (i in 1:length(x4)) lowers4[i] <- quantile(errors4.t[,i], probs = 0.025) # Determining if heterozygosity at indiviudal MHC Loci influences breeding success DBBBreedSuc <- cbind(Q$Offspring, Q$Copulations - Q$Offspring) Model5 <- glm(DBBBreedSuc ~ Year + DBB.Heterozygosity + DCB.Heterozygosity + DAB.Heterozygosity, family = "binomial", data = Q) summary(Model5) s.Model5 <- standardize(Model5) summary(s.Model5) M5 <- data.frame(summary(s.Model5)$coefficients) # Determining if heterozygosity at individual MHC Loci influences offspring success DBBOffSuc <- cbind(Q$Surviving.Offspring, Q$Offspring - Q$Surviving.Offspring) Model6 <- glm(DBBOffSuc ~ Year + DBB.Heterozygosity, family = "binomial", data = Q) summary(Model6) s.Model6 <- standardize(Model6) summary(s.Model6) M6 <- data.frame(summary(s.Model6)$coefficients) #Plotting all the graphs for MHC heterozygosity par(mfrow = c(3,4), family = "sans", cex.axis = 1.2, cex.lab = 1.5, mar = c(2, 2, 2, 2), oma = c(4,4,4,0), xpd = NA) plot(x = Q$Overall.MHC.Hs, y = Q$Copulations/Q$Pairings*100, cex = (Q$Pairings)/20, xlab ="", xaxt = "n", ylab = "Copulation Success", ylim = c(0,100), xlim = c(0,1.5)) axis (side = 1, at = seq(0,1.5,0.5), labels = F) points(x = x1, y = y1, type = "l", lwd = 2) points(x = x1, y = uppers1*100, type = "l", lty = "dotted") points(x = x1, y = lowers1*100, type = "l", lty = "dotted") plot(x = Q$DBB.Heterozygosity, y = Q$Copulation.Success, cex = (Q$Pairings)/20, xlab = "", xaxt = "n", ylab = "", yaxt = "n", xlim = c(-0.1,1.1), ylim = c(0,100)) axis(side = 1, at = c(0,1), labels = F) axis(side = 2, at = seq(0,100,20), labels = F) plot(x = Q$DCB.Heterozygosity, y = Q$Copulation.Success, cex = (Q$Pairings)/20, xlab = "", xaxt = "n", ylab = "", yaxt = "n", xlim = c(-0.1,1.1), ylim = c(0,100)) axis(side = 1, at = c(0,1), labels = F) axis(side = 2, at = seq(0,100,20), labels = F) plot(x = Q$DAB.Heterozygosity, y = Q$Copulation.Success, cex = (Q$Pairings)/20, xlab = "", xaxt = "n", ylab = "", yaxt = "n", xlim = c(-0.1,1.1), ylim = c(0,100)) axis(side = 1, at = c(0,1), labels = F) axis(side = 2, at = seq(0,100,20), labels = F) points(x = x4, y = y4, type = "l", lwd = 2) points(x = x4, y = uppers4*100, type = "l", lty = "dotted") points(x = x4, y = lowers4*100, type = "l", lty = "dotted") plot(x = Q$Overall.MHC.Hs, y = Q$Offspring/Q$Copulations*100, cex = (Q$Copulations)/10, xlab ="", xaxt = "n", ylab = "Breeding Success", ylim = c(0,100), xlim = c(0,1.5)) axis (side = 1, at = seq(0,1.5,0.5), labels = F) points(x = x2, y = y2, type = "l", lwd = 2) points(x = x2, y = uppers2*100, type = "l", lty = "dotted") points(x = x2, y = lowers2*100, type = "l", lty = "dotted") plot(x = Q$DBB.Heterozygosity, y = Q$Breeding.Success, cex = (Q$Copulations)/10, xlab = "", xaxt = "n", ylab = "", yaxt = "n", xlim = c(-0.1,1.1), ylim = c(0,100)) axis(side = 1, at = c(0,1), labels = F) axis(side = 2, at = seq(0,100,20), labels = F) plot(x = Q$DCB.Heterozygosity, y = Q$Breeding.Success, cex = (Q$Copulations)/10, xlab = "", xaxt = "n", ylab = "", yaxt = "n", xlim = c(-0.1,1.1), ylim = c(0,100)) axis(side = 1, at = c(0,1), labels = F) axis(side = 2, at = seq(0,100,20), labels = F) plot(x = Q$DAB.Heterozygosity, y = Q$Breeding.Success, cex = (Q$Copulations)/10, xlab = "", xaxt = "n", ylab = "", yaxt = "n", xlim = c(-0.1,1.1), ylim = c(0,100)) axis(side = 1, at = c(0,1), labels = F) axis(side = 2, at = seq(0,100,20), labels = F) plot(x = Q$Overall.MHC.Hs, y = Q$Surviving.Offspring/Q$Offspring*100, cex = (Q$Offspring)/5, xlab ="Overall", ylab = "Offspring Success", ylim = c(0,100), xlim = c(0,1.5), xaxt = "n") axis (side = 1, at = seq(0,1.5,0.5)) plot(x = Q$DBB.Heterozygosity, y = Q$Offspring.Success, cex = (Q$Offspring)/5, xlab = "DBB", xaxt = "n", ylab = "", yaxt = "n", xlim = c(-0.1,1.1), ylim = c(0,100)) axis(side = 1, at = c(0,1)) axis(side = 2, at = seq(0,100,20), labels = F) plot(x = Q$DCB.Heterozygosity, y = Q$Offspring.Success, cex = (Q$Offspring)/5, xlab = "DCB",xaxt = "n", ylab = "", yaxt = "n", xlim = c(-0.1,1.1), ylim = c(0,100)) axis(side = 1, at = c(0,1)) axis(side = 2, at = seq(0,100,20), labels = F) plot(x = Q$DAB.Heterozygosity, y = Q$Offspring.Success, cex = (Q$Offspring)/5, xlab = "DAB", xaxt = "n", ylab = "", yaxt = "n", xlim = c(-0.1,1.1), ylim = c(0,100)) axis(side = 1, at = c(0,1)) axis(side = 2, at = seq(0,100,20), labels = F) # Determining if genome-wide heterozygosity influences copulation success overall CopSuc <- cbind(Q$Copulations, Q$Pairings - Q$Copulations) Model1 <- glm(CopSuc ~ Year + Age + Non.MHC.Hs, family = "binomial", data = Q) summary(Model1) s.Model1 <- standardize(Model1) summary(s.Model1) M7 <- data.frame(summary(s.Model1)$coefficients) # Determining if genome-wide heterozygosity influences breeding success overall BreedSuc <- cbind(Q$Offspring, Q$Copulations - Q$Offspring) Model2 <- glm(BreedSuc ~ Year + Non.MHC.Hs, family = "binomial", data = Q) summary(Model2) s.Model2 <- standardize(Model2) summary(s.Model2) M8 <- data.frame(summary(s.Model2)$coefficients) # Determining if genome-wide heterozygosity influences offspring success overall OffSuc <- cbind(Q$Surviving.Offspring, Q$Offspring - Q$Surviving.Offspring) Model3 <- glm(OffSuc ~ Year + Non.MHC.Hs, family = "binomial", data = Q) summary(Model3) s.Model3 <- standardize(Model3) summary(s.Model3) M9 <- data.frame(summary(s.Model3)$coefficients) #Plotting all the graphs for genome-wide heterozygosity par(mfrow = c(3,1), family = "sans") plot(x = Q$Non.MHC.Hs, y = Q$Copulations/Q$Pairings*100, cex = (Q$Pairings)/20, xlab ="", xaxt = "n", ylab = "Copulation Success") axis (side = 1, labels = F) plot(x = Q$Overall.MHC.Hs, y = Q$Offspring/Q$Copulations*100, cex = (Q$Copulations)/10, xlab ="", xaxt = "n", ylab = "Breeding Success") axis (side = 1, labels = F) plot(x = Q$Non.MHC.Hs, y = Q$Surviving.Offspring/Q$Offspring*100, cex = (Q$Offspring)/5, xlab ="Genome-wide Heterozygosity", ylab = "Offspring Success") # Extracting the results Results <- rbind(M1, M2, M3, M4, M5, M6, M7, M8, M9) write.csv(data.frame(Results), file = "QuantityResults.csv")