## ------------ Loading and formatting the data data <- read.csv("data.csv") ## ------------ Simulation parameters K <- 100 # Number of simulations vec_ind <- data[["Female.ID"]] # List of female ID inds <- unique(vec_ind) # List of unique ID I <- length(vec_ind) # Number of individuals N <- nrow(data) # Number of records mu <- 80 # Large intercept va <- 0.5 * 140 # Additive genetic variance vpe <- 1 * 140 # Permanent environment variance vmate <- 1 * 140 # Mate effect vyear <- 1 * 140 # Year effect vr <- 1.5 * 140 # Residual variance vec_mate <- data[["Male.ID"]] vec_mate[is.na(vec_mate)] <- paste0("mate", 1:length(which(is.na(vec_mate)))) vec_year <- as.character(data[["Year"]]) ## ------------ Simulating the data phens <- matrix(0, nrow = N, ncol = K) print("Simulating the phenotypes...") for (k in 1:K) { # Simulating the breeding values a <- rbv(pedigree, va) # Simulating the maternal effects pe_eff <- sapply(unique(vec_ind), function(string) {rnorm(1, 0, sqrt(vpe))}) # Simulating the social sire effects mate_eff <- sapply(unique(vec_mate), function(string) {rnorm(1, 0, sqrt(vmate))}) # Simulating the year effect y_eff <- sapply(unique(vec_year), function(string) {rnorm(1, 0, sqrt(vyear))}) # Simulating the phenotype phens[ , k] <- mu + a[vec_ind, ] + pe_eff[vec_ind] + mate_eff[vec_mate] + y_eff[vec_year] + rnorm(N, 0, sqrt(vr)) } ## ------------ Analysing the data colnames(pedigree) <- c("animal", "dam", "sire") # Defining the prior prior <- list(R = list(V=1, nu=1), G = list(G1 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000), G2 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000), G3 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000), G4 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000))) # Creating the output out <- data.frame(sim = 1:K, va.mode = NA, va.median = NA, va.low = NA, va.high = NA, vpe.mode = NA, vpe.median = NA, vpe.low = NA, vpe.high = NA, vmate.mode = NA, vmate.median = NA, vmate.low = NA, vmate.high = NA, vy.mode = NA, vy.median = NA, vy.low = NA, vy.high = NA, vr.mode = NA, vr.median = NA, vr.low = NA, vr.high = NA, h2.mode = NA, h2.median = NA, h2.low = NA, h2.high = NA) # Looping over the models for (k in 1:K) { p rint(paste0("Performing model number ", k, "...")) # Creating a dataset data <- data.frame(phen = phens[ , k], animal = vec_ind, id = vec_ind, mate = vec_mate, year = vec_year) # Running the model mod <- MCMCglmm(phen ~ 1, random = ~ animal + id + mate + year, data = data, pedigree = pedigree, prior = prior, nitt = 13000, thin = 10, burnin = 3000, verbose = FALSE) # Computing the estimates h2 <- mod[["VCV"]][,"animal"] / rowSums(mod[["VCV"]]) out[k, ] <- data.frame(sim = k, va.mode = posterior.mode(mod[["VCV"]][,"animal"]), va.median = median(mod[["VCV"]][,"animal"]), va.low = HPDinterval(mod[["VCV"]][,"animal"])[1], va.high = HPDinterval(mod[["VCV"]][,"animal"])[2], vpe.mode = posterior.mode(mod[["VCV"]][,"id"]), vpe.median = median(mod[["VCV"]][,"id"]), vpe.low = HPDinterval(mod[["VCV"]][,"id"])[1], vpe.high = HPDinterval(mod[["VCV"]][,"id"])[2], vmate.mode = posterior.mode(mod[["VCV"]][,"mate"]), vmate.median = median(mod[["VCV"]][,"mate"]), vmate.low = HPDinterval(mod[["VCV"]][,"mate"])[1], vmate.high = HPDinterval(mod[["VCV"]][,"mate"])[2], vy.mode = posterior.mode(mod[["VCV"]][,"year"]), vy.median = median(mod[["VCV"]][,"year"]), vy.low = HPDinterval(mod[["VCV"]][,"year"])[1], vy.high = HPDinterval(mod[["VCV"]][,"year"])[2], vr.mode = posterior.mode(mod[["VCV"]][,"units"]), vr.median = median(mod[["VCV"]][,"units"]), vr.low = HPDinterval(mod[["VCV"]][,"units"])[1], vr.high = HPDinterval(mod[["VCV"]][,"units"])[2], h2.mode = posterior.mode(h2), h2.median = median(h2), h2.low = HPDinterval(h2)[1], h2.high = HPDinterval(h2)[2]) } save(out, file = "power_analysis.Rdata")