In pvrqualitasag/MeatValueIndex:

Disclaimer

This notebook can be used for testing functions to compute asymptotic selection responses.

Debugging of function asym_sel_res.R

Number of traits (ccc cca cfc cfa cwc cwa)

#s_a <- "../inst/extdata/economic_values/AN_economic_values_Silvan"
s_a <- "../inst/extdata/economic_values/LM_economic_values_Silvan"

df_a <- read.table(file = s_a, sep=";", header=TRUE)
mat_a <- as.matrix(df_a)
pvec_economic_values <- as.vector(mat_a)


n_r_trait <- length(pvec_economic_values)
n_r_trait

Proportion of calves number of offspring per category for male and female

pvec_proportion_progeny_adult <- 0.5
pvec_number_progeny <- c(20,5)

proportion_calves <- 1 -  pvec_proportion_progeny_adult
proportion_calves

pvec_proportion_selected <- c(0.1,0.5)
pvec_generation_interval <- c(5.7,5.4)

number of offspring per category for male and female

male_number_adults <- floor(pvec_proportion_progeny_adult*pvec_number_progeny[[1]])
male_number_adults
male_number_calves <- pvec_number_progeny[[1]] - male_number_adults
male_number_calves

-----------------------------------------

Starting with male selection candidate

Design matrix Z

  ### ## Building Design Matrix Z
  ### ## Design Matrix for one calf
  male_designMatrix_calf <- cbind(c(1,0,0),c(0,0,0),c(0,1,0),c(0,0,0),c(0,0,1),c(0,0,0))
  male_designMatrix_calf
  ### ## Design Matrix for one adult
  male_designMatrix_adult <- cbind(c(0,0,0),c(1,0,0),c(0,0,0),c(0,1,0),c(0,0,0),c(0,0,1))
  male_designMatrix_adult
  ### ## Design Matrix for calves
  male_designMatrix_calves <- diag(1,male_number_calves)%x%male_designMatrix_calf
  dim(male_designMatrix_calves)
  male_designMatrix_calves <- cbind(male_designMatrix_calves,
                                    matrix(0, nrow=nrow(male_designMatrix_calves),
                                              ncol=male_number_adults*n_r_trait))
  dim(male_designMatrix_calves)
  ### ## Design Matrix for adults
  male_designMatrix_adults <- diag(1,male_number_adults)%x%male_designMatrix_adult
  dim(male_designMatrix_adults)
  male_designMatrix_adults <- cbind(matrix(0, nrow=nrow(male_designMatrix_adults),
                                           ncol=male_number_calves*n_r_trait),
                                    male_designMatrix_adults)
  dim(male_designMatrix_adults)
  ### ## Bind Design Matrix for calves and adults
  male_designMatrixZ <- rbind(male_designMatrix_calves, male_designMatrix_adults)
  ### ## Add selection candidate
  male_designMatrixZ<-cbind(matrix(0, nrow=pvec_number_progeny[[1]]*(n_r_trait/2),ncol=n_r_trait),male_designMatrixZ)

Residual Matrix

## residual
s_residual_var_cov <- "../inst/extdata/variance_component_values/adults_calves_residual_variances_covariances"
df_residual_var_cov <- read.table(file = s_residual_var_cov, sep=";", header=TRUE)
mat_residual_var_cov <- as.matrix(df_residual_var_cov)
pmat_residual_varcov <- mat_residual_var_cov
pmat_residual_varcov

mat_residual_var_cov_calves <- pmat_residual_varcov[c(1,3,5),c(1,3,5)]
mat_residual_var_cov_calves
mat_residual_var_cov_adults <- pmat_residual_varcov[c(2,4,6),c(2,4,6)]
mat_residual_var_cov_adults

male_calves_kronecker_residual <- diag(male_number_calves) %x% mat_residual_var_cov_calves
male_adults_kronecker_residual <- diag(male_number_adults) %x% mat_residual_var_cov_adults
male_calves_kronecker_residual_extended <- cbind(male_calves_kronecker_residual,
                                                 matrix(0,nrow=nrow(male_calves_kronecker_residual),
                                                          ncol=ncol(male_adults_kronecker_residual)))
male_adults_kronecker_residual_extended <- cbind(matrix(0,nrow=nrow(male_adults_kronecker_residual),
                                                           ncol=ncol(male_calves_kronecker_residual)),
                                                male_adults_kronecker_residual)
male_ResidualMatrix <- rbind(male_calves_kronecker_residual_extended,male_adults_kronecker_residual_extended)

Relationship Matrix $A_y$ and Inverse $A_y^{-1}$

  mnumb <- pvec_number_progeny[[1]]+1
  male_Pedigree <- pedigreemm::pedigree(sire = c(NA,rep(1,times = pvec_number_progeny[[1]])),
                                        dam =c(NA,rep(NA,times = pvec_number_progeny[[1]])),
                                        label = 1:mnumb)
  ### ## Compute inverse of A
  male_AInv <- as.matrix(pedigreemm::getAInv(male_Pedigree))
  dim(male_AInv)

Genetic Variance Covariance Matrix G

s_genetic_var_cov <- "../inst/extdata/variance_component_values/adults_calves_genetic_variances_covariances"
df_genetic_var_cov <- read.table(file = s_genetic_var_cov, sep=";", header=TRUE)
pmat_genetic_varcov <- as.matrix(df_genetic_var_cov)
pmat_genetic_varcov

Computation of PEV

  male_PEV <- solve(t(male_designMatrixZ)%*%solve(male_ResidualMatrix)%*%male_designMatrixZ+male_AInv%x%solve(pmat_genetic_varcov))
  ### ## Computation of PEV for selection candidate
  male_PEV <- male_PEV[1:n_r_trait,1:n_r_trait]

Computation of variance covariance matrix for predicting breeding values

  mat_male_var_cov_pbv <- pmat_genetic_varcov - male_PEV

-----------------------------------------

Following with female selection candidate

Number of offspring per category for female

female_number_adults <- floor(pvec_proportion_progeny_adult*pvec_number_progeny[[2]])
female_number_calves <- pvec_number_progeny[[2]] - female_number_adults

Building Design Matrix Z

  ### ## Design Matrix for one calf
  female_designMatrix_calf <- cbind(c(1,0,0),c(0,0,0),c(0,1,0),c(0,0,0),c(0,0,1),c(0,0,0))
  ### ## Design Matrix for one adult
  female_designMatrix_adult <- cbind(c(0,0,0),c(1,0,0),c(0,0,0),c(0,1,0),c(0,0,0),c(0,0,1))
  ### ## Design Matrix for calves
  female_designMatrix_calves <- diag(1,female_number_calves)%x%female_designMatrix_calf
  female_designMatrix_calves <- cbind(female_designMatrix_calves,
                                    matrix(0, nrow=nrow(female_designMatrix_calves),
                                           ncol=female_number_adults*n_r_trait))
  ### ## Design Matrix for adults
  female_designMatrix_adults <- diag(1,female_number_adults)%x%female_designMatrix_adult
  female_designMatrix_adults <- cbind(matrix(0, nrow=nrow(female_designMatrix_adults),
                                           ncol=female_number_calves*n_r_trait),
                                    female_designMatrix_adults)
  ### ## Bind Design Matrix for calves and adults
  female_designMatrixZ <- rbind(female_designMatrix_calves, female_designMatrix_adults)
  ### ## Add selection candidate
  female_designMatrixZ<-cbind(matrix(0, nrow=pvec_number_progeny[[2]]*(n_r_trait/2),ncol=n_r_trait),female_designMatrixZ)

Residual Matrix R

  female_calves_kronecker_residual <- diag(female_number_calves) %x% mat_residual_var_cov_calves
  female_adults_kronecker_residual <- diag(female_number_adults) %x% mat_residual_var_cov_adults
  female_calves_kronecker_residual_extended <- cbind(female_calves_kronecker_residual,
                                                   matrix(0,nrow=nrow(female_calves_kronecker_residual),
                                                          ncol=ncol(female_adults_kronecker_residual)))
  female_adults_kronecker_residual_extended <- cbind(matrix(0,nrow=nrow(female_adults_kronecker_residual),
                                                          ncol=ncol(female_calves_kronecker_residual)),
                                                   female_adults_kronecker_residual)
  female_ResidualMatrix <- rbind(female_calves_kronecker_residual_extended,female_adults_kronecker_residual_extended)

Relationship Matrix A

  fnumb <- pvec_number_progeny[[2]]+1
  female_Pedigree <- pedigreemm::pedigree(sire = c(NA,rep(NA,times = pvec_number_progeny[[2]])),
                                          dam = c(NA,rep(1,times = pvec_number_progeny[[2]])),
                                         label = 1:fnumb)
  ### ## Compute inverse of A
  female_AInv <- as.matrix(pedigreemm::getAInv(female_Pedigree))

Computation of PEV

  female_PEV <- solve(t(female_designMatrixZ)%*%solve(female_ResidualMatrix)%*%female_designMatrixZ+female_AInv%x%solve(pmat_genetic_varcov))

  ### ## Computation of PEV for selection candidate
  female_PEV <- female_PEV[1:n_r_trait,1:n_r_trait]

Computation of variance covariance matrix for predicting breeding values

  mat_female_var_cov_pbv <- pmat_genetic_varcov - female_PEV

-----------------------------------------

Selection intensity and Factor of variance reduction k

  male_i <- dnorm(qnorm(1-pvec_proportion_selected[[1]]))/pvec_proportion_selected[[1]]
  male_x <- qnorm(pvec_proportion_selected[[1]], lower.tail = FALSE)
  male_k <- male_i*(male_i-male_x)

  female_i <- dnorm(qnorm(1-pvec_proportion_selected[[2]]))/pvec_proportion_selected[[2]]
  female_x <- qnorm(pvec_proportion_selected[[2]], lower.tail = FALSE)
  female_k <- female_i*(female_i-female_x)

Computation of asymptotic variance covariance matrix

  mat_kInv <- solve(matrix(c(1+0.5*male_k, 0.5*female_k,0.5*male_k, 1+0.5*female_k),nrow = 2, ncol = 2, byrow=TRUE))
  mat_var_fact_reduction <- mat_kInv%x%diag(1,n_r_trait)

  mat_var_cov_pbv <- rbind(mat_male_var_cov_pbv, mat_female_var_cov_pbv)
  mat_asympt_var_cov_pbv <- mat_var_fact_reduction%*%mat_var_cov_pbv

-----------------------------------------

Asymptotic genetic gain of the estimated overall breeding value (SZENARIO 2)

  ### ## Separate asymptotic variance covariance matrix for male and female
  mat_male_asympt_var_cov_pbv <- mat_asympt_var_cov_pbv[1:n_r_trait,]
  mat_female_asympt_var_cov_pbv <- mat_asympt_var_cov_pbv[(n_r_trait+1):nrow(mat_asympt_var_cov_pbv),]

  result_asym_sel_res <- (male_i*sqrt(t(as.matrix(pvec_economic_values))%*%mat_male_asympt_var_cov_pbv%*%as.matrix(pvec_economic_values))+
                          female_i*sqrt(t(as.matrix(pvec_economic_values))%*%mat_female_asympt_var_cov_pbv%*%as.matrix(pvec_economic_values)))/
                          sum(pvec_generation_interval)

result_asym_sel_res

Comparison with output function asym_sel_res.R

(res_LM_asym_select_response <- MeatValueIndex::compute_asym_selection_response(pvec_economic_values = pvec_economic_values,
                                            pvec_proportion_selected = pvec_proportion_selected,
                                            pvec_generation_interval = pvec_generation_interval,
                                            pvec_number_progeny = pvec_number_progeny,
                                            pvec_proportion_progeny_adult = pvec_proportion_progeny_adult,
                                            pmat_genetic_varcov = pmat_genetic_varcov,
                                            pmat_residual_varcov = pmat_residual_varcov,
                                            pb_out = FALSE))

OB - Szenario B

s_a_OB_B <- "../inst/extdata/economic_values/OB_economic_values_SzenarioB"

df_a_OB_B <- read.table(file = s_a_OB_B, sep=";", header=TRUE)
mat_a_OB_B <- as.matrix(df_a_OB_B)
pvec_economic_values_OB_szenarioB <- as.vector(mat_a_OB_B)

(res_OB_szenarioB_asym_select_response <- MeatValueIndex::compute_asym_selection_response(pvec_economic_values = pvec_economic_values_OB_szenarioB,
                                            pvec_proportion_selected = pvec_proportion_selected,
                                            pvec_generation_interval = pvec_generation_interval,
                                            pvec_number_progeny = pvec_number_progeny,
                                            pvec_proportion_progeny_adult = pvec_proportion_progeny_adult,
                                            pmat_genetic_varcov = pmat_genetic_varcov,
                                            pmat_residual_varcov = pmat_residual_varcov,
                                            pb_out = FALSE))

OB - Szenario B*

s_a_OB_BStern <- "../inst/extdata/economic_values/OB_economic_values_SzenarioB*"

df_a_OB_BStern <- read.table(file = s_a_OB_BStern, sep=";", header=TRUE)
mat_a_OB_BStern <- as.matrix(df_a_OB_BStern)
pvec_economic_values_OB_szenarioBStern <- as.vector(mat_a_OB_BStern)

(res_OB_szenarioB_asym_select_response <- MeatValueIndex::compute_asym_selection_response(pvec_economic_values = pvec_economic_values_OB_szenarioBStern,
                                            pvec_proportion_selected = pvec_proportion_selected,
                                            pvec_generation_interval = pvec_generation_interval,
                                            pvec_number_progeny = pvec_number_progeny,
                                            pvec_proportion_progeny_adult = pvec_proportion_progeny_adult,
                                            pmat_genetic_varcov = pmat_genetic_varcov,
                                            pmat_residual_varcov = pmat_residual_varcov,
                                            pb_out = FALSE))

pvrqualitasag/MeatValueIndex documentation built on May 13, 2019, 4:44 p.m.