Resolução das equações de modelos mistos
com grupos de pais
desconhecidos (UPG)
1) Resolução direta (R)
2) Método iterativo (BLUPF90+)
Autor: Eula Carrara (eulacarrara@gmail.com)
Criado em: 22-Feb-2025
Modificado em: 15-Mar-2025
Teoria: notas de aula Dr. Daniela Lourenco, UGA,
2025
Exemplo numérico retirado de tutorial_blupf90.pdf
Pacotes R exigidos: igraph, MASS, optiSel,
tidyverse
install.packages(c("igraph", "MASS", "optiSel", "tidyverse"))Considerando um individuo a, filho de s e d,
\[\begin{array}{c} s_i \quad \quad d_i \\ \backslash \quad \, / \\ \, a_i \end{array}\]o seu valor genético esperado será:
se ambos pais conhecidos: \(a_i = \frac{a_{s_i} + a_{d_i}}{2}\)
se apenas pai conhecido: \(a_i = \frac{a_{s_i} + 0}{2}\)
se apenas mãe conhecida: \(a_i = \frac{0 + a_{d_i}}{2}\)
se ambos pais desconhecidos: \(a_i = \frac{0 + 0}{2}\)
Ou seja, o valor genético de um pai/mãe desconhecido é igual a zero.
O que não condiz com a realidade.
Uma maneira de preencher essas lacunas na genealogia é construir grupos
de pais desconhecidos (UPG; ou pais fantasmas, ou grupo genético):
apenas pai conhecido: \(a_i = \frac{a_{s_i} + \boldsymbol{\mathit{UPG_i}}}{2}\)
apenas mãe conhecida: \(a_i = \frac{\boldsymbol{\mathit{UPG_i}} + a_{d_i}}{2}\)
ambos pais desconhecidos: \(a_i = \frac{\boldsymbol{\mathit{UPG_i}} + \boldsymbol{\mathit{UPG_i}}}{2}\)
Como formar UPGs?
Utilizando informacoes de raça, ano de nascimento, sexo, país,
etc…
UPGs não são animais! São tratados, na maioria das vezes, como um efeito
fixo no modelo.
Modelo animal misto
O modelo animal misto pode ser descrito como:
\(y = Xb + Za + e\)
em que:
- y é o vetor de fenótipos,
- b é o vetor de efeitos fixos,
- a é o vetor de efeitos genéticos aditivos diretos (valor
genético),
- X e Z relacionam b e
a à y,
- e é o vetor de resíduos.
Modelo animal misto, mas com UPG:
\(y = Xb + \mathbf{\mathit{ZQg}} + Za + e\)
em que:
- Q é a matriz que relaciona os animais aos UPGs
- g é o efeito dos UPG.
Então precisamos construir Q. Como?
A matriz Q é construída com as frações de
contribuição de cada UPG nos valores genéticos esperados dos
indivíduos relacionados aos UPG.
Exemplo de construção da matriz Q:
Genealogia:
## ID Sire Dam
## 1 a1 g2 a4
## 2 a2 g3 a1
## 3 a3 g3 a1
## 4 a4 g1 g1
## 5 a5 a3 a6
## 6 a6 g3 g3
Primeiro, precisamos colocar todos os valores genéticos esperados em função dos UPG:
Definindo os valores genéticos esperados em função dos UPG \[\begin{aligned} E(a_4) &= \frac{1}{2} \cdot g_1 + \frac{1}{2} \cdot g_1 = g_1 \\ \\[0.1cm] E(a_1) &= \frac{1}{2} \cdot g_2 + \frac{1}{2} \cdot E(a_4) = \frac{1}{2} g_1 + \frac{1}{2} g_2 \\ \\[0.1cm] E(a_3) &= \frac{1}{2} \cdot g_3 + \frac{1}{2} \cdot E(a_1) = \frac{1}{4} g_1 + \frac{1}{4} g_2 + \frac{1}{2} g_3 \\ \\[0.1cm] E(a_6) &= \frac{1}{2} \cdot g_3 + \frac{1}{2} \cdot g_3 = g_3 \\ \\[0.1cm] E(a_5) &= \frac{1}{2} \cdot E(a_3) + \frac{1}{2} \cdot E(a_6) = \frac{1}{8} g_1 + \frac{1}{8} g_2 + \frac{3}{4} g_3 \\ \\[0.1cm] E(a_2) &= \frac{1}{2} \cdot g_3 + \frac{1}{2} \cdot E(a_1) = \frac{1}{4} g_1 + \frac{1}{4} g_2 + \frac{1}{2} g_3 \end{aligned}\] Reorganizando e colocando zero nos grupos faltantes (apenas para melhor visualização) \[\begin{aligned} E(a_1) &= \frac{1}{2} \cdot g_1 + \frac{1}{2} \cdot g_2 + 0 \cdot g_3 \\ \\[0.1cm] E(a_2) &= \frac{1}{4} \cdot g_1 + \frac{1}{4} \cdot g_2 + \frac{1}{2} \cdot g_3 \\ \\[0.1cm] E(a_3) &= \frac{1}{4} \cdot g_1 + \frac{1}{4} \cdot g_2 + \frac{1}{2} \cdot g_3 \\ \\[0.1cm] E(a_4) &= 1 \cdot g_1 + 0 \cdot g_2 + 0 \cdot g_3 \\ \\[0.1cm] E(a_5) &= \frac{1}{8} \cdot g_1 + \frac{1}{8} \cdot g_2 + \frac{3}{4} \cdot g_3 \\ \\[0.1cm] E(a_6) &= 0 \cdot g_1 + 0 \cdot g_2 + 1 \cdot g_3 \end{aligned}\]Matriz Q com as frações de cada UPG \[ \begin{bmatrix} 1/2 & 1/2 & 0 \\ 1/4 & 1/4 & 1/2 \\ 1/4 & 1/4 & 1/2 \\ 1 & 0 & 0 \\ 1/8 & 1/8 & 3/4 \\ 0 & 0 & 1 \end{bmatrix} \]
Essa é a matriz Q que será incorporada nas equações de
modelos mistos.
Exemplo da resolução das equações de modelos mistos considerando UPG
Resolução via programa R, passo a passo
Dados
data1 <- data.frame(
ID = c('ID006', 'ID009', 'ID012', 'ID007', 'ID010', 'ID013', 'ID008', 'ID011', 'ID014', 'ID015'),
A = c('A', 'A', 'A', 'B', 'B', 'B', 'C', 'C', 'C', 'C'),
S = c(1, 2, 1, 2, 1, 2, 1, 2, 1, 2),
cov = c(1.0, 1.0, 2.0, 2.0, 1.0, 2.0, 2.0, 1.0, 1.0, 2.0),
obs = c(3.0, 2.0, 4.0, 6.0, 3.0, 6.0, 6.0, 6.0, 8.0, 4.0)
)
print(data1)## ID A S cov obs
## 1 ID006 A 1 1 3
## 2 ID009 A 2 1 2
## 3 ID012 A 1 2 4
## 4 ID007 B 2 2 6
## 5 ID010 B 1 1 3
## 6 ID013 B 2 2 6
## 7 ID008 C 1 2 6
## 8 ID011 C 2 1 6
## 9 ID014 C 1 1 8
## 10 ID015 C 2 2 4Genealogia sem UPG
pedigree <- data.frame(
ID = c('ID001', 'ID002', 'ID003', 'ID004', 'ID005', 'ID006', 'ID007', 'ID008', 'ID009', 'ID010', 'ID011', 'ID012', 'ID013', 'ID014', 'ID015'),
Sire = c(NA, NA, NA, NA, NA, NA, 'ID002', 'ID001', 'ID002', 'ID007', 'ID007', 'ID011', 'ID011', 'ID009', 'ID011'),
Dam = c(NA, NA, NA, NA, NA, NA, 'ID005', 'ID004', 'ID003', 'ID006', 'ID004', 'ID008', 'ID010', 'ID013', 'ID010')
)
print(pedigree)## ID Sire Dam
## 1 ID001 <NA> <NA>
## 2 ID002 <NA> <NA>
## 3 ID003 <NA> <NA>
## 4 ID004 <NA> <NA>
## 5 ID005 <NA> <NA>
## 6 ID006 <NA> <NA>
## 7 ID007 ID002 ID005
## 8 ID008 ID001 ID004
## 9 ID009 ID002 ID003
## 10 ID010 ID007 ID006
## 11 ID011 ID007 ID004
## 12 ID012 ID011 ID008
## 13 ID013 ID011 ID010
## 14 ID014 ID009 ID013
## 15 ID015 ID011 ID010Genealogia com UPG
pedigree_upg <- data.frame(
ID = c('ID001', 'ID002', 'ID003', 'ID004', 'ID005', 'ID006', 'ID007', 'ID008', 'ID009', 'ID010', 'ID011', 'ID012', 'ID013', 'ID014', 'ID015'),
Sire = c('g1', 'g2', 'g1', 'g2', 'g2', 'g1', 'ID002', 'ID001', 'ID002', 'ID007', 'ID007', 'ID011', 'ID011', 'ID009', 'ID011'),
Dam = c('g4', 'g3', 'g3', 'g3', 'g4', 'g3', 'ID005', 'ID004', 'ID003', 'ID006', 'ID004', 'ID008', 'ID010', 'ID013', 'ID010')
)
print(pedigree_upg)## ID Sire Dam
## 1 ID001 g1 g4
## 2 ID002 g2 g3
## 3 ID003 g1 g3
## 4 ID004 g2 g3
## 5 ID005 g2 g4
## 6 ID006 g1 g3
## 7 ID007 ID002 ID005
## 8 ID008 ID001 ID004
## 9 ID009 ID002 ID003
## 10 ID010 ID007 ID006
## 11 ID011 ID007 ID004
## 12 ID012 ID011 ID008
## 13 ID013 ID011 ID010
## 14 ID014 ID009 ID013
## 15 ID015 ID011 ID010# Ordenar o pedigree (sem UPG) pelo ID do animal
pedigree <- pedigree[order(pedigree$ID), ]
# Ordenar os dados pelo ID do animal
data1 <- data1[order(data1$ID), ]Modelo
\(y = Xb + ZQg + Za + e\)
Equações de modelos mistos
\[ \begin{bmatrix} X'X & 0 & X'Z \\ 0 & Q'A^{-1}Q\alpha & -Q'A^{-1}\alpha \\ Z'X & -A^{-1}Q\alpha & Z'Z + A^{-1}\alpha \end{bmatrix} * \begin{bmatrix} b\\ g\\ Qg + a \end{bmatrix} = \begin{bmatrix} X'y \\ 0 \\ Z'y \end{bmatrix} \]
Elementos das equações de modelos mistos
Vetor y
y <- as.matrix(data1$obs)
print(y)## [,1]
## [1,] 3
## [2,] 6
## [3,] 6
## [4,] 2
## [5,] 3
## [6,] 6
## [7,] 4
## [8,] 6
## [9,] 8
## [10,] 4Matriz X
library(dplyr)
X <- data1 %>%
mutate(A_A = ifelse(A == "A", 1, 0),
A_B = ifelse(A == "B", 1, 0),
A_C = ifelse(A == "C", 1, 0),
S_1 = ifelse(S == 1, 1, 0),
S_2 = ifelse(S == 2, 1, 0)) %>%
select(A_A, A_B, A_C, S_1, S_2, cov)
X <- as.matrix(X)
print(X)## A_A A_B A_C S_1 S_2 cov
## 1 1 0 0 1 0 1
## 4 0 1 0 0 1 2
## 7 0 0 1 1 0 2
## 2 1 0 0 0 1 1
## 5 0 1 0 1 0 1
## 8 0 0 1 0 1 1
## 3 1 0 0 1 0 2
## 6 0 1 0 0 1 2
## 9 0 0 1 1 0 1
## 10 0 0 1 0 1 2Matriz Z (expandindo a Z para todos os animais para somar com a matriz dos numeradores dos relacionamentos \(A^{-1}\))
animals <- pedigree$ID # Lista de animais na genealogia
N <- nrow(data1) # Número de observações
Np <- length(animals) # Número de animais na genealogia
Z <- matrix(0, nrow = N, ncol = Np)
for (i in 1:N) {
animal_index <- match(data1$ID[i], animals)
Z[i, animal_index] <- 1
}
print(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 1Construindo as matrizes A e A inversa (note que a A inversa é construída utilizando o pedigree sem UPG)
library(optiSel)
# Matriz A
A <- makeA(pedigree)# Mostrando apenas as 5 primeiras linhas/colunas
print(A[1:5, 1:5], digits = 4)## ID001 ID002 ID003 ID004 ID005
## ID001 1 0 0 0 0
## ID002 0 1 0 0 0
## ID003 0 0 1 0 0
## ID004 0 0 0 1 0
## ID005 0 0 0 0 1# Matriz A inversa
Ainv <- solve(A)# Mostrando apenas as 5 primeiras linhas/colunas
print(Ainv[1:5, 1:5], digits = 4)## ID001 ID002 ID003 ID004 ID005
## ID001 1.5 0.0 0.0 0.5 0.0
## ID002 0.0 2.0 0.5 0.0 0.5
## ID003 0.0 0.5 1.5 0.0 0.0
## ID004 0.5 0.0 0.0 2.0 0.0
## ID005 0.0 0.5 0.0 0.0 1.5Matriz Q
- Função para computar a matriz Q:
n_UPG = 4
compute_expectation <- function(pedigree, N) {
expectations <- list()
base_groups <- unique(c(pedigree$Sire, pedigree$Dam))
base_groups <- base_groups[grepl("^g[0-9]+$", base_groups)]
for (group in base_groups) {
expectations[[group]] <- setNames(as.list(rep(0, N)), paste0("g", 1:N))
expectations[[group]][[group]] <- 1
}
get_expectation <- function(id) {
if (!is.null(expectations[[id]])) {
return(expectations[[id]])
}
row <- pedigree[pedigree$ID == id, ]
if (nrow(row) == 0) {
return(setNames(as.list(rep(0, N)), paste0("g", 1:N)))
}
sire_exp <- get_expectation(row$Sire)
dam_exp <- get_expectation(row$Dam)
expectations[[id]] <- setNames(lapply(1:N, function(i) {
(1/2) * sire_exp[[paste0("g", i)]] + (1/2) * dam_exp[[paste0("g", i)]]
}), paste0("g", 1:N))
return(expectations[[id]])
}
animals <- pedigree$ID[grepl("^ID[0-9]+$", pedigree$ID)]
for (animal in animals) {
expectations[[animal]] <- get_expectation(animal)
}
return(expectations)
}Aplicando a função na genealogia com UPG (formato: ID, Pai, Mãe, Pai/Mãe com UPG)
expectations <- compute_expectation(pedigree_upg, n_UPG)
coefficients <- data.frame(
ID = names(expectations),
g1 = sapply(expectations, function(x) x$g1),
g2 = sapply(expectations, function(x) x$g2),
g3 = sapply(expectations, function(x) x$g3),
g4 = sapply(expectations, function(x) x$g4))Matriz Q final:
Q <- as.matrix(coefficients[-(1:4), -1])
print(Q)## g1 g2 g3 g4
## ID001 0.5000 0.0000 0.0000 0.5000
## ID002 0.0000 0.5000 0.5000 0.0000
## ID003 0.5000 0.0000 0.5000 0.0000
## ID004 0.0000 0.5000 0.5000 0.0000
## ID005 0.0000 0.5000 0.0000 0.5000
## ID006 0.5000 0.0000 0.5000 0.0000
## ID007 0.0000 0.5000 0.2500 0.2500
## ID008 0.2500 0.2500 0.2500 0.2500
## ID009 0.2500 0.2500 0.5000 0.0000
## ID010 0.2500 0.2500 0.3750 0.1250
## ID011 0.0000 0.5000 0.3750 0.1250
## ID012 0.1250 0.3750 0.3125 0.1875
## ID013 0.1250 0.3750 0.3750 0.1250
## ID014 0.1875 0.3125 0.4375 0.0625
## ID015 0.1250 0.3750 0.3750 0.1250Vamos considerar que os componentes de variância são conhecidos:
Variância aditiva direta: \(0.5\)
vara <- 0.5Variância residual: \(2.0\)
vare <- 2.0Portanto, alpha é igual a:
alpha <- vare / vara
print(alpha)## [1] 4Temos os elementos individuais das equações de modelos mistos, ou seja, \(X\), \(Z\), \(Q\), \(A^{-1}\), \(\alpha\) e \(y\).
Agora, precisamos calcular cada “bloco”:
\[ \begin{bmatrix} X'X & 0 & X'Z \\ 0 & Q'A^{-1}Q\alpha & -Q'A^{-1}\alpha \\ Z'X & -A^{-1}Q\alpha & Z'Z + A^{-1}\alpha \end{bmatrix} * \begin{bmatrix} b \\ g \\ Qg + a \end{bmatrix} = \begin{bmatrix} X'y \\ 0 \\ Z'y \end{bmatrix} \]
\[ \begin{bmatrix} \text{block11} & \text{block12} & \text{block13} \\ \text{block21} & \text{block22} & \text{block23} \\ \text{block31} & \text{block32} & \text{block33} \end{bmatrix} * \begin{bmatrix} b \\ g \\ Qg + a \end{bmatrix} = \begin{bmatrix} \text{block1} \\ \text{block2} \\ \text{block3} \end{bmatrix} \]
Lembrando que LHS * sol = RHS
LHS - lado esquerdo das equações de modelos mistos
block_11 <- t(X) %*% X
block_12 <- matrix(0, nrow = nrow(block_11), ncol = ncol(Q))
block_13 <- t(X) %*% Z
block_21 <- matrix(0, nrow = ncol(Q), ncol = nrow(block_11))
block_22 <- t(Q) %*% Ainv %*% Q * alpha
block_23 <- -t(Q) %*% Ainv * alpha
block_31 <- t(Z) %*% X
block_32 <- -Ainv %*% Q * alpha
block_33 <- t(Z) %*% Z + Ainv * alpha
LHS <- rbind(
cbind(block_11, block_12, block_13),
cbind(block_21, block_22, block_23),
cbind(block_31, block_32, block_33))
dimnames(LHS) <- NULL# Mostrando apenas as 5 primeiras linhas/colunas
print(LHS[1:5, 1:5], digits = 4)## [,1] [,2] [,3] [,4] [,5]
## [1,] 3 0 0 2 1
## [2,] 0 3 0 1 2
## [3,] 0 0 4 2 2
## [4,] 2 1 2 5 0
## [5,] 1 2 2 0 5RHS - lado direito das equações de modelos mistos
block_1 <- t(X) %*% y
block_2 <- matrix(0, nrow = ncol(Q), ncol = ncol(y))
block_3 <- t(Z) %*% y
RHS <- rbind(block_1, block_2, block_3)
dimnames(RHS) <- NULL# Mostrando apenas as 5 primeiras linhas/colunas
print(RHS[1:5, 1], digits = 4)## [1] 9 15 24 24 24Nós precisamos inverter o LHS, pois:
\(\text{sol} = \text{LHS}^{-1} \times
\text{RHS}\)
solve(LHS)! system is computationally singular: reciprocal condition number = 1.58711e-18
Percebam que o LHS não é uma matriz positiva
definida (system is computationally singular),
então não possui inversa.
Precisamos calcular a inversa generalizada…
Vamos usar a função ginv() do pacote MASS,
que calcula a inversa generalizada de Moore-Penrose (uma das mais
conhecidas)
library(MASS)
LHSinv <- ginv(LHS)# Mostrando apenas as 5 primeiras linhas/colunas
print(LHSinv[1:5, 1:5], digits = 4)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.4778 0.4631 0.6591 0.9670 1.633
## [2,] 0.4631 0.7888 0.4141 0.6595 1.007
## [3,] 0.6591 0.4141 0.8001 0.6504 1.223
## [4,] 0.9670 0.6595 0.6504 1.0605 1.216
## [5,] 1.6330 1.0065 1.2230 1.2164 2.646Multiplicamos a inversa do LHS pelo RHS e teremos nosso vetor de soluções.
Soluções finais (R):
solutions_R <- LHSinv %*% RHS
print(solutions_R)## [,1]
## [1,] -0.53552906
## [2,] 1.38930712
## [3,] 2.13580021
## [4,] 2.38873040
## [5,] 0.60084788
## [6,] 1.02099858
## [7,] -4.60399744
## [8,] 2.90828845
## [9,] 3.30986193
## [10,] -5.00557093
## [11,] -4.79790276
## [12,] 3.12605649
## [13,] -0.62320503
## [14,] 3.09897532
## [15,] -1.05552267
## [16,] -0.67781191
## [17,] 1.02838548
## [18,] -0.84258229
## [19,] 1.27528846
## [20,] -0.05636038
## [21,] 2.04669910
## [22,] 0.54680764
## [23,] 1.16463449
## [24,] 1.34945546
## [25,] 0.79807885A ordem das soluções é: \[
\begin{bmatrix}
b \\
g \\
Qg + a
\end{bmatrix}
\] em que:
- b é o vetor de soluções para os efeitos fixos
- g é o vetor de soluções para os efeitos dos UPG (também
fixos aqui)
- Qg + a é o vetor de valores genéticos.
Reparem que o valor genético do animal será o
efeito do UPG (Qg) + valor genético (a).
ou seja:
Vetor de efeitos fixos b
beta_hat <- solutions_R[1:ncol(X)]
b <- as.matrix(beta_hat)
print(b)## [,1]
## [1,] -0.5355291
## [2,] 1.3893071
## [3,] 2.1358002
## [4,] 2.3887304
## [5,] 0.6008479
## [6,] 1.0209986Vetor de efeitos dos UPG (g)
upg_hat <- solutions_R[(ncol(X) + 1):(ncol(X) + ncol(Q))]
g <- as.matrix(upg_hat)
print(g)## [,1]
## [1,] -4.603997
## [2,] 2.908288
## [3,] 3.309862
## [4,] -5.005571Vetor de efeitos genéticos aditivos diretos - valor genético (Qg+a)
animal_hat <- solutions_R[(ncol(X) + ncol(Q) + 1):(ncol(X) + ncol(Q) + ncol(Z))]
Qg_a <- as.matrix(animal_hat)
print(Qg_a)## [,1]
## [1,] -4.79790276
## [2,] 3.12605649
## [3,] -0.62320503
## [4,] 3.09897532
## [5,] -1.05552267
## [6,] -0.67781191
## [7,] 1.02838548
## [8,] -0.84258229
## [9,] 1.27528846
## [10,] -0.05636038
## [11,] 2.04669910
## [12,] 0.54680764
## [13,] 1.16463449
## [14,] 1.34945546
## [15,] 0.79807885E essas foram as nossas soluções finais!
y predito
Agora, vamos calcular o vetor y (y predito) a partir das
soluções que encontramos.
Pra isso, vamos isolar o y da primeira e da terceira linha
das equações de modelos mistos, pois a segunda linha do RHS é igual a
zero.
Da primeira linha:
\[
X'y = X'Xb + X'Z(Qg+a)
\]
rhs1 <- t(X) %*% X %*% b + t(X) %*% Z %*% Qg_aDa terceira linha:
\[
Z'y = Z'Xb - A^{-1}Q\alpha g + (Z'Z + A^{-1} \alpha)(Qg + a)
\]
rhs3 <- t(Z) %*% X %*% b - Ainv %*% Q %*% (alpha * g) + (t(Z) %*% Z + Ainv * alpha) %*% Qg_aVamos criar um sistema com essas duas equações, para resolver
y
\[
\left\{
\begin{aligned}
X'y &= X'bX + X'Z(Qg+a) \\
Z'y &= Z'Xb - A^{-1}Q\alpha g + (Z'Z + A^{-1} \alpha)(Qg
+ a)
\end{aligned}
\right.
\]
Qualquer que seja o modelo animal, a parte que multiplica
y será sempre uma matriz mxn. Vamos
chamá-la de M.
Da mesma forma, toda a parte depois do sinal de igual, será um
vetor nx1. Vamos chamá-lo de c.
Então, nos temos que \(My=c\).
Matriz M
M <- rbind(t(X), t(Z))Vetor c
c <- rbind(rhs1, rhs3)Resolvendo y…
Se tentarmos inverter a matriz M, veremos que ela não é
positiva definida, então precisamos de uma inversa generalizada:
ginv().
y_hat <- ginv(M) %*% c
final <- data.frame(y_predito = y_hat, y_original = y)
print(final)## y_predito y_original
## 1 3 3
## 2 6 6
## 3 6 6
## 4 2 2
## 5 3 3
## 6 6 6
## 7 4 4
## 8 6 6
## 9 8 8
## 10 4 4O y predito é idêntico ao y original.
Isso é fácil para 15 animais. Mas e se tivermos 1 milhão deles?
Precisamos de um programa mais eficiente.
Resolução via programa BLUPF90+
Vamos resolver esse mesmo exemplo utilizando os programas da família
BLUPF90
(Misztal I., Tsuruta S., Lourenco D.A.L., Aguilar I., Legarra A.,
and Vitezica Z. 2014. Manual for BLUPF90 family of programs.).
Data
ID006 A 1 1.0 3.0
ID009 A 2 1.0 2.0
ID012 A 1 2.0 4.0
ID007 B 2 2.0 6.0
ID010 B 1 1.0 3.0
ID013 B 2 2.0 6.0
ID008 C 1 2.0 6.0
ID011 C 2 1.0 6.0
ID014 C 1 1.0 8.0
ID015 C 2 2.0 4.0
Pedigree com UPG
ID001 0 0 -1 -4
ID002 0 0 -2 -3
ID003 0 0 -1 -3
ID004 0 0 -2 -3
ID005 0 0 -2 -4
ID006 0 0 -1 -3
ID007 ID002 ID005 ID002 ID005
ID008 ID001 ID004 ID001 ID004
ID009 ID002 ID003 ID002 ID003
ID010 ID007 ID006 ID007 ID006
ID011 ID007 ID004 ID007 ID004
ID012 ID011 ID008 ID011 ID008
ID013 ID011 ID010 ID011 ID010
ID014 ID009 ID013 ID009 ID013
ID015 ID011 ID010 ID011 ID010
Aqui nós codificamos os UPG com números negativos ao invés de letras. É uma exigência do programa BLUPF90+!
Arquivo de parâmetros
DATAFILE
data1.txt
TRAITS
5
FIELDS_PASSED TO OUTPUT
WEIGHT(S)
RESIDUAL_VARIANCE
2.0
EFFECT
2 cross alpha
EFFECT
3 cross alpha
EFFECT
4 cov
EFFECT
1 cross alpha
RANDOM
animal
FILE
ped2.txt
FILE_POS
1 4 5 0 0 # id, pai, mae - com codigos upg
UPG_TYPE
in_ped
INBREEDING
no-inbreeding # desconsidera a endogamia
(CO)VARIANCES
0.5
Rodamos o programa RENUMF90
e então o programa BLUPF90+.
As soluções (arquivo solutions) serão:
sol_blup <- data.frame(
trait = rep(1, 25),
effect = c(rep(1, 3), rep(2, 2), 3, rep(4, 19)),
level = c(1,2,3, 1,2, 1, 1:19),
solution = c(
-0.09929111, 1.82312339, 2.56875029, 2.43588849, 0.65341111,
1.01898764, 0.31554213, -1.15906490, 0.54803039, -1.31841690,
-0.53693713, 0.78918582, 1.56194108, 0.06661379, 0.68188768,
0.87340751, -1.52607734, 2.63594054, 2.60934675, -5.25998297,
-1.10486630, -4.38182541, 3.11359067, 2.12479542, -6.15194294
))
print(sol_blup)## trait effect level solution
## 1 1 1 1 -0.09929111
## 2 1 1 2 1.82312339
## 3 1 1 3 2.56875029
## 4 1 2 1 2.43588849
## 5 1 2 2 0.65341111
## 6 1 3 1 1.01898764
## 7 1 4 1 0.31554213
## 8 1 4 2 -1.15906490
## 9 1 4 3 0.54803039
## 10 1 4 4 -1.31841690
## 11 1 4 5 -0.53693713
## 12 1 4 6 0.78918582
## 13 1 4 7 1.56194108
## 14 1 4 8 0.06661379
## 15 1 4 9 0.68188768
## 16 1 4 10 0.87340751
## 17 1 4 11 -1.52607734
## 18 1 4 12 2.63594054
## 19 1 4 13 2.60934675
## 20 1 4 14 -5.25998297
## 21 1 4 15 -1.10486630
## 22 1 4 16 -4.38182541
## 23 1 4 17 3.11359067
## 24 1 4 18 2.12479542
## 25 1 4 19 -6.15194294As soluções do programa BLUPF90+ possuem a ordem:
[vetor de efeitos fixos], [valores genéticos],
[efeitos dos UPG].
Vamos colocar na mesma ordem das soluções que obtivemos no programa R,
para compararmos.
Você pode encontrar o ID original do animal no arquivo
renadd04.ped e corresponder com o arquivo
solutions para identificar a ordem.
sol_blup$order <- c("A", "B", "C", "S1", "S2", "cov",
"ID015", "ID006", "ID007", "ID008", "ID010", "ID009", "ID011", "ID012",
"ID013", "ID014", "ID005", "ID002", "ID004", "ID001", "ID003",
"g1", "g2", "g3", "g4")
R_blup_order <- c("A", "B", "C", "S1", "S2", "cov",
"g1", "g2", "g3", "g4",
"ID001", "ID002", "ID003", "ID004", "ID005", "ID006", "ID007", "ID008",
"ID009", "ID010", "ID011", "ID012", "ID013", "ID014", "ID015")Reordenando
sol_blup$order <- factor(sol_blup$order, levels = R_blup_order)
sol_blup2 <- sol_blup %>% arrange(order)Soluções finais (BLUPF90+)
solutions_blup <- as.matrix(sol_blup2$solution)
print(solutions_blup)## [,1]
## [1,] -0.09929111
## [2,] 1.82312339
## [3,] 2.56875029
## [4,] 2.43588849
## [5,] 0.65341111
## [6,] 1.01898764
## [7,] -4.38182541
## [8,] 3.11359067
## [9,] 2.12479542
## [10,] -6.15194294
## [11,] -5.25998297
## [12,] 2.63594054
## [13,] -1.10486630
## [14,] 2.60934675
## [15,] -1.52607734
## [16,] -1.15906490
## [17,] 0.54803039
## [18,] -1.31841690
## [19,] 0.78918582
## [20,] -0.53693713
## [21,] 1.56194108
## [22,] 0.06661379
## [23,] 0.68188768
## [24,] 0.87340751
## [25,] 0.31554213Separando os vetores
# Efeitos fixos
beta_hat <- solutions_blup[1:ncol(X)]; b <- as.matrix(beta_hat)
print(beta_hat)## [1] -0.09929111 1.82312339 2.56875029 2.43588849 0.65341111 1.01898764# Efeitos UPG
upg_hat <- solutions_blup[(ncol(X) + 1):(ncol(X) + ncol(Q))]; g <- as.matrix(upg_hat)
print(beta_hat)## [1] -0.09929111 1.82312339 2.56875029 2.43588849 0.65341111 1.01898764# Valores genéticos
animal_hat <- solutions_blup[(ncol(X) + ncol(Q) + 1):(ncol(X) + ncol(Q) + ncol(Z))]; Qg_a <- as.matrix(animal_hat)
print(beta_hat)## [1] -0.09929111 1.82312339 2.56875029 2.43588849 0.65341111 1.01898764As soluções são diferentes das obtidas diretamente,
pois utilizamos a inversa generalizada.
Ou seja, teremos diferentes soluções, dependendo de qual inversa
usarmos.
Contudo, o nosso y predito deve ser igual em ambas
predições, pois funções estimáveis são invariantes a escolha da
inversa generalizada (isso é tópico para outro post).
Vamos obter o vetor y predito, tal qual fizemos para as
soluções R.
rhs1 <- t(X) %*% X %*% b + t(X) %*% Z %*% Qg_a
rhs3 <- t(Z) %*% X %*% b - Ainv %*% Q %*% (alpha * g) + (t(Z) %*% Z + Ainv * alpha) %*% Qg_a\(My=c\)
Matriz M
M <- rbind(t(X), t(Z))Vetor c
c <- rbind(rhs1, rhs3)Resolvendo y
y_hat <- ginv(M) %*% c
final_blup <- data.frame(y_predito = y_hat, y_original = y)
print(final_blup)## y_predito y_original
## 1 2.999027 3
## 2 6.021592 6
## 3 5.976333 6
## 4 1.987161 2
## 5 2.991506 3
## 6 5.996045 6
## 7 4.004215 4
## 8 5.984827 6
## 9 8.044676 8
## 10 4.001232 4Vamos comparar os resultados via R e via blupf90!
final_solutions <- data.frame(
ID = data1$ID,
y_original = data1$obs,
y_predito_R = final$y_predito,
y_predito_blupf90 = round(final_blup$y_predito,3),
residual_blupf90 = round((final_blup$y_predito - data1$obs),4))
print(final_solutions)## ID y_original y_predito_R y_predito_blupf90 residual_blupf90
## 1 ID006 3 3 2.999 -0.0010
## 2 ID007 6 6 6.022 0.0216
## 3 ID008 6 6 5.976 -0.0237
## 4 ID009 2 2 1.987 -0.0128
## 5 ID010 3 3 2.992 -0.0085
## 6 ID011 6 6 5.996 -0.0040
## 7 ID012 4 4 4.004 0.0042
## 8 ID013 6 6 5.985 -0.0152
## 9 ID014 8 8 8.045 0.0447
## 10 ID015 4 4 4.001 0.0012Percebam que o y predito obtido tanto diretamente
(R) ou iterativamente (BLUPF90+) conferem com
o y original.
As soluções obtidas via BLUPF90+ não são exatas, pois foram
geradas por um processo iterativo. Contudo, o resíduo é mínimo.