Pedigree BLUP

# RCBD
# 4 gens
# 3 blocks
# 12 observations
data <- read.csv("../data/example_1.csv") |>
  mutate(gen = as.factor(gen), block = as.factor(block))
head(data)

  block gen yield
1     1  g1   7.4
2     2  g1   6.5
3     3  g1   5.6
4     1  g2   9.8
5     2  g2   6.8
6     3  g2   6.2

n <- 12
n_blks <- 3
n_gens <- 4
g_lvls <- levels(data$gen)

# Pedigree
ped <- data.frame(
  id = g_lvls,
  p1 = c(NA, NA, "g2", "g3"),
  p2 = c(NA, NA, "g1", "g1")
)

Additive relationship matrix

A <- A_mat(ped, ploidy = 2)
A

     g1   g2   g3   g4
g1 1.00 0.00 0.50 0.75
g2 0.00 1.00 0.50 0.25
g3 0.50 0.50 1.00 0.75
g4 0.75 0.25 0.75 1.25

Variance components

Online License checked out Tue Feb 17 21:45:48 2026

asreml.options(Cfixed = TRUE, maxit = 50, trace = 0)
mme <- asreml(fixed = yield ~ 1, random = ~ vm(gen, A) + block, data = data)
mme$coefficients$random

                  effect
vm(gen, A)_g1 -0.4594011
vm(gen, A)_g2  0.4594011
vm(gen, A)_g3 -0.3201192
vm(gen, A)_g4  1.0022245
block_1        1.1439320
block_2       -0.3660582
block_3       -0.7778737
attr(,"terms")
                tname n
vm(gen, A) vm(gen, A) 4
block           block 3

var_comps <- summary(mme)$varcomp
var_g <- var_comps[2, 1]
var_b <- var_comps[1, 1]
var_e <- var_comps[3, 1]

Reconstructing MME

X <- model.matrix(yield ~ 1, data)
Zb <- model.matrix(yield ~ -1 + block, data)
Zg <- model.matrix(yield ~ -1 + gen, data)
Z <- cbind(Zb, Zg)
y <- matrix(data[, "yield"])

Gb <- diag(x = var_b, nrow = n_blks)
Gg <- A * var_g
G <- bdiag(Gb, Gg)
R <- diag(x = var_e, nrow = n)
V <- Z %*% G %*% t(Z) + R

# Mixed Model Equations
C11 <- t(X) %*% chol2inv(chol(R)) %*% X
C12 <- t(X) %*% chol2inv(chol(R)) %*% Z
C21 <- t(Z) %*% chol2inv(chol(R)) %*% X
C22 <- t(Z) %*% chol2inv(chol(R)) %*% Z + solve(G)

# Coefficient matrix (LHS)
C <- as.matrix(
  rbind(
    cbind(C11, C12),
    cbind(C21, C22)
  )
)

# RHS
rhs <- rbind(
  t(X) %*% solve(R) %*% y,
  t(Z) %*% solve(R) %*% y
)

# Solution
C_inv <- ginv(C)
rownames(C_inv) <- colnames(C_inv) <- rownames(C)
print(round(C_inv, 3))

            (Intercept) block1 block2 block3  geng1  geng2  geng3  geng4
(Intercept)       1.321 -0.373 -0.373 -0.373 -0.885 -0.883 -0.941 -0.943
block1           -0.373  0.437  0.342  0.342  0.000  0.000  0.000  0.000
block2           -0.373  0.342  0.437  0.342  0.000  0.000  0.000  0.000
block3           -0.373  0.342  0.342  0.437  0.000  0.000  0.000  0.000
geng1            -0.885  0.000  0.000  0.000  0.947  0.821  0.883  0.889
geng2            -0.883  0.000  0.000  0.000  0.821  0.947  0.885  0.879
geng3            -0.941  0.000  0.000  0.000  0.883  0.885  1.051  0.945
geng4            -0.943  0.000  0.000  0.000  0.889  0.879  0.945  1.059

ans <- C_inv %*% rhs
ans

                  [,1]
(Intercept)  7.0794023
block1       1.1440635
block2      -0.3661003
block3      -0.7779632
geng1       -0.4596044
geng2        0.4596044
geng3       -0.3203845
geng4        1.0027751

# PEV = C22_g
C22_g <- C_inv[5:8, 5:8]
C22_g

          geng1     geng2     geng3     geng4
geng1 0.9472216 0.8208755 0.8832245 0.8892003
geng2 0.8208755 0.9472216 0.8848727 0.8788968
geng3 0.8832245 0.8848727 1.0513126 0.9445423
geng4 0.8892003 0.8788968 0.9445423 1.0589367

Reliability

PEV_gi <- diag(C22_g)
var_gi <- diag(Gg)
r2_g <- 1 - PEV_gi / var_gi
r2_g

    geng1     geng2     geng3     geng4 
0.4642706 0.4642706 0.4053989 0.5208695 

# Linear combination w = Lb + Mu
var_uhat <- G - C_inv[-1, -1]
var_beta <- C_inv[1, 1]
cov_ubeta <- -matrix(C_inv[-1, 1], nrow = 7, ncol = 1)

L <- cbind(matrix(1, nrow = n_gens, ncol = 1)) # Intercept
M <- cbind(
  matrix(0, nrow = 4, ncol = 3),               # No block
  diag(n_gens)                                 # Genotype 
)
var_what <- L %*% var_beta %*% t(L) + M %*% var_uhat %*% t(M)
var_w <- M %*% G %*% t(M)
cov_wwhat <- M %*% cov_ubeta %*% t(L) + M %*% var_uhat %*% t(M)
r2_w <- diag(cov_wwhat)^2 / (diag(var_w) * diag(var_what))
r2_w

[1] 0.7685843 0.7666359 0.7628126 0.8026305

Predictions

# g
pp <- predict(mme, classify = "vm(gen, A)", only = "vm(gen, A)", sed = TRUE, vcov = TRUE)
as.data.frame(pp$pvals) # Solution

  vm.gen..A. predicted.value std.error    status
1         g1      -0.4596291 0.9729385 Estimable
2         g2       0.4596291 0.9729385 Estimable
3         g3      -0.3204166 1.0249924 Estimable
4         g4       1.0028419 1.0287003 Estimable

pp$vcov # C22_g

4 x 4 Matrix of class "dspMatrix"
          [,1]      [,2]      [,3]      [,4]
[1,] 0.9466094 0.8203812 0.8826721 0.8886408
[2,] 0.8203812 0.9466094 0.8843185 0.8783498
[3,] 0.8826721 0.8843185 1.0506093 0.9439335
[4,] 0.8886408 0.8783498 0.9439335 1.0582244

pp$sed^2 # vd_BLUP_mat

4 x 4 Matrix of class "dspMatrix"
          [,1]      [,2]      [,3]      [,4]
[1,]        NA 0.2524565 0.2318745 0.2275522
[2,] 0.2524565        NA 0.2285818 0.2481342
[3,] 0.2318745 0.2285818        NA 0.2209668
[4,] 0.2275522 0.2481342 0.2209668        NA

# mu + g
lc <- predict(mme, classify = "gen", sed = TRUE, vcov = TRUE)
as.data.frame(lc$pvals)

  gen predicted.value std.error    status
1  g1        6.619765 0.7052467 Estimable
2  g2        7.539023 0.7083045 Estimable
3  g3        6.758977 0.6998271 Estimable
4  g4        8.082236 0.7025421 Estimable

lc$vcov

4 x 4 Matrix of class "dspMatrix"
          [,1]      [,2]      [,3]      [,4]
[1,] 0.4973729 0.3733059 0.3776282 0.3816931
[2,] 0.3733059 0.5016952 0.3814357 0.3735632
[3,] 0.3776282 0.3814357 0.4897579 0.3811783
[4,] 0.3816931 0.3735632 0.3811783 0.4935654

lc$sed^2

4 x 4 Matrix of class "dspMatrix"
          [,1]      [,2]      [,3]      [,4]
[1,]        NA 0.2524565 0.2318745 0.2275522
[2,] 0.2524565        NA 0.2285818 0.2481342
[3,] 0.2318745 0.2285818        NA 0.2209668
[4,] 0.2275522 0.2481342 0.2209668        NA