Common Environment Effect Model

An example from the Linear Models for the Prediction of Animal Breeding Values, 3rd Edition ( Example 4.3 p67)

The matrix notation of the models is:
y = Xb + Za + Wc + e
where c is the common environment effect and W relates the records to the common environment effects.

Start by installing or loading the MASS packages if they are already installed.

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

require(MASS)

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)
}

Set the working directory

setwd("F:/R/RMME")

Print working directory

getwd()

[1] "F:/R/RMME"

The next step is to load the pedgree and data files. Read the pedigree file and sort by the animal column as follows;

ped <- read.csv("F:/R/RMME/ped.txt", sep="")
ped <- ped[order(ped[,1]),]

Read the data file

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

The variance components are given below:

sigma_2_a = 20
sigma_2_e = 65
sigma_2_c = 15
sigma_2_y = 100
alpha_1   = sigma_2_e/sigma_2_a
alpha_2   = sigma_2_e/sigma_2_c
alpha_1

[1] 3.25

alpha_2

[1] 4.333333

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

X = desgn(data[,4])
X

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

t(X)

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    0    0    0    1    0    0    1    0     1
[2,]    0    1    1    1    0    1    1    0    1     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] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
 [1,]    0    0    0    0    0    1    0    0    0     0     0     0     0     0     0
 [2,]    0    0    0    0    0    0    1    0    0     0     0     0     0     0     0
 [3,]    0    0    0    0    0    0    0    1    0     0     0     0     0     0     0
 [4,]    0    0    0    0    0    0    0    0    1     0     0     0     0     0     0
 [5,]    0    0    0    0    0    0    0    0    0     1     0     0     0     0     0
 [6,]    0    0    0    0    0    0    0    0    0     0     1     0     0     0     0
 [7,]    0    0    0    0    0    0    0    0    0     0     0     1     0     0     0
 [8,]    0    0    0    0    0    0    0    0    0     0     0     0     1     0     0
 [9,]    0    0    0    0    0    0    0    0    0     0     0     0     0     1     0
[10,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     1

The number of animals is:

no_animals = ncol(Z)
no_animals

[1] 15

• Design matrix for the common environmental effect. This is unique to the dams with offspring because the common environment effects are calculated for full-siblings or offspring of the same dam.

W = desgn(data[,3])[,unique(data[,3])]
W

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

The number of animals (levels) will be:

no_animals_ce = ncol(W)
no_animals_ce

[1] 3

The matrix of observations(y) is given by:

y = data[,5]
y

 [1]  90  70  65  98 106  60  80 100  85  68

The inverse of the numerator relationship matrix is obtained from the function ainv as:

A = ainv(ped)
A

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
 [1,]  4.0  1.5    0    0  1.5   -1   -1   -1    0     0     0     0    -1    -1    -1
 [2,]  1.5  2.5    0    0  0.0   -1   -1   -1    0     0     0     0     0     0     0
 [3,]  0.0  0.0    3    2  0.0    0    0    0   -1    -1    -1    -1     0     0     0
 [4,]  0.0  0.0    2    3  0.0    0    0    0   -1    -1    -1    -1     0     0     0
 [5,]  1.5  0.0    0    0  2.5    0    0    0    0     0     0     0    -1    -1    -1
 [6,] -1.0 -1.0    0    0  0.0    2    0    0    0     0     0     0     0     0     0
 [7,] -1.0 -1.0    0    0  0.0    0    2    0    0     0     0     0     0     0     0
 [8,] -1.0 -1.0    0    0  0.0    0    0    2    0     0     0     0     0     0     0
 [9,]  0.0  0.0   -1   -1  0.0    0    0    0    2     0     0     0     0     0     0
[10,]  0.0  0.0   -1   -1  0.0    0    0    0    0     2     0     0     0     0     0
[11,]  0.0  0.0   -1   -1  0.0    0    0    0    0     0     2     0     0     0     0
[12,]  0.0  0.0   -1   -1  0.0    0    0    0    0     0     0     2     0     0     0
[13,] -1.0  0.0    0    0 -1.0    0    0    0    0     0     0     0     2     0     0
[14,] -1.0  0.0    0    0 -1.0    0    0    0    0     0     0     0     0     2     0
[15,] -1.0  0.0    0    0 -1.0    0    0    0    0     0     0     0     0     0     2

Construction of the Left Hand Side (LHS)

LHS = rbind(
  cbind(t(X) %*% X, t(X) %*% Z,  t(X) %*% W), 
  cbind(t(Z) %*% X, t(Z) %*% Z + A * alpha_1, t(Z) %*% W),
  cbind(t(W) %*% X, t(W) %*% Z, t(W) %*% W + diag(no_animals_ce) * alpha_2))
LHS

      [,1] [,2]   [,3]   [,4]  [,5]  [,6]   [,7]  [,8]  [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17]    [,18]    [,19]    [,20]
 [1,]    4    0  0.000  0.000  0.00  0.00  0.000  1.00  0.00  0.00  0.00  1.00  0.00  0.00  1.00  0.00  1.00 1.000000 1.000000 2.000000
 [2,]    0    6  0.000  0.000  0.00  0.00  0.000  0.00  1.00  1.00  1.00  0.00  1.00  1.00  0.00  1.00  0.00 2.000000 3.000000 1.000000
 [3,]    0    0 13.000  4.875  0.00  0.00  4.875 -3.25 -3.25 -3.25  0.00  0.00  0.00  0.00 -3.25 -3.25 -3.25 0.000000 0.000000 0.000000
 [4,]    0    0  4.875  8.125  0.00  0.00  0.000 -3.25 -3.25 -3.25  0.00  0.00  0.00  0.00  0.00  0.00  0.00 0.000000 0.000000 0.000000
 [5,]    0    0  0.000  0.000  9.75  6.50  0.000  0.00  0.00  0.00 -3.25 -3.25 -3.25 -3.25  0.00  0.00  0.00 0.000000 0.000000 0.000000
 [6,]    0    0  0.000  0.000  6.50  9.75  0.000  0.00  0.00  0.00 -3.25 -3.25 -3.25 -3.25  0.00  0.00  0.00 0.000000 0.000000 0.000000
 [7,]    0    0  4.875  0.000  0.00  0.00  8.125  0.00  0.00  0.00  0.00  0.00  0.00  0.00 -3.25 -3.25 -3.25 0.000000 0.000000 0.000000
 [8,]    1    0 -3.250 -3.250  0.00  0.00  0.000  7.50  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 1.000000 0.000000 0.000000
 [9,]    0    1 -3.250 -3.250  0.00  0.00  0.000  0.00  7.50  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 1.000000 0.000000 0.000000
[10,]    0    1 -3.250 -3.250  0.00  0.00  0.000  0.00  0.00  7.50  0.00  0.00  0.00  0.00  0.00  0.00  0.00 1.000000 0.000000 0.000000
[11,]    0    1  0.000  0.000 -3.25 -3.25  0.000  0.00  0.00  0.00  7.50  0.00  0.00  0.00  0.00  0.00  0.00 0.000000 1.000000 0.000000
[12,]    1    0  0.000  0.000 -3.25 -3.25  0.000  0.00  0.00  0.00  0.00  7.50  0.00  0.00  0.00  0.00  0.00 0.000000 1.000000 0.000000
[13,]    0    1  0.000  0.000 -3.25 -3.25  0.000  0.00  0.00  0.00  0.00  0.00  7.50  0.00  0.00  0.00  0.00 0.000000 1.000000 0.000000
[14,]    0    1  0.000  0.000 -3.25 -3.25  0.000  0.00  0.00  0.00  0.00  0.00  0.00  7.50  0.00  0.00  0.00 0.000000 1.000000 0.000000
[15,]    1    0 -3.250  0.000  0.00  0.00 -3.250  0.00  0.00  0.00  0.00  0.00  0.00  0.00  7.50  0.00  0.00 0.000000 0.000000 1.000000
[16,]    0    1 -3.250  0.000  0.00  0.00 -3.250  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  7.50  0.00 0.000000 0.000000 1.000000
[17,]    1    0 -3.250  0.000  0.00  0.00 -3.250  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  7.50 0.000000 0.000000 1.000000
[18,]    1    2  0.000  0.000  0.00  0.00  0.000  1.00  1.00  1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 7.333333 0.000000 0.000000
[19,]    1    3  0.000  0.000  0.00  0.00  0.000  0.00  0.00  0.00  1.00  1.00  1.00  1.00  0.00  0.00  0.00 0.000000 8.333333 0.000000
[20,]    2    1  0.000  0.000  0.00  0.00  0.000  0.00  0.00  0.00  0.00  0.00  0.00  0.00  1.00  1.00  1.00 0.000000 0.000000 7.333333

RHS construction

RHS = rbind(t(X) %*% y, t(Z) %*% y, t(W) %*% y)
RHS

      [,1]
 [1,]  364
 [2,]  458
 [3,]    0
 [4,]    0
 [5,]    0
 [6,]    0
 [7,]    0
 [8,]   90
 [9,]   70
[10,]   65
[11,]   98
[12,]  106
[13,]   60
[14,]   80
[15,]  100
[16,]   85
[17,]   68
[18,]  225
[19,]  344
[20,]  253

Solutions to the Mixed Model Equations

• generalized inverse of LHS

gi_LHS = ginv(LHS)
gi_LHS

             [,1]        [,2]         [,3]         [,4]         [,5]         [,6]         [,7]        [,8]         [,9]        [,10]       [,11]       [,12]
 [1,]  0.44292085  0.13589744 -0.105239688 -0.041248606 -0.048606466 -0.048606466 -0.063991081 -0.11428465 -0.073348198 -0.073348198 -0.05052397 -0.09146042
 [2,]  0.13589744  0.34358974 -0.087179487 -0.051282051 -0.066666667 -0.066666667 -0.035897436 -0.06786325 -0.095555556 -0.095555556 -0.09025641 -0.06256410
 [3,] -0.10523969 -0.08717949  0.286733556 -0.011148272  0.020958751  0.020958751 -0.009810479  0.13568190  0.133273876  0.133273876  0.02559643  0.02800446
 [4,] -0.04124861 -0.05128205 -0.011148272  0.285395764  0.011148272  0.011148272  0.011148272  0.12879970  0.130137495  0.130137495  0.01426979  0.01293200
 [5,] -0.04860647 -0.06666667  0.020958751  0.011148272  0.286733556 -0.020958751  0.009810479  0.01816425  0.020572278  0.020572278  0.12824972  0.12584169
 [6,] -0.04860647 -0.06666667  0.020958751  0.011148272 -0.020958751  0.286733556  0.009810479  0.01816425  0.020572278  0.020572278  0.12824972  0.12584169
 [7,] -0.06399108 -0.03589744 -0.009810479  0.011148272  0.009810479  0.009810479  0.286733556  0.00688220  0.003136381  0.003136381  0.01132664  0.01507246
 [8,] -0.11428465 -0.06786325  0.135681903  0.128799703  0.018164251  0.018164251  0.006882200  0.26818927  0.128666419  0.128666419  0.02115793  0.02734745
 [9,] -0.07334820 -0.09555556  0.133273876  0.130137495  0.020572278  0.020572278  0.003136381  0.12866642  0.264960733  0.131627400  0.02645559  0.02349461
[10,] -0.07334820 -0.09555556  0.133273876  0.130137495  0.020572278  0.020572278  0.003136381  0.12866642  0.131627400  0.264960733  0.02645559  0.02349461
[11,] -0.05052397 -0.09025641  0.025596433  0.014269788  0.128249721  0.128249721  0.011326644  0.02115793  0.026455593  0.026455593  0.26163657  0.12300557
[12,] -0.09146042 -0.06256410  0.028004459  0.012931996  0.125841695  0.125841695  0.015072464  0.02734745  0.023494612  0.023494612  0.12300557  0.26019175
[13,] -0.05052397 -0.09025641  0.025596433  0.014269788  0.128249721  0.128249721  0.011326644  0.02115793  0.026455593  0.026455593  0.12830323  0.12300557
[14,] -0.05052397 -0.09025641  0.025596433  0.014269788  0.128249721  0.128249721  0.011326644  0.02115793  0.026455593  0.026455593  0.12830323  0.12300557
[15,] -0.11959123 -0.06427350  0.135994054  0.003270160  0.017852100  0.017852100  0.132723894  0.07563929  0.068263595  0.068263595  0.02047120  0.02784690
[16,] -0.07865478 -0.09196581  0.133586027  0.004607952  0.020260126  0.020260126  0.128978075  0.06944977  0.071224576  0.071224576  0.02576886  0.02399405
[17,] -0.11959123 -0.06427350  0.135994054  0.003270160  0.017852100  0.017852100  0.132723894  0.07563929  0.068263595  0.068263595  0.02047120  0.02784690
[18,] -0.06187291 -0.07692308 -0.016722408 -0.033444816  0.016722408  0.016722408  0.016722408 -0.03756968 -0.035562988 -0.035562988  0.02140468  0.01939799
[19,] -0.07290970 -0.10000000  0.031438127  0.016722408 -0.031438127 -0.031438127  0.014715719  0.02724638  0.030858417  0.030858417 -0.03839465 -0.04200669
[20,] -0.09598662 -0.05384615 -0.014715719  0.016722408  0.014715719  0.014715719 -0.031438127  0.01032330  0.004704571  0.004704571  0.01698997  0.02260870
            [,13]       [,14]       [,15]        [,16]       [,17]        [,18]       [,19]        [,20]
 [1,] -0.05052397 -0.05052397 -0.11959123 -0.078654775 -0.11959123 -0.061872910 -0.07290970 -0.095986622
 [2,] -0.09025641 -0.09025641 -0.06427350 -0.091965812 -0.06427350 -0.076923077 -0.10000000 -0.053846154
 [3,]  0.02559643  0.02559643  0.13599405  0.133586027  0.13599405 -0.016722408  0.03143813 -0.014715719
 [4,]  0.01426979  0.01426979  0.00327016  0.004607952  0.00327016 -0.033444816  0.01672241  0.016722408
 [5,]  0.12824972  0.12824972  0.01785210  0.020260126  0.01785210  0.016722408 -0.03143813  0.014715719
 [6,]  0.12824972  0.12824972  0.01785210  0.020260126  0.01785210  0.016722408 -0.03143813  0.014715719
 [7,]  0.01132664  0.01132664  0.13272389  0.128978075  0.13272389  0.016722408  0.01471572 -0.031438127
 [8,]  0.02115793  0.02115793  0.07563929  0.069449771  0.07563929 -0.037569677  0.02724638  0.010323300
 [9,]  0.02645559  0.02645559  0.06826359  0.071224576  0.06826359 -0.035562988  0.03085842  0.004704571
[10,]  0.02645559  0.02645559  0.06826359  0.071224576  0.06826359 -0.035562988  0.03085842  0.004704571
[11,]  0.12830323  0.12830323  0.02047120  0.025768859  0.02047120  0.021404682 -0.03839465  0.016989967
[12,]  0.12300557  0.12300557  0.02784690  0.023994054  0.02784690  0.019397993 -0.04200669  0.022608696
[13,]  0.26163657  0.12830323  0.02047120  0.025768859  0.02047120  0.021404682 -0.03839465  0.016989967
[14,]  0.12830323  0.26163657  0.02047120  0.025768859  0.02047120  0.021404682 -0.03839465  0.016989967
[15,]  0.02047120  0.02047120  0.26994773  0.129238697  0.13661439  0.004905240  0.02677815 -0.031683389
[16,]  0.02576886  0.02576886  0.12923870  0.264346835  0.12923870  0.006911929  0.03039019 -0.037302118
[17,]  0.02047120  0.02047120  0.13661439  0.129238697  0.26994773  0.004905240  0.02677815 -0.031683389
[18,]  0.02140468  0.02140468  0.00490524  0.006911929  0.00490524  0.180602007  0.02508361  0.025083612
[19,] -0.03839465 -0.03839465  0.02677815  0.030390190  0.02677815  0.025083612  0.18361204  0.022073579
[20,]  0.01698997  0.01698997 -0.03168339 -0.037302118 -0.03168339  0.025083612  0.02207358  0.183612040

solution

sol = gi_LHS %*% RHS
sol

            [,1]
 [1,] 91.4931401
 [2,] 75.7644444
 [3,] -1.4407729
 [4,] -1.1748792
 [5,]  1.4407729
 [6,]  1.4407729
 [7,] -0.2658937
 [8,] -1.0975588
 [9,] -1.6670660
[10,] -2.3337327
[11,]  3.9252560
[12,]  2.8947633
[13,] -1.1414106
[14,]  1.5252560
[15,]  0.4478712
[16,]  0.5450306
[17,] -3.8187955
[18,] -1.7623188
[19,]  2.1611594
[20,] -0.3988406

round(sol,3)

        [,1]
 [1,] 91.493
 [2,] 75.764
 [3,] -1.441
 [4,] -1.175
 [5,]  1.441
 [6,]  1.441
 [7,] -0.266
 [8,] -1.098
 [9,] -1.667
[10,] -2.334
[11,]  3.925
[12,]  2.895
[13,] -1.141
[14,]  1.525
[15,]  0.448
[16,]  0.545
[17,] -3.819
[18,] -1.762
[19,]  2.161
[20,] -0.399

sol_sex = sol[1:no_levels_sex]
sol_sex

[1] 91.49314 75.76444

sol_animal = sol[(no_levels_sex+1):(no_levels_sex+no_animals)]
sol_animal

 [1] -1.4407729 -1.1748792  1.4407729  1.4407729 -0.2658937 -1.0975588 -1.6670660 -2.3337327  3.9252560  2.8947633 -1.1414106  1.5252560  0.4478712  0.5450306
[15] -3.8187955

sol_ce = sol[(no_levels_sex+no_animals+1):(no_levels_sex+no_animals+no_animals_ce)]
sol_ce

[1] -1.7623188  2.1611594 -0.3988406

Reliability(r2), Accuracy(r), Standard Error of Prediction(SEP)

The diagonal elements of the generalized inverse of LHS for animal equation (D)

D = diag(gi_LHS[(no_levels_sex + 1) : (no_levels_sex + no_animals), (no_levels_sex + 1) : (no_levels_sex + no_animals)]) 
D

 [1] 0.2867336 0.2853958 0.2867336 0.2867336 0.2867336 0.2681893 0.2649607 0.2649607 0.2616366 0.2601918 0.2616366 0.2616366 0.2699477 0.2643468 0.2699477

The Reliability is calculated as:

r2 = 1 - D * alpha_1 
r2

 [1] 0.06811594 0.07246377 0.06811594 0.06811594 0.06811594 0.12838486 0.13887762 0.13887762 0.14968116 0.15437681 0.14968116 0.14968116 0.12266989 0.14087279
[15] 0.12266989

The accuracy is calculated as:

r = sqrt(r2) 

The standard error of prediction(SEP)

sep = sqrt(D * sigma_2_e) 
sep

 [1] 4.317138 4.307055 4.317138 4.317138 4.317138 4.175201 4.149994 4.149994 4.123879 4.112477 4.123879 4.123879 4.188866 4.145183 4.188866

Partitioning of Breeding values

Yield Deviations for the 10 animals with records

yd1 = ginv(t(Z) %*% Z) %*% t(Z) %*% (y - X %*% sol_sex - W %*% sol_ce)
yd1

             [,1]
 [1,]   0.0000000
 [2,]   0.0000000
 [3,]   0.0000000
 [4,]   0.0000000
 [5,]   0.0000000
 [6,]   0.2691787
 [7,]  -4.0021256
 [8,]  -9.0021256
 [9,]  20.0743961
[10,]  12.3457005
[11,] -17.9256039
[12,]   2.0743961
[13,]   8.9057005
[14,]   9.6343961
[15,] -23.0942995

The numerator of n2 (z’Z-number of records the animal has)

a2 = diag(t(Z) %*% Z)

Parents Average (PA) and Progeny Contribution (PC)

parents avearge

pa1 = rep(0, no_animals) 
# progeny contribution numerator 
pc0 = rep(0, no_animals) 
# numerator of n3, denominator of progeny contribution 
a3 = rep(0, no_animals)
# numerator of n1 
a1 = rep(0, no_animals) 

Lopping through the pedigree

for (i in 1 : no_animals) { 
  sire = ped[i, 2] 
  dam = ped[i, 3] 
  if (sire == 0 & dam == 0) {
    # both parents unknown 
    # PA 
    a1[i] = 1 * alpha_1 
  } else if (sire != 0 & dam == 0) {

    # PA 
    a1[i] = 4/3 * alpha_1 
    pa1[i] = sol_animal[sire]/2 
    # PC for sire 
    a3[sire] = a3[sire] + 0.5 * alpha_1 * (2/3) 
    pc0[sire] = pc0[sire] + 0.5 * alpha_1 * (2/3) * (2 * sol_animal[i]) 
  } else if (sire == 0 & dam != 0) { 
    # dam known 
    # PA 
    a1[i] = 4/3 * alpha_1 
    pa1[i] = sol_animal[dam]/2 
    # PC for dam 
    a3[dam] = a3[dam] + 0.5 * alpha_1 * (2/3) 
    pc0[dam] = pc0[dam] + 0.5 * alpha_1 * (2/3) * (2 * sol_animal[i]) 
  } else {
    # both parents known 
        # PA 
    a1[i] = 2 * alpha_1
    pa1[i] = (sol_animal[sire] + sol_animal[dam])/2
    # PC for sire 
    a3[sire] = a3[sire] + 0.5 * alpha_1 
    pc0[sire] = pc0[sire] + 0.5 * alpha_1 * (2 * sol_animal[i] - sol_animal[dam]) 
    # PC for dam 
    a3[dam] = a3[dam] + 0.5 * alpha_1
    pc0[dam] = pc0[dam] + 0.5 * alpha_1 * (2 * sol_animal[i] - sol_animal[sire]) 
  }
}

a1

 [1] 3.25 3.25 3.25 3.25 3.25 6.50 6.50 6.50 6.50 6.50 6.50 6.50 6.50 6.50 6.50

pa1

 [1]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000 -1.3078261 -1.3078261 -1.3078261  1.4407729  1.4407729  1.4407729  1.4407729 -0.8533333 -0.8533333
[15] -0.8533333

a3

 [1] 9.750 4.875 6.500 6.500 4.875 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

pc0

 [1] -18.730048  -9.545894  14.047536  14.047536  -2.160386   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
[15]   0.000000

The sum of the denominators of n1, n2, n3 is:

n_D = a1 +a2 + a3
n_D

 [1] 13.000  8.125  9.750  9.750  8.125  7.500  7.500  7.500  7.500  7.500  7.500  7.500  7.500  7.500  7.500

The PA fraction of the breeeding values

PA = pa1 * (a1/n_D)
PA

 [1]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000 -1.1334493 -1.1334493 -1.1334493  1.2486699  1.2486699  1.2486699  1.2486699 -0.7395556 -0.7395556
[15] -0.7395556

Yield Deviation (YD) fraction of the breeding values

YD = yd1 * (a2/n_D)
YD

            [,1]
 [1,]  0.0000000
 [2,]  0.0000000
 [3,]  0.0000000
 [4,]  0.0000000
 [5,]  0.0000000
 [6,]  0.0358905
 [7,] -0.5336167
 [8,] -1.2002834
 [9,]  2.6765862
[10,]  1.6460934
[11,] -2.3900805
[12,]  0.2765862
[13,]  1.1874267
[14,]  1.2845862
[15,] -3.0792399

Progeny Contribution (PC)

pc1 = pc0/a3
pc1[is.nan(pc1)==TRUE]=0
pc1

 [1] -1.9210306 -1.9581320  2.1611594  2.1611594 -0.4431562  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[15]  0.0000000

PC = pc1 * (a3/n_D)
PC

 [1] -1.4407729 -1.1748792  1.4407729  1.4407729 -0.2658937  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[15]  0.0000000

Progeny(Daughter) Yield Deviation(PYD, DYD)

# n2 of the offspring
n2prog = a2/(a1+a2)
n2prog

 [1] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1333333 0.1333333 0.1333333 0.1333333 0.1333333 0.1333333 0.1333333 0.1333333 0.1333333 0.1333333

The numerator of the DYD:
Summation of u of the progeny * n2 progeny * (2DYD-BV of the mate)
thus;

dyd_n = rep(0,no_animals)
# denominator of dyd : summation of u of progeny * n2 of progeny
dyd_d = rep(0, no_animals)

Looping through the pedigree gives:

 for (i in 1 : no_animals) { 
   sire = ped[i, 2] 
   dam = ped[i, 3] 
   if (sire == 0 & dam == 0) { 
     # both parents unknown 
    } else if (sire != 0 & dam == 0) { 
     # dyd_n
       dyd_n[sire] = dyd_n[sire] + n2prog[i] * 2 / 3 * 2 * yd1[i] 
       # dyd_d 
       dyd_d[sire] = dyd_d[sire] + n2prog[i] * 2 / 3
    } else if (sire == 0 & dam != 0) { 
         # dam known 
         # dyd_n 
         dyd_n[dam] = dyd_n[dam] + n2prog[i] * 2 / 3 * 2 * yd1[i] 
         # dyd_d 
         dyd_d[dam] = dyd_d[dam] + n2prog[i] * 2 / 3
    } else { 
        # both parents known 
         # dyd_n 
             dyd_n[sire] = dyd_n[sire] + n2prog[i] * (2 * yd1[i] - sol_animal[dam]) 
             dyd_n[dam] = dyd_n[dam] + n2prog[i] * (2 * yd1[i] - sol_animal[sire]) 
            # dyd_d 
             dyd_d[sire] = dyd_d[sire] + n2prog[i] 
             dyd_d[dam] = dyd_d[dam] + n2prog[i] 
    }
 }
dyd_n

 [1] -4.0341643 -2.8197101  3.6499581  3.6499581 -0.6381449  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[15]  0.0000000

dyd_d

 [1] 0.8000000 0.4000000 0.5333333 0.5333333 0.4000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000

DYD = dyd_n/dyd_d
DYD[is.nan(DYD) == TRUE] = 0
DYD

 [1] -5.042705 -7.049275  6.843671  6.843671 -1.595362  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000

Breeeding values and the partitions

BV_result  = data.frame(animal = ped[,1], animal_bv = sol_animal, r2 = r2, r
           = r, SEP = sep, PA = PA, YD = YD, PC = PC, sum_of_fr = PA + YD+ PC, 
           DYD = DYD)
round(BV_result,3)

BV_result