Methods to reduce the dimension of Multivariate Models

1. Canonical Transformation

Canonical transformations involve the use of special matrices in the transformation of related traits into unrelated traits and thus using single trait models for evaluation( transformation of polymorphic model into monomorphic models). Using the transformation matrix Q and P we are able to transform the models.
Limitations for practical use occur when the models are different and some values missing.
A methodology by Ducrocq and Besbes(1993):
An example (R.A Mrode 2014 Linear Models for the Prediction of Animal Breeding values Exp 5.1 section 5.2.2.) to calculate the transformation Matrix and its inverse.

if (!("MASS" %in% installed.packages())) { install.packages("MASS", dependencies = TRUE) }
if (!("dplyr" %in% installed.packages())) { install.packages("plyr", dependencies = TRUE) }
require(MASS)

Loading required package: MASS

require(dplyr)

Loading required package: dplyr

Attaching package: 㤼㸱dplyr㤼㸲

The following object is masked from 㤼㸱package:MASS㤼㸲:

    select

The following objects are masked from 㤼㸱package:stats㤼㸲:

    filter, lag

The following objects are masked from 㤼㸱package:base㤼㸲:

    intersect, setdiff, setequal, union

consider the G and R matrices below:

G <- matrix(c(20, 18, 18, 40), nrow = 2)
G

     [,1] [,2]
[1,]   20   18
[2,]   18   40

R <- matrix(c(40, 11, 11, 30), nrow = 2)
R

     [,1] [,2]
[1,]   40   11
[2,]   11   30

1. First we calculate the eigenvalues (B) and eigenvectors (U) for
   the residual covariance matrix (R).
   R = UBU’
   B is the diagonal matrix of the eigenvalues of (R)
   U is the matrix of the eigenvectors of (R)

Eigen_B <- eigen(R)
Eigen_B

eigen() decomposition
$values
[1] 47.08305 22.91695

$vectors
           [,1]       [,2]
[1,] -0.8407743  0.5413857
[2,] -0.5413857 -0.8407743

B <- diag(Eigen_B$values) 
B

         [,1]     [,2]
[1,] 47.08305  0.00000
[2,]  0.00000 22.91695

U <- (Eigen_B$vectors)
U

           [,1]       [,2]
[1,] -0.8407743  0.5413857
[2,] -0.5413857 -0.8407743

2. Calculate P and PGP’

\(\mathbf{P} = \mathbf{U}\sqrt{B^{-1}}\mathbf{{U^\prime}}\)

P <- U %*% sqrt(ginv(B))%*%(t(U))
P

            [,1]        [,2]
[1,]  0.16424710 -0.02874736
[2,] -0.02874736  0.19038107

\(\mathbf{PGP^\prime}\)

PGP_prime <- P %*% G %*% t(P)
PGP_prime

          [,1]      [,2]
[1,] 0.4026185 0.2643755
[2,] 0.2643755 1.2692999

3. Then calculate the eigenvalues (W) and eignevectors (L) of PGP’
   \(\mathbf{PGP^\prime} = \mathbf{LWL^\prime}\)

Eigen_PGP <- eigen(PGP_prime)
Eigen_PGP

eigen() decomposition
$values
[1] 1.3435798 0.3283387

$vectors
          [,1]       [,2]
[1,] 0.2704897 -0.9627229
[2,] 0.9627229  0.2704897

W <- diag(Eigen_PGP$values)
W

        [,1]      [,2]
[1,] 1.34358 0.0000000
[2,] 0.00000 0.3283387

L <- Eigen_PGP$vectors
L

          [,1]       [,2]
[1,] 0.2704897 -0.9627229
[2,] 0.9627229  0.2704897

4. The transformation matrix Q and its inverse can be obtained as:
   \(\mathbf{Q = L^\prime P}\)

Q <- t(L) %*% P
Q

            [,1]       [,2]
[1,]  0.01675141 0.17550834
[2,] -0.16590031 0.07917187

Qinv <- ginv(Q)
Qinv

         [,1]       [,2]
[1,] 2.600648 -5.7651217
[2,] 5.449516  0.5502527

After obtaining Q and Q’ we can go ahead and transform example data using the transformation matrix Q.
Load example data

# set working directory
setwd("F:/R/RMME")
getwd()

[1] "F:/R/RMME"

# Load data
data <- read.csv("F:/R/RMME/WGdata.txt", sep="")
data

Transformaing the data gives:

# data number
data_no <- nrow(data)
data_no

[1] 5

# data transformation
for (i in 1:data_no){
  data[i,7:8]= Q %*% as.numeric(data[i,5:6])
}
colnames(data)

[1] "Calves" "Sex"    "Sire"   "Dam"    "WWG"    "PWG"    "V7"     "V8"    

names(data)[names(data) == "V7"] <- "tWWG"
names(data)[names(data) == "V8"] <- "tPWG"
data

# the pedigree file will be
ped1 <- data[,c(1,3,4)]
ped1

ped2 <- (cbind(Calves=seq(1:3),Sire=rep(0,3),Dam=rep(0,3)))
ped2

     Calves Sire Dam
[1,]      1    0   0
[2,]      2    0   0
[3,]      3    0   0

ped  <- rbind(ped2,ped1)
ped  <- ped[order(ped[,1]),]
colnames(ped) <- c("animal","sire","dam")
ped

Load pedigree file

Two functions will be required;
(i) To create design matrices (desgn)
(ii) To create the inverse of the Numerator Relationship Matrix(ainv)

Function (desgn)

desgn <- function(v)
  { 
  if (is.numeric(v)) { 
     va = v 
     mrow = length(va) 
     mcol = max(va) 
  } 
  if (is.character(v)) 
    { vf = factor(v) # save the index value
  va = as.numeric(vf) 
  mrow = length(va) 
  mcol = length(levels(vf))
  }

  X = matrix(data = c(0), nrow = mrow, ncol = mcol) 
  for (i in 1:mrow) { 
    ic = va[i] 
    X[i, ic] = 1 } 
  return(X) 
 }

Function(ainv)

ainv = function(ped) {
  n = nrow(ped)
  Ainv = matrix(c(0), nrow = n, ncol = n)
  
  for (i in 1:n) {
    animal = ped[i, 1]
    sire = ped[i, 2]
    dam = ped[i, 3]
    
    if (sire == 0 & dam == 0) {
      # both parents unknown
      alpha = 1
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
    } else if (sire != 0 & dam == 0) {
      # sire known
      alpha = 4/3
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
      Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
      Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
    } else if (sire == 0 & dam != 0) {
      # dam known
      alpha = 4/3
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
      Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
      Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
    } else {
      # both parents known
      alpha = 2
      Ainv[animal, animal] = alpha + Ainv[animal, animal]
      Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
      Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
      Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
      Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
      Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
      Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
      Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
      Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
    }
  }
  return(Ainv)
}

Setting up the MME (mixed model equations)

Variance components from G and R matrices give alpha_1 and alpha_2 for WWG and PWG
respectively.(see page 73 example 5.2.2) After transformation \(\mathbf{\alpha_i=1/w_{ii}}\)

alpha_1 = 1/0.3283387
alpha_1

[1] 3.045635

alpha_2 = 1/1.34358
alpha_2  

[1] 0.7442802

The NRM matrix and its inverse (A)

A <- ainv(ped)
A

           [,1] [,2] [,3]       [,4] [,5] [,6] [,7] [,8]
[1,]  1.8333333  0.5  0.0 -0.6666667  0.0 -1.0    0    0
[2,]  0.5000000  2.0  0.5  0.0000000 -1.0 -1.0    0    0
[3,]  0.0000000  0.5  2.0  0.0000000 -1.0  0.5    0   -1
[4,] -0.6666667  0.0  0.0  1.8333333  0.5  0.0   -1    0
[5,]  0.0000000 -1.0 -1.0  0.5000000  2.5  0.0   -1    0
[6,] -1.0000000 -1.0  0.5  0.0000000  0.0  2.5    0   -1
[7,]  0.0000000  0.0  0.0 -1.0000000 -1.0  0.0    2    0
[8,]  0.0000000  0.0 -1.0  0.0000000  0.0 -1.0    0    2

Formulation fo the design matrices
- design matrix for the sex effect

X = desgn(data[,2])
X

     [,1] [,2]
[1,]    1    0
[2,]    0    1
[3,]    0    1
[4,]    1    0
[5,]    1    0

t(X)

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    1    1
[2,]    0    1    1    0    0

The number of levels for the sex effect will be:

no_levels_sex = ncol(X)
no_levels_sex

[1] 2

• design matrix for the animal effect

Z = desgn(data[,1])
Z

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    0    0    0    1    0    0    0    0
[2,]    0    0    0    0    1    0    0    0
[3,]    0    0    0    0    0    1    0    0
[4,]    0    0    0    0    0    0    1    0
[5,]    0    0    0    0    0    0    0    1

no_animals = ncol(Z)
no_animals

[1] 8

The matrices of observations(y1 and y2 for WWG and PWG respectively) are given by:

y1 = data[,7]
y1

[1] 1.2688381 0.9261208 1.2587872 1.1116800 1.4000696

y2 = data[,8]
y2

[1] -0.20818267 -0.08525154 -0.10864249 -0.10561986 -0.23571251

Construction of the Left Hand Side (LHS)
Lets first work on WWG

LHS1 = rbind(
  cbind(t(X) %*% X, t(X) %*% Z), 
  cbind(t(Z) %*% X, t(Z) %*% Z + A * alpha_2))
LHS1

      [,1] [,2]       [,3]       [,4]       [,5]       [,6]       [,7]       [,8]       [,9]      [,10]
 [1,]    3    0  0.0000000  0.0000000  0.0000000  1.0000000  0.0000000  0.0000000  1.0000000  1.0000000
 [2,]    0    2  0.0000000  0.0000000  0.0000000  0.0000000  1.0000000  1.0000000  0.0000000  0.0000000
 [3,]    0    0  1.3645137  0.3721401  0.0000000 -0.4961868  0.0000000 -0.7442802  0.0000000  0.0000000
 [4,]    0    0  0.3721401  1.4885604  0.3721401  0.0000000 -0.7442802 -0.7442802  0.0000000  0.0000000
 [5,]    0    0  0.0000000  0.3721401  1.4885604  0.0000000 -0.7442802  0.3721401  0.0000000 -0.7442802
 [6,]    1    0 -0.4961868  0.0000000  0.0000000  2.3645137  0.3721401  0.0000000 -0.7442802  0.0000000
 [7,]    0    1  0.0000000 -0.7442802 -0.7442802  0.3721401  2.8607005  0.0000000 -0.7442802  0.0000000
 [8,]    0    1 -0.7442802 -0.7442802  0.3721401  0.0000000  0.0000000  2.8607005  0.0000000 -0.7442802
 [9,]    1    0  0.0000000  0.0000000  0.0000000 -0.7442802 -0.7442802  0.0000000  2.4885604  0.0000000
[10,]    1    0  0.0000000  0.0000000 -0.7442802  0.0000000  0.0000000 -0.7442802  0.0000000  2.4885604

for PWG

LHS2 = rbind(
  cbind(t(X) %*% X, t(X) %*% Z), 
  cbind(t(Z) %*% X, t(Z) %*% Z + A * alpha_1))
LHS2

      [,1] [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]      [,9]     [,10]
 [1,]    3    0  0.000000  0.000000  0.000000  1.000000  0.000000  0.000000  1.000000  1.000000
 [2,]    0    2  0.000000  0.000000  0.000000  0.000000  1.000000  1.000000  0.000000  0.000000
 [3,]    0    0  5.583665  1.522818  0.000000 -2.030424  0.000000 -3.045635  0.000000  0.000000
 [4,]    0    0  1.522818  6.091271  1.522818  0.000000 -3.045635 -3.045635  0.000000  0.000000
 [5,]    0    0  0.000000  1.522818  6.091271  0.000000 -3.045635  1.522818  0.000000 -3.045635
 [6,]    1    0 -2.030424  0.000000  0.000000  6.583665  1.522818  0.000000 -3.045635  0.000000
 [7,]    0    1  0.000000 -3.045635 -3.045635  1.522818  8.614089  0.000000 -3.045635  0.000000
 [8,]    0    1 -3.045635 -3.045635  1.522818  0.000000  0.000000  8.614089  0.000000 -3.045635
 [9,]    1    0  0.000000  0.000000  0.000000 -3.045635 -3.045635  0.000000  7.091271  0.000000
[10,]    1    0  0.000000  0.000000 -3.045635  0.000000  0.000000 -3.045635  0.000000  7.091271

RHS construction for WWG

RHS1 = rbind(t(X) %*% y1, t(Z) %*% y1)
RHS1

           [,1]
 [1,] 3.7805877
 [2,] 2.1849080
 [3,] 0.0000000
 [4,] 0.0000000
 [5,] 0.0000000
 [6,] 1.2688381
 [7,] 0.9261208
 [8,] 1.2587872
 [9,] 1.1116800
[10,] 1.4000696

RHS construction for PWG

RHS2 = rbind(t(X) %*% y2, t(Z) %*% y2)
RHS2

             [,1]
 [1,] -0.54951504
 [2,] -0.19389403
 [3,]  0.00000000
 [4,]  0.00000000
 [5,]  0.00000000
 [6,] -0.20818267
 [7,] -0.08525154
 [8,] -0.10864249
 [9,] -0.10561986
[10,] -0.23571251

The solution for the MMEs will be

# generalized inverse of LHS for WWG
gi_LHS1 = ginv(LHS1)
gi_LHS1

            [,1]       [,2]        [,3]        [,4]        [,5]       [,6]       [,7]       [,8]       [,9]      [,10]
 [1,]  1.0339157  0.4232500 -0.43499518 -0.22362124 -0.36351988 -0.6996631 -0.4036403 -0.4428597 -0.7454441 -0.6566399
 [2,]  0.4232500  1.2887131 -0.37182316 -0.62660706 -0.27660090 -0.2887349 -0.7749368 -0.8024894 -0.4882018 -0.4928132
 [3,] -0.4349952 -0.3718232  1.20572204  0.03132429  0.15627342  0.5282514  0.2132414  0.5304050  0.3965640  0.3801701
 [4,] -0.2236212 -0.6266071  0.03132429  1.30211999 -0.04621605  0.1001898  0.6115194  0.6416947  0.3027181  0.2679558
 [5,] -0.3635199 -0.2766009  0.15627342 -0.04621605  1.13499436  0.2259291  0.4404279  0.1127738  0.3453709  0.5192597
 [6,] -0.6996631 -0.2887349  0.52825138  0.10018979  0.22592911  0.9979403  0.2242076  0.3532623  0.6466725  0.4543767
 [7,] -0.4036403 -0.7749368  0.21324136  0.61151944  0.44042795  0.2242076  1.0026430  0.5472305  0.5291256  0.4575878
 [8,] -0.4428597 -0.8024894  0.53040496  0.64169467  0.11277384  0.3532623  0.5472305  1.0577483  0.4472781  0.5280387
 [9,] -0.7454441 -0.4882018  0.39656403  0.30271812  0.34537086  0.6466725  0.5291256  0.4472781  1.0530455  0.5366143
[10,] -0.6566399 -0.4928132  0.38017012  0.26795580  0.51925967  0.4543767  0.4575878  0.5280387  0.5366143  0.9789287

# solution
sol1 = gi_LHS1 %*% RHS1
sol1

              [,1]
 [1,]  1.266489432
 [2,]  1.088949914
 [3,]  0.051599829
 [4,] -0.001593477
 [5,] -0.031209521
 [6,] -0.002395332
 [7,] -0.088391637
 [8,]  0.095399837
 [9,] -0.089361056
[10,]  0.072875775

# generalized inverse of LHS for PWG
gi_LHS2 = ginv(LHS2)
gi_LHS2

             [,1]        [,2]         [,3]         [,4]         [,5]         [,6]        [,7]        [,8]        [,9]       [,10]
 [1,]  0.50583878  0.10315741 -0.108239013 -0.054909474 -0.084752712 -0.174775757 -0.09693905 -0.10937576 -0.18803160 -0.15470898
 [2,]  0.10315741  0.70045727 -0.085914006 -0.160046511 -0.075729022 -0.052930744 -0.19866943 -0.20224510 -0.12260705 -0.13393442
 [3,] -0.10823901 -0.08591401  0.314411868  0.003383116  0.015711692  0.149453970  0.02168938  0.15013864  0.08876819  0.08649488
 [4,] -0.05490947 -0.16004651  0.003383116  0.324631479 -0.005045589  0.009746099  0.15807230  0.16202072  0.07981970  0.07516262
 [5,] -0.08475271 -0.07572902  0.015711692 -0.005045589  0.307386755  0.023284610  0.14069914  0.01075890  0.08238121  0.14859232
 [6,] -0.17477576 -0.05293074  0.149453970  0.009746099  0.023284610  0.293508017  0.02266150  0.08319998  0.16043850  0.07038076
 [7,] -0.09693905 -0.19866943  0.021689376  0.158072301  0.140699144  0.022661505  0.29364653  0.10369234  0.14952158  0.11863408
 [8,] -0.10937576 -0.20224510  0.150138636  0.162020720  0.010758901  0.083199983  0.10369234  0.30079787  0.09569253  0.14923476
 [9,] -0.18803160 -0.12260705  0.088768186  0.079819704  0.082381206  0.160438499  0.14952158  0.09569253  0.30065936  0.10299695
[10,] -0.15470898 -0.13393442  0.086494884  0.075162620  0.148592321  0.070380756  0.11863408  0.14923476  0.10299695  0.29074925

# solution
sol2 = gi_LHS2 %*% RHS2
sol2

               [,1]
 [1,] -0.1851085678
 [2,] -0.0980527192
 [3,] -0.0029006449
 [4,]  0.0019511186
 [5,] -0.0004810105
 [6,]  0.0006954827
 [7,]  0.0070173543
 [8,] -0.0048059447
 [9,]  0.0145219665
[10,] -0.0094067883

# combined solution 
sol = cbind(sol1,sol2)
colnames(sol) <- c("tsol_WWG","tsol_PWG")
sol

          tsol_WWG      tsol_PWG
 [1,]  1.266489432 -0.1851085678
 [2,]  1.088949914 -0.0980527192
 [3,]  0.051599829 -0.0029006449
 [4,] -0.001593477  0.0019511186
 [5,] -0.031209521 -0.0004810105
 [6,] -0.002395332  0.0006954827
 [7,] -0.088391637  0.0070173543
 [8,]  0.095399837 -0.0048059447
 [9,] -0.089361056  0.0145219665
[10,]  0.072875775 -0.0094067883

nsol <- nrow(sol)
nsol

[1] 10

After obtaining the canonical solutions we need to transform then using \(\mathbf{Q^{-1}}\) to ordinary solutions as below:

sol_combined =cbind(sol,sol_WWG=rep(0,nsol),sol_PWG=rep(0,nsol))
                    
for(i in 1: nsol){
  sol_combined[i,3:4] = Qinv %*% as.numeric(sol[i,1:2])
}
sol_combined

          tsol_WWG      tsol_PWG     sol_WWG      sol_PWG
 [1,]  1.266489432 -0.1851085678  4.36086700  6.799897623
 [2,]  1.088949914 -0.0980527192  3.39726159  5.880295936
 [3,]  0.051599829 -0.0029006449  0.15091558  0.279597996
 [4,] -0.001593477  0.0019511186 -0.01539251 -0.007610072
 [5,] -0.031209521 -0.0004810105 -0.07839190 -0.170341453
 [6,] -0.002395332  0.0006954827 -0.01023896 -0.012670709
 [7,] -0.088391637  0.0070173543 -0.27033146 -0.477830301
 [8,]  0.095399837 -0.0048059447  0.27580828  0.517238429
 [9,] -0.089361056  0.0145219665 -0.31611758 -0.478983734
[10,]  0.072875775 -0.0094067883  0.24375554  0.391961574

round(sol_combined,3)

      tsol_WWG tsol_PWG sol_WWG sol_PWG
 [1,]    1.266   -0.185   4.361   6.800
 [2,]    1.089   -0.098   3.397   5.880
 [3,]    0.052   -0.003   0.151   0.280
 [4,]   -0.002    0.002  -0.015  -0.008
 [5,]   -0.031    0.000  -0.078  -0.170
 [6,]   -0.002    0.001  -0.010  -0.013
 [7,]   -0.088    0.007  -0.270  -0.478
 [8,]    0.095   -0.005   0.276   0.517
 [9,]   -0.089    0.015  -0.316  -0.479
[10,]    0.073   -0.009   0.244   0.392

# Solution for the Male sex effect
Male_sol <- sol_combined[1,]
Male_sol

  tsol_WWG   tsol_PWG    sol_WWG    sol_PWG 
 1.2664894 -0.1851086  4.3608670  6.7998976 

# Solution for the female sex effect
Female_sol <- sol_combined[2,]
Female_sol

   tsol_WWG    tsol_PWG     sol_WWG     sol_PWG 
 1.08894991 -0.09805272  3.39726159  5.88029594 

# Solution for the animal effect
Animal_sol <- sol_combined[3:10]
Animal_sol

[1]  0.051599829 -0.001593477 -0.031209521 -0.002395332 -0.088391637  0.095399837 -0.089361056  0.072875775

NOTE: It should be noted that in the Mrode book, W ( the eigenvectors were interchanged) was interchanged and thus in order to show similar results i interchanged alpha1 and alpha2 from W.