In charlotte-ngs/asmss2022: Applied Statistical Methods in Animal Science - Spring Semester 2022

knitr::opts_chunk$set(echo = TRUE)

Problem 1: Milk Dataset

Use the dataset milk from package pedigreemm and fit a sire model to each of the response variables (milk, fat, prot and scs) in the data. The dataset can be loaded using the command pedigreemm::milk. The other variables like lact and herd can be used as fixed effects. The sire column is used as a random effect. For this analysis, we assume that sires are unrelated.

Your Tasks

• Analyse the milk dataset from package pedigreemm using the function lme4::lmer() for all given response variables. You can use the same model for each of the responses.
• Compute the estimated heritability for each response variable, using the fact that the heritability $h^2$ can be computed from the variance $\sigma_s^2$ of the sire effects and the phenotypic variance $\sigma_p^2$ with the formula

$$h^2 = \frac{4* \sigma_s^2}{\sigma_p^2} $$

• Compute the summary statistic using the function summary() of all the predicted sire breeding values. Solutions for the sire breeding values are obtained using the function ranef()

Solution

• Assign milk dataset to a tibble

tbl_milk <- tibble::as_tibble(pedigreemm::milk)

• Analyse the data for each of the responses. The responses are pre-defined and must appear as column names

vec_resp_milk <- c("milk", "fat", "prot", "scs")
# check
if (!all(is.element(vec_resp_milk, colnames(tbl_milk)))){
  stop(" *** ERROR: Response not in dataset: ", 
       vec_resp_milk[!is.element(vec_resp_milk, colnames(tbl_milk))])
}

For each of the responses run an analysis

vec_fix_fact_milk <- c("lact", "herd")
vec_fix_cov_milk <- c("dim")
s_rand_eff_milk <- "sire"
s_rhs <- paste0(c(paste0(vec_fix_fact_milk, collapse = " + "),
              paste0(vec_fix_cov_milk, collapse = " + "),
              paste0("(1|", s_rand_eff_milk, ")", collapse = "")),
              collapse = " + ")

Convert fixed effects to factors

# convert fixed effects to factors
for (f in vec_fix_fact_milk){
  tbl_milk[[f]] <- as.factor(tbl_milk[[f]])
}

Construct a function that runs an analysis for a given response

run_lmer_single_response <- function(ps_resp, ps_rhs, ptbl_data){
  # check
  # formula
  s_formula <- paste(ps_resp, ps_rhs, sep = " ~ ")
  frma_lmem <- as.formula(s_formula)
  return(lme4::lmer(frma_lmem, data = ptbl_data))
}

Let the analysis be done for all responses

l_lmer_results <- lapply(vec_resp_milk, 
                         run_lmer_single_response, 
                         ps_rhs = s_rhs, 
                         ptbl_data = tbl_milk)

Print all the summaries as results

for (l in l_lmer_results){
  print(summary(l))
}

• Compute estimated heritability. For the sire model, the heritability $h^2$ is defined as

$$h^2 = \frac{4*\sigma_s^2}{\sigma_s^2 + \sigma_e^2}$$ The estimates for the sire variance $\sigma_s^2$ and $\sigma_e^2$ can be extracted from the summary of the lmer result.

for (l in l_lmer_results){
  smry_l <- summary(l)
  # getting sire std. dev
  n_sire_var <- attr(smry_l$varcor[[1]], "stddev")[[1]]
  n_res_var <- smry_l$sigma
  n_h2 <- 4*n_sire_var / (n_sire_var+n_res_var)
  s_trait <- unlist(strsplit(as.character(smry_l$call)[2], split = " ~ ", fixed = TRUE))[1]
  cat("Heritability for trait: ", s_trait, ": ", n_h2, "\n")
}

• Obtain summary statistics for predicted sire breeding values. The breeding values can be obtained by the function raneff()

for (l in l_lmer_results){
  smry_l <- summary(l)
  s_trait <- unlist(strsplit(as.character(smry_l$call)[2], split = " ~ ", fixed = TRUE))[1]
  cat("Summary statistics of breeding values for trait: ", s_trait, "\n")
  print(summary(lme4::ranef(l)$sire[[1]]))
}

charlotte-ngs/asmss2022 documentation built on June 7, 2022, 1:33 p.m.