In pvrqualitasag/rvcetools: Post-Processing Of Variance Components Results

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(rvcetools)

Disclaimer

The package rvcetools is used to post-process and to work with results from programs that do variance components estimation (VCE). Initially, we limit our focus for result files from the program called vce.

Package Usage

One possible use case for this package is to produce an input parameter file for the prediction of breeding values. Because we are using MiX99 to predict breeding values, the next subsection describes how to prepare a parameter input file for MiX99.

Input For MiX99

A parameter input file for MiX99 can be produced using the function parameter_varCovar_mix99(). The function requires at least an input file with the results from the variance components estimation and an output file in which the parameters for the prediction of the breeding values are stored. Based on that a minimal call looks as follows

sInputFile <- system.file("extdata","VCE_results.csv", package = "rvcetools")
sMix99File <- file.path(tempdir(), 'mix99_output.txt')
parameter_varCovar_mix99(ps_input_file = sInputFile, ps_output_file = sMix99File)

The above function call reads the variance-covariance matrices, checks whether they are positive definite or not. If they are not, the matrices are bent using the method described in the book by L. Schaeffer in chapter 7.4. The produced parameter file has the following structure.

vec_result <- readLines(con = file(sMix99File))
file.remove(sMix99File)
head(vec_result)

There are two more bending methods available in rvcetools.

1. bending based on an upper bound of the ratio of the largest to the smallest eigenvalue and
2. weighted and unweighted bending described in a paper by Jorjani et al. 2003

These two methods can be chosen by changing the parameters of the call to parameter_varCovar_mix99(). The eigenvalue ratio based method is applied with the following function call

parameter_varCovar_mix99(ps_input_file = sInputFile, ps_output_file = sMix99File, pn_ratio = 100)

The additional parameter pn_ratio specifies the maximal allowed ratio of the largest to the smallest eigenvalue. The structure of the output file is the same

vec_result <- readLines(con = file(sMix99File))
file.remove(sMix99File)
head(vec_result)

Unweighted bending is used when specifying the parameter pn_eps which can be used to give a fixed value to the smallest eigenvalue of the resulting matrix.

parameter_varCovar_mix99(ps_input_file = sInputFile, ps_output_file = sMix99File, pn_eps = 1e-4)

file.remove(sMix99File)

The weighted version of the method can be used when specifying a weight matrix using the function parameter pmat_weight. The specified weight matrix must have the same dimension as the variance-covariance matrices. The current version allows only to specify one weight matrix for all variance-covariance matrices of the different random effects.

Background

Variance components estimation is the process of estimating scale parameters such as variance components from data. The parameter estimation uses linear mixed models to define the set of different random effects to which the observed variation should be attributed to.

Linear mixed effects models recently have gained some popularity outside of the area of animal breeding. But still there are not standard software packages around that can estimate variance components for very large datasets. Therefore, specialized softare programs are used for this task. These programs do not provide any features outside of the parameter estimation functionality. As a consequence of that results from the specialised programs must be read into systems such as R to do bending, plotting or other post-processing tasks.

Features

The most central feature is to read in raw outfiles from the different variance components estimation programs. In the most basic case there is only one output file from one analysis. In more advanced analyses, there might be output files from repeated analyses of different sample data sets.

The package rvcetools contains a number of features already and they will hopefully grow over time. For the purpose of a genetic evaluation the most important feature is to read estimated variance components and to produce an input file for the software program that does the prediction of breeding values. This function is called parameter_varCovar_mix99(). All it needs is an input file and an output file and all steps in between will be done automatically. The functions envolved in all the steps are shown in the dependency graph below.

if (require(CodeDepends)){
  gg <- CodeDepends::makeCallGraph("package:rvcetools")
  if(require(Rgraphviz)) {
    gg = layoutGraph(gg, layoutType = "circo")
    graph.par(list(nodes = list(fontsize=55)))
    renderGraph(gg) ## could also call plot directly
  } 
}

Conversion of Correlation to Covariance

When it comes to processing outputs of variance components estimation, it is useful to transform variance-covariance matrices into correlation matrices and the other way around. In most practical calses the matrices are small and hence a solution via iterative loops should be ok. An example is shown here

(mat_vcov <- matrix(c(104,75,18,75,56,12,18,12,7), nrow = 3, byrow = TRUE))

In base-R there is the function cov2cor which does this.

cov2cor(mat_vcov)

The package contains a wrapper that does the same thing

(mat_cor <- cov_to_cor(mat_vcov))

In rare cases there might also be the interest of going the other way round. Hence given a correlation matrix and a vector of variances, we might be interested in re-building the original variance-covariance matrix. This is done with

cor_to_cov(pmat_cor = mat_cor, pvec_var = diag(mat_vcov))

Session Info

sessionInfo()

Latest Update

r paste(Sys.time(),paste0("(", Sys.info()[["user"]],")" ))

pvrqualitasag/rvcetools documentation built on Dec. 31, 2021, 2:09 p.m.