#Set working directory
setwd("~/Desktop/Analysis")

#Save corresponding excel sheets as csv files in working directory and add to global environment
A <- read.csv("Alleles.csv")
library(arm)

# Base only Models
CS <- cbind(A$Copulations, A$Pairings - A$Copulations) #Copulation success response variable
BM1 <- glm(CS ~ A$Year + A$Age +  1, family = "binomial")
summary(BM1)

BS <- cbind(A$Offspring, A$Copulations - A$Offspring) #Breeding success response variable
BM2 <- glm(BS ~ A$Year + 1, family = "binomial")
summary(BM2)

OS <- cbind(A$Surviving.Offspring, A$Offspring - A$Surviving.Offspring) #Offspring success response variable
BM3 <- glm(OS ~ A$Year + 1, family = "binomial")
summary(BM3)

# Determining if females prefer to copulate with males that have specific DBB alleles
CS <- cbind(A$Copulations, A$Pairings - A$Copulations)
M1 <- glm(CS ~ A$Year + A$Age +  A$DBB277, family = "binomial")
summary(M1)
M2 <- glm(CS ~ A$Year + A$Age +  A$DBB287, family = "binomial")
summary(M2)
M3 <- glm(CS ~ A$Year + A$Age +  A$DBB289, family = "binomial")
summary(M3)
M4 <- glm(CS ~ A$Year + A$Age +  A$DBB297, family = "binomial")
summary(M4)

# Determining if females prefer to copulate with males that have specific DCB alleles
CS <- cbind(A$Copulations, A$Pairings - A$Copulations)
M5 <- glm(CS ~ A$Year + A$Age +  A$DCB220, family = "binomial")
summary(M5)
M6 <- glm(CS ~ A$Year + A$Age +  A$DCB226, family = "binomial")
summary(M6)
M7 <- glm(CS ~ A$Year + A$Age +  A$DCB228, family = "binomial")
summary(M7)
M8 <- glm(CS ~ A$Year + A$Age +  A$DCB250, family = "binomial")
summary(M8)
M9 <- glm(CS ~ A$Year + A$Age +  A$DCB252, family = "binomial")
summary(M9)
M10 <- glm(CS ~ A$Year + A$Age +  A$DCB254, family = "binomial")
summary(M10)
M11 <- glm(CS ~ A$Year + A$Age +  A$DCB256, family = "binomial")
summary(M11)
M12 <- glm(CS ~ A$Year + A$Age +  A$DCB260, family = "binomial")
summary(M12)
M13 <- glm(CS ~ A$Year + A$Age +  A$DCB266, family = "binomial")
summary(M13)

# Determining if females prefer to copulate with males that have specific DAB alleles
CS <- cbind(A$Copulations, A$Pairings - A$Copulations)
M14 <- glm(CS ~ A$Year + A$Age +  A$DAB281, family = "binomial")
summary(M14)
M15 <- glm(CS ~ A$Year + A$Age +  A$DAB285, family = "binomial")
summary(M15)
M16 <- glm(CS ~ A$Year + A$Age +  A$DAB287, family = "binomial")
summary(M16)
M17 <- glm(CS ~ A$Year + A$Age +  A$DAB289, family = "binomial")
summary(M17)
M18 <- glm(CS ~ A$Year + A$Age +  A$DAB291, family = "binomial")
summary(M18)
M19 <- glm(CS ~ A$Year + A$Age +  A$DAB293, family = "binomial")
summary(M19)
M20 <- glm(CS ~ A$Year + A$Age +  A$DAB297, family = "binomial")
summary(M20)

# Saving copulation results
results <- array(NA, dim = c(20, 4))
results <- as.data.frame(results)
colnames(results) <- c("Allele", "Slope", "SE", "AIC")
results$Allele <- c("DBB277", "DBB287", "DBB289", "DBB297", 
                    "DCB220", "DCB226", "DCB228", "DCB250", "DCB252", "DCB254", "DCB256", "DCB260", "DCB266",
                    "DCB281", "DCB285", "DCB287", "DCB289", "DCB291", "DCB293", "DCB297")

model.list <- list(M1, M2, M3, M4, M5, M6, M7, M8, M9, M10, M11, M12, M13, M14, M15, M16, M17, M18, M19, M20)
AIC.data <- lapply(X = model.list, FUN = AIC)
AIC.data <- as.numeric(AIC.data)
results$AIC <- AIC.data

results$Slope <- c(summary(M1)$coef [4,1],summary(M2)$coef [4,1],summary(M3)$coef [4,1],summary(M4)$coef [4,1],summary(M5)$coef [4,1],
                   summary(M6)$coef [4,1],summary(M7)$coef [4,1],summary(M8)$coef [4,1],summary(M9)$coef [4,1],summary(M10)$coef [4,1],
                   summary(M11)$coef [4,1],summary(M12)$coef [4,1],summary(M13)$coef [4,1],summary(M14)$coef [4,1],summary(M15)$coef [4,1],
                   summary(M16)$coef [4,1],summary(M17)$coef [4,1],summary(M18)$coef [4,1],summary(M19)$coef [4,1],summary(M20)$coef [4,1])

results$SE <- c(summary(M1)$coef [4,2],summary(M2)$coef [4,2],summary(M3)$coef [4,2],summary(M4)$coef [4,2],summary(M5)$coef [4,2],
                summary(M6)$coef [4,2],summary(M7)$coef [4,2],summary(M8)$coef [4,2],summary(M9)$coef [4,2],summary(M10)$coef [4,2],
                summary(M11)$coef [4,2],summary(M12)$coef [4,2],summary(M13)$coef [4,2],summary(M14)$coef [4,2],summary(M15)$coef [4,2],
                summary(M16)$coef [4,2],summary(M17)$coef [4,2],summary(M18)$coef [4,2],summary(M19)$coef [4,2],summary(M20)$coef [4,2])

write.csv(results, file = "AllelesCopResults.csv")

# Determining if females prefer to breed with males that have specific DBB alleles
BS <- cbind(A$Offspring, A$Copulations - A$Offspring)
M2.1 <- glm(BS ~ A$Year + A$DBB277, family = "binomial")
summary(M2.1)
M2.2 <- glm(BS ~ A$Year + A$DBB287, family = "binomial")
summary(M2.2)
M2.3 <- glm(BS ~ A$Year + A$DBB289, family = "binomial")
summary(M2.3)
M2.4 <- glm(BS ~ A$Year + A$DBB297, family = "binomial")
summary(M2.4)

# Determining if females prefer to breed with males that have specific DCB alleles
BS <- cbind(A$Offspring, A$Copulations - A$Offspring)
M2.5 <- glm(BS ~ A$Year +  A$DCB220, family = "binomial")
summary(M2.5)
M2.6 <- glm(BS ~ A$Year +  A$DCB226, family = "binomial")
summary(M2.6)
M2.7 <- glm(BS ~ A$Year +  A$DCB228, family = "binomial")
summary(M2.7)
M2.8 <- glm(BS ~ A$Year +  A$DCB250, family = "binomial")
summary(M2.8)
M2.9 <- glm(BS ~ A$Year +  A$DCB252, family = "binomial")
summary(M2.9)
M2.10 <- glm(BS ~ A$Year +  A$DCB254, family = "binomial")
summary(M2.10)
M2.11 <- glm(BS ~ A$Year +  A$DCB256, family = "binomial")
summary(M2.11)
M2.12 <- glm(BS ~ A$Year +  A$DCB260, family = "binomial")
summary(M2.12)
M2.13 <- glm(BS ~ A$Year +  A$DCB266, family = "binomial")
summary(M2.13)

# Determining if females prefer to breed with males that have specific DAB alleles
BS <- cbind(A$Offspring, A$Copulations - A$Offspring)
M2.14 <- glm(BS ~ A$Year +  A$DAB281, family = "binomial")
summary(M2.14)
M2.15 <- glm(BS ~ A$Year +  A$DAB285, family = "binomial")
summary(M2.15)
M2.16 <- glm(BS ~ A$Year +  A$DAB287, family = "binomial")
summary(M2.16)
M2.17 <- glm(BS ~ A$Year +  A$DAB289, family = "binomial")
summary(M2.17)
M2.18 <- glm(BS ~ A$Year +  A$DAB291, family = "binomial")
summary(M2.18)
M2.19 <- glm(BS ~ A$Year +  A$DAB293, family = "binomial")
summary(M2.19)
M2.20 <- glm(BS ~ A$Year +  A$DAB297, family = "binomial")
summary(M2.20)

# Saving breeding results
results2 <- array(NA, dim = c(20, 4))
results2 <- as.data.frame(results2)
colnames(results2) <- c("Allele", "Slope", "SE", "AIC")
results2$Allele <- c("DBB277", "DBB287", "DBB289", "DBB297", 
                     "DCB220", "DCB226", "DCB228", "DCB250", "DCB252", "DCB254", "DCB256", "DCB260", "DCB266",
                     "DCB281", "DCB285", "DCB287", "DCB289", "DCB291", "DCB293", "DCB297")

model.list2 <- list(M2.1, M2.2, M2.3, M2.4, M2.5, M2.6, M2.7, M2.8, M2.9, M2.10, M2.11, M2.12, M2.13, M2.14, M2.15, M2.16, M2.17, M2.18, M2.19, M2.20)
AIC.data2 <- lapply(X = model.list2, FUN = AIC)
AIC.data2 <- as.numeric(AIC.data2)
results2$AIC <- AIC.data2

results2$Slope <- c(summary(M2.1)$coef [3,1],summary(M2.2)$coef [3,1],summary(M2.3)$coef [3,1],summary(M2.4)$coef [3,1],summary(M2.5)$coef [3,1],
                    summary(M2.6)$coef [3,1],summary(M2.7)$coef [3,1],summary(M2.8)$coef [3,1],summary(M2.9)$coef [3,1],summary(M2.10)$coef [3,1],
                    summary(M2.11)$coef [3,1],summary(M2.12)$coef [3,1],summary(M2.13)$coef [3,1],summary(M2.14)$coef [3,1],summary(M2.15)$coef [3,1],
                    summary(M2.16)$coef [3,1],summary(M2.17)$coef [3,1],summary(M2.18)$coef [3,1],summary(M2.19)$coef [3,1],summary(M2.20)$coef [3,1])

results2$SE <- c(summary(M2.1)$coef [3,2],summary(M2.2)$coef [3,2],summary(M2.3)$coef [3,2],summary(M2.4)$coef [3,2],summary(M2.5)$coef [3,2],
                 summary(M2.6)$coef [3,2],summary(M2.7)$coef [3,2],summary(M2.8)$coef [3,2],summary(M2.9)$coef [3,2],summary(M2.10)$coef [3,2],
                 summary(M2.11)$coef [3,2],summary(M2.12)$coef [3,2],summary(M2.13)$coef [3,2],summary(M2.14)$coef [3,2],summary(M2.15)$coef [3,2],
                 summary(M2.16)$coef [3,2],summary(M2.17)$coef [3,2],summary(M2.18)$coef [3,2],summary(M2.19)$coef [3,2],summary(M2.20)$coef [3,2])

write.csv(results2, file = "AllelesBreedResults.csv")

# Determining if offspring survive better from males that have specific DBB alleles
OS <- cbind(A$Surviving.Offspring, A$Offspring - A$Surviving.Offspring)
M3.1 <- glm(OS ~ A$Year + A$DBB277, family = "binomial")
summary(M3.1)
M3.2 <- glm(OS ~ A$Year + A$DBB287, family = "binomial")
summary(M3.2)
M3.3 <- glm(OS ~ A$Year + A$DBB289, family = "binomial")
summary(M3.3)
M3.4 <- glm(OS ~ A$Year + A$DBB297, family = "binomial")
summary(M3.4)


# Determining if offspring survive better from males that have specific DCB alleles
OS <- cbind(A$Surviving.Offspring, A$Offspring - A$Surviving.Offspring)
M3.5 <- glm(OS ~ A$Year +  A$DCB220, family = "binomial")
summary(M3.5)
M3.6 <- glm(OS ~ A$Year +  A$DCB226, family = "binomial")
summary(M3.6)
M3.7 <- glm(OS ~ A$Year +  A$DCB228, family = "binomial")
summary(M3.7)
M3.8 <- glm(OS ~ A$Year +  A$DCB250, family = "binomial")
summary(M3.8)
M3.9 <- glm(OS ~ A$Year +  A$DCB252, family = "binomial")
summary(M3.9)
M3.10 <- glm(OS ~ A$Year +  A$DCB254, family = "binomial")
summary(M3.10)
M3.11 <- glm(OS ~ A$Year +  A$DCB256, family = "binomial")
summary(M3.11)
M3.12 <- glm(OS ~ A$Year +  A$DCB260, family = "binomial")
summary(M3.12)
M3.13 <- glm(OS ~ A$Year +  A$DCB266, family = "binomial")
summary(M3.13)

# Determining if offspring survive better from males that have specific DAB alleles
OS <- cbind(A$Surviving.Offspring, A$Offspring - A$Surviving.Offspring)
M3.14 <- glm(OS ~ A$Year +  A$DAB281, family = "binomial")
summary(M3.14)
M3.15 <- glm(OS ~ A$Year +  A$DAB285, family = "binomial")
summary(M3.15)
M3.16 <- glm(OS ~ A$Year +  A$DAB287, family = "binomial")
summary(M3.16)
M3.17 <- glm(OS ~ A$Year +  A$DAB289, family = "binomial")
summary(M3.17)
M3.18 <- glm(OS ~ A$Year +  A$DAB291, family = "binomial")
summary(M3.18)
M3.19 <- glm(OS ~ A$Year +  A$DAB293, family = "binomial")
summary(M3.19)
M3.20 <- glm(OS ~ A$Year +  A$DAB297, family = "binomial")
summary(M3.20)

# Saving offspring results
results3 <- array(NA, dim = c(20, 4))
results3 <- as.data.frame(results3)
colnames(results3) <- c("Allele", "Slope", "SE", "AIC")
results3$Allele <- c("DBB277", "DBB287", "DBB289", "DBB297", 
                     "DCB220", "DCB226", "DCB228", "DCB250", "DCB252", "DCB254", "DCB256", "DCB260", "DCB266",
                     "DCB281", "DCB285", "DCB287", "DCB289", "DCB291", "DCB293", "DCB297")

model.list3 <- list(M3.1, M3.2, M3.3, M3.4, M3.5, M3.6, M3.7, M3.8, M3.9, M3.10, M3.11, M3.12, M3.13, M3.14, M3.15, M3.16, M3.17, M3.18, M3.19, M3.20)
AIC.data3 <- lapply(X = model.list3, FUN = AIC)
AIC.data3 <- as.numeric(AIC.data3)
results3$AIC <- AIC.data3

results3$Slope <- c(summary(M3.1)$coef [3,1],summary(M3.2)$coef [3,1],summary(M3.3)$coef [3,1],summary(M3.4)$coef [3,1],summary(M3.5)$coef [3,1],
                    summary(M3.6)$coef [3,1],summary(M3.7)$coef [3,1],summary(M3.8)$coef [3,1],summary(M3.9)$coef [3,1],summary(M3.10)$coef [3,1],
                    summary(M3.11)$coef [3,1],summary(M3.12)$coef [3,1],summary(M3.13)$coef [3,1],summary(M3.14)$coef [3,1],summary(M3.15)$coef [3,1],
                    summary(M3.16)$coef [3,1],summary(M3.17)$coef [3,1],summary(M3.18)$coef [3,1],summary(M3.19)$coef [3,1],summary(M3.20)$coef [3,1])

results3$SE <- c(summary(M3.1)$coef [3,2],summary(M3.2)$coef [3,2],summary(M3.3)$coef [3,2],summary(M3.4)$coef [3,2],summary(M3.5)$coef [3,2],
                 summary(M3.6)$coef [3,2],summary(M3.7)$coef [3,2],summary(M3.8)$coef [3,2],summary(M3.9)$coef [3,2],summary(M3.10)$coef [3,2],
                 summary(M3.11)$coef [3,2],summary(M3.12)$coef [3,2],summary(M3.13)$coef [3,2],summary(M3.14)$coef [3,2],summary(M3.15)$coef [3,2],
                 summary(M3.16)$coef [3,2],summary(M3.17)$coef [3,2],summary(M3.18)$coef [3,2],summary(M3.19)$coef [3,2],summary(M3.20)$coef [3,2])

write.csv(results3, file = "AllelesOffResults.csv")