Odhad Var(A) a plemenných hodnot
• ANOVA odhaduje komponenty variance při jednom typu
příbuznosti a vybalancovaných datech
• Reálné soubory jsou s nevybalancovanými daty a širokou
rodokmenovou strukturou
• Téměř většina šlechtění zvířat je založena na modelech:
• REML (restricted maximum likelihood) pro odhad variancí
• BLUP (best linear unbiased predictors) pro předpověď
plemenných hodnot
Smíšený model
y, X , Z – zaznamenané hodnoty
Odhad pevných efektů b
Odhad náhodných efektů u, e
Příklad Sire modelu
• Chceme odhadnout plemenné hodnoty tří otců (sire),
každý byl pářen s náhodnými matkami (dam) a každý
měl dva potomky, vyvíjející se ve dvou různých
prostředích (stájích).
Pozorování y Otec Stáj
y111 9 1 1
y121 12 1 2
y211 11 2 1
y212 6 2 1
y311 7 3 1
y321 14 3 2
Základní model
• Vektory a matice:
v programu R
y <- matrix(c(9,12,11,6,7,14),6,1)
X <- matrix(c(1,0,1,1,1,0,
0,1,0,0,0,1),6,2)
Z <- matrix(c(1,1,0,0,0,0,
0,0,1,1,0,0,
0,0,0,0,1,1),6,3)
Vytvoří sloupcový vektor y, s užitkovostmi dcer, s 6ti
řádky a 1 sloupcem
Vytvoří designovou matici X6,2 (pro pevný efekt
stáje)
Vytvoří designovou matici Z6,3 (pro náhodný efekt
otce)
Průměry a variance pro y = Xb + Zu + e
• Průměry:
– E(u) = E(e) = 0
– E(y) = Xb
• Variance:
– R je VCV matice pro rezidua (prostředí); předpoklad že
R = σ2
e I
– G je VCV matice pro plemenné hodnoty
– VCV matice pro y je: V = ZGZ` + R
Odhady pevných efektů a předpovědi
náhodných efektů
• Ve smíšeném modelu jsou pozorovány y, X a Z
• b, u, R a G jsou obecně neznámé
• Provádí se dva současné odhady:
BLUE pro pevné efekty b:
BLUP pro náhodné efekty u:
(Henderson, 1963)
V = ZGZ` + R
• Předpokládáme, že rezidua nejsou korelována a R = σ2
e I
(6×6)
– σ2
e = 6 -> R = 6I
• VCV matice G: předpoklad, že otcové jsou nepříbuzní a G
je diagonální matice (3×3) s prvky σ2
G = variance otců,
kde
σ2
G = σ2
A /4
– σ2
A = 8 -> G = 8/4*I
• V = ZGZ` + R
v programu R
I3 <- diag(c(1,1,1))
G <- 8/4*I3
R <- 6* diag(c(1,1,1,1,1,1))
V <- Z%*%G%*%t(Z) + R
invV <- solve(V)
Vytvoří jednotkovou matici 3×3 (~ 3 otci)
Vytvoří genetickou variančně kovarianční matici G3,3
Vytvoří prostřeďovou variančně kovarianční matici R6,6
Vytvoří fenotypovou variančně kovarianční matici V6,6
Vytvoří inverzi matice V
b <- solve(t(X)%*%invV%*%X) %*% (t(X)%*%invV%*%y)
u <- G%*%t(Z)%*%invV %*% (y-(X%*%b))
v programu R
v programu R
• u ~ (0, G), e (0, R), cov(u, e) = 0
Odvozená soustava normálních rovnic smíšeného modelu (Henderson):
v programu R
XRX <- t(X)%*%solve(R)%*%X
XRy <- t(X)%*%solve(R)%*%y
XRZ <- t(X)%*%solve(R)%*%Z
ZRX <- t(Z)%*%solve(R)%*%X
ZRZG <- t(Z)%*%solve(R)%*%Z + solve(G)
ZRy <- t(Z)%*%solve(R)%*%y
v programu R
LS1 <- cbind(XRX, XRZ)
LS2 <- cbind(ZRX, ZRZG)
LS <- rbind(LS1, LS2)
PS <- rbind(XRy, ZRy)
Vytvoří velkou matici levé strany (LS)
Vytvoří velkou matici pravé strany (PS)
bu <- solve(LS)%*%PS
Vektor řešení bu
Animal Model – viz přednáška č. 11
Předpokládáme, že naměřená užitkovost krávy je ovlivněna jen stádem, ve kterém je
chována, věkem (tj. pořadím laktace) a genotypem (tj. jedincem se svou jedinečnou
genetickou výbavou).
jedince stádo laktace užitkovost
1 1 1 4500
2* 1 1 5000
3 1 2 6500
4 2 2 8000
5* 2 1 7000
V 1. stádě jsou tři dojnice, z toho dvě jsou na první laktaci
a jedna na druhé laktaci. Ve 2. stádě jsou dvě krávy, jedna
na první a dvě na druhé laktaci. Podle původu víme, že
dojnice č. 2 a 5 mají společného otce* – jsou tedy
polosestry. Jiné příbuzenské vztahy nejsou známy.
V populaci byl odhadnuta hodnota h2 = 0,25.
modelová rovnice: yijkl = Si + Lj + uk + eijkl
maticový zápis: y = Xb +Zu + e
Odvozená soustava normálních rovnic smíšeného modelu:
Řešení
• Matice A
Jedince Otec
1 0
2* 6
3 0
4 0
5* 6
aii = 1 + 0,5(aom)
aij = 0,5(ajo + ajm)
1 2 3 4 5 otec
1 1 0 0 0 0 0
2 0 1 0 0 0,25 0,5
3 0 0 1 0 0 0
4 0 0 0 1 0 0
5 0 0,25 0 0 1 0,5
otec 0 0,5 0 0 0,5 1
Designová matice X
X1 <- matrix(c(1,1,1,0,0,0,0,0,1,1),5)
X1
[,1] [,2]
[1,] 1 0
[2,] 1 0
[3,] 1 0
[4,] 0 1
[5,] 0 1
X2 <- matrix(c(1,1,0,0,1,0,0,1,1,0),5)
X2
[,1] [,2]
[1,] 1 0
[2,] 1 0
[3,] 0 1
[4,] 0 1
[5,] 1 0
X <- cbind(X1,X2)
X
[,1] [,2] [,3] [,4]
[1,] 1 0 1 0
[2,] 1 0 1 0
[3,] 1 0 0 1
[4,] 0 1 0 1
[5,] 0 1 1 0
Pro pevný efekt stádo
Pro pevný efekt pořadí
laktace
Spojení do
jedné
designové
matice X
Matice aditivně genetické příbuznosti A
A <- matrix(c(1,0,0,0,0,0, 0,1,0,0,0.25,0.5, 0,0,1,0,0,0,
0,0,0,1,0,0, 0,0.25,0,0,1,0.5, 0,0.5,0,0,0.5,1),6)
A
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0.00 0 0 0.00 0.0
[2,] 0 1.00 0 0 0.25 0.5
[3,] 0 0.00 1 0 0.00 0.0
[4,] 0 0.00 0 1 0.00 0.0
[5,] 0 0.25 0 0 1.00 0.5
[6,] 0 0.50 0 0 0.50 1.0
Vektor užitkovostí y, designová matice Z
y <- matrix(c(4500,5000,6500,8000,7000),5,1)
y
[,1]
[1,] 4500
[2,] 5000
[3,] 6500
[4,] 8000
[5,] 7000
h2 <- 0.25
K <- (1-h2)/h2
K
[1] 3
Z <- diag(1,5)
Z
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 0
[2,] 0 1 0 0 0
[3,] 0 0 1 0 0
[4,] 0 0 0 1 0
[5,] 0 0 0 0 1
Vytvoření matice soustavy normálních
rovnic
XX <- t(X)%*%X
XZ <- t(X)%*%Z
ZX <- t(Z)%*%X
ZZ <- t(Z)%*%Z
AK <- round(K*solve(A)) současně i zaokrouhlí
/* Abychom mohli spojit matice ZZ(5x5) a AK (6x6) musíme přidat řádek a sloupec nul do matice ZZ, aby
vznikla matice o rozměrech 6x6
Stejně musíme upravit i matice XZ (+ 1 sloupec nul) a ZX (+ 1 řádek nul)
Přidáním nul se nic nemění – jen je pak možné spojit tyto submatice do matice levé strany LS */
Z0 <- matrix(c(0,0,0,0,0))
ZZ1 <- cbind(ZZ,Z0)
Z01 <- matrix(c(0,0,0,0,0,0),1,6)
ZZ2 <- rbind(ZZ1,Z01)
XZ0 <- cbind(XZ,matrix(c(0,0,0,0)))
ZX0 <- rbind(ZX,matrix(c(0,0,0,0),1,4))
ZZAK <- ZZ2+AK
Spojení dílčích matic do jedné matice levé
strany (LS)
LS1 <- cbind(XX,XZ0)
LS2 <- cbind(ZX0,ZZAK)
LS <- rbind(LS1,LS2)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 3 0 2 1 1 1 1 0 0 0
[2,] 0 2 1 1 0 0 0 1 1 0
[3,] 2 1 3 0 1 1 0 0 1 0
[4,] 1 1 0 2 0 0 1 1 0 0
[5,] 1 0 1 0 4 0 0 0 0 0
[6,] 1 0 1 0 0 5 0 0 0 -2
[7,] 1 0 0 1 0 0 4 0 0 0
[8,] 0 1 0 1 0 0 0 4 0 0
[9,] 0 1 1 0 0 0 0 0 5 -2
[10,] 0 0 0 0 0 -2 0 0 -2 5
Konstrukce matice pravé strany (PS)
Xy <- t(X)%*%y
Zy <- t(Z)%*%y
PS <- rbind(Xy,Zy,matrix(c(0)))
/* rovněž u matice pravé strany musíme > PS přidat nulu, aby
vznikl vektor o 10 řádcích */
[,1]
[1,] 16000
[2,] 15000
[3,] 16500
[4,] 14500
[5,] 4500
[6,] 5000
[7,] 6500
[8,] 8000
[9,] 7000
[10,] 0
Určení determinantu LS a zobecněná
inverze
det <- round(det(LS))
> det
[1] 0
> bu <- solve(LS)%*%PS
Error in solve.default(LS) :
system is computationally singular: reciprocal condition number = 1.33628e-17
/* Protože determinant matice LS je roven nule, je tato matice singulární a nelze ji
invertovat -> jedním z řešení je použít zobecněnou inverzi…
Je nutné si nahrát balíček MASS z nabídky: Packages -> Load Packages */
bu <- ginv(LS)%*%PS
> bu
[,1]
[1,] 2302.23138 ~ odhadnutá odchylka stáda 1
[2,] 4229.74702 ~ odhadnutá odchylka stáda 2 (odchylka 1927)
[3,] 2547.96760 ~ odhadnutá odchylka 1. laktace
[4,] 3984.01080 ~ odhadnutá odchylka 2. laktace (odchylka 1436)
[5,] -87.54974 ~ OPH krávy č. 1
[6,] 47.47015 ~ OPH krávy č. 2
[7,] 53.43945 ~ OPH krávy č. 3
[8,] -53.43945 ~ OPH krávy č. 4
[9,] 61.96703 ~ OPH krávy č. 5
[10,] 43.77487 ~ OPH otce krav č. 2 a 5
Závěry:
• Stáda se liší v chovatelské péči o 1972 kg mléka, druhá laktace převyšuje první o 1436 kg mléka.
Nejlepší kráva je č. 5 (OPH = +62 kg) a nejhorší je kráva č. 1 (OPH = -88 kg). Genetický rozdíl mezi
nimi je 150kg mléka.
• Kráva č. 4 je druhá nejhorší s plemennou hodnotou -53 kg mléka, přestože v rámci ledovaného
souboru dosahuje nejvyšší užitkovost (8000 kg mléka). Při pozornějším ledování však zjistíme, že je
na druhé laktaci, tzn., že její vysoká užitkovost je dána vyšším stupněm tělesné dospělosti (+ 1436
kg mléka) a je ve stádě s lepší chovatelskou péčí (+ 1927 kg mléka). Jestliže o tyto položky, které
jsou dány technikou chovu, se praví její užitkovost, dostane se na podprůměrnou úroveň.
• Naopak její vrstevnice – kráva č. 5 – je na první laktaci a na druhé laktaci lze tedy u ní očekávat
užitkovost 7000 + 1436 = 8436 kg mléka. Kráva č. 5 je proto po korekci +436 kg mléka lepší než
kráva č. 4, což činí v plemenné hodnotě rozdíl (v odhadu rozdílu genetického založení) 115 kg
mléka (62 + 53).
• U krav č. 2 a č. 5 jsou při odhadu plemenné hodnoty využity vlastní užitkovosti zároveň vzájemný
příbuzenský vztah zásluhou společného otce. Jejich plemenné hodnoty jsou proto stanoveny
přesněji než u ostatních krav. Plemenná hodnota otce e stanovena na základě užitkovosti těchto
dcer a činí +44 kg mléka.
• Jak ukazuje příklad, nelze se při výběru do plemenitby řídit naměřenou užitkovostí, neboť ta je
ovlivněna několika činiteli.
• Na základě odhadu plemenných hodnot dáme přednost zařazení do plemenitby krávám podle
tohoto pořadí:
Rozdíly v užitkovostech působené chovatelským prostředím jsou mnohem větší, než genetické
rozdíly mezi zvířaty. Naměřená užitkovost je ovlivněna větším počtem významných faktorů, a
proto jsou soustavy rovnic složitější a zahrnují více efektů.
Porovnejte s příkladem na konci přednášky č. 9
pořadí kráva OPH užitkovost
1 5 +62 7000
2 3 +53 6500
3 2 +47 5000
4 4 -53 8000
5 1 -88 4500