AplikGenZ Příklad Animal Modelu LS2014
TGU Stránka 1
Příklad animal modelu dle přednášky č. 10
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
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).
modelová rovnice: yijkl = Si + Lj + uk + eijkl
maticový zápis: y = Xb + Zu + e
Odvozená soustava normálních rovnic smíšeného modelu:
Do okna programu R (www.r-project.cz) vkládáme jednotlivé příkazy (znak >
na začátku každého příkazu je automaticky přidáván programem):
> y <- matrix(c(4500,5000,6500,8000,7000),5,1)
> y
[,1]
[1,] 4500
[2,] 5000
[3,] 6500
[4,] 8000
[5,] 7000
> X <- matrix(c(1,1,1,0,0,
+ 0,0,0,1,1,
+ 1,1,0,0,1,
+ 0,0,1,1,0),5,4)






+′′
′′
K-1
AZZXZ
ZXXX






u
b






′
′
=
yZ
yX
AplikGenZ Příklad Animal Modelu LS2014
TGU Stránka 2
>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
> 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
> h2 = 0.25
> K = (1-h2)/h2
> K
[1] 3
> 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,6)
AplikGenZ Příklad Animal Modelu LS2014
TGU Stránka 3
> 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
> XX <- t(X)%*%X
> XX
[,1] [,2] [,3] [,4]
[1,] 3 0 2 1
[2,] 0 2 1 1
[3,] 2 1 3 0
[4,] 1 1 0 2
> XZ <- t(X)%*%Z
> XZ
[,1] [,2] [,3] [,4] [,5]
[1,] 1 1 1 0 0
[2,] 0 0 0 1 1
[3,] 1 1 0 0 1
[4,] 0 0 1 1 0
> ZX <- t(Z)%*%X
> ZX
[,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
AplikGenZ Příklad Animal Modelu LS2014
TGU Stránka 4
> ZZ <- t(Z)%*%Z
> ZZ
[,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
> Ainv <- solve(A)
> Ainv
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0.000000e+00 0 0 0.000000e+00 0.0000000
[2,] 0 1.333333e+00 0 0 5.551115e-17 -0.6666667
[3,] 0 0.000000e+00 1 0 0.000000e+00 0.0000000
[4,] 0 0.000000e+00 0 1 0.000000e+00 0.0000000
[5,] 0 2.960595e-17 0 0 1.333333e+00 -0.6666667
[6,] 0 -6.666667e-01 0 0 -6.666667e-01 1.6666667
> Ak <- Ainv*K
> Ak
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3 0.000000e+00 0 0 0.000000e+00 0
[2,] 0 4.000000e+00 0 0 1.665335e-16 -2
[3,] 0 0.000000e+00 3 0 0.000000e+00 0
[4,] 0 0.000000e+00 0 3 0.000000e+00 0
[5,] 0 8.881784e-17 0 0 4.000000e+00 -2
[6,] 0 -2.000000e+00 0 0 -2.000000e+00 5
/* 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 */
AplikGenZ Příklad Animal Modelu LS2014
TGU Stránka 5
> Z0 <- matrix(c(0,0,0,0,0))
> Z0
[,1]
[1,] 0
[2,] 0
[3,] 0
[4,] 0
[5,] 0
> ZX0 <- cbind(ZX,Z0)
> ZX0
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 1 0 0
[2,] 1 0 1 0 0
[3,] 1 0 0 1 0
[4,] 0 1 0 1 0
[5,] 0 1 1 0 0
> Z0 <- matrix(c(0,0,0,0),1,4)
> Z0
[,1] [,2] [,3] [,4]
[1,] 0 0 0 0
> ZX0 <- rbind(ZX,Z0)
> ZX0
[,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
[6,] 0 0 0 0
AplikGenZ Příklad Animal Modelu LS2014
TGU Stránka 6
> ZZ0 <- matrix(c(0,0,0,0,0),1,5)
> ZZ00 <- matrix(c(0,0,0,0,0,0),6,1)
> ZZa <- rbind(ZZ,ZZ0)
> ZZa
[,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
[6,] 0 0 0 0 0
> ZZZ <- cbind(ZZa,ZZ00)
> ZZZ
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0 0
[2,] 0 1 0 0 0 0
[3,] 0 0 1 0 0 0
[4,] 0 0 0 1 0 0
[5,] 0 0 0 0 1 0
[6,] 0 0 0 0 0 0
> ZZAK <- ZZZ + Ak
> ZZAK
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 4 0.000000e+00 0 0 0.000000e+00 0
[2,] 0 5.000000e+00 0 0 1.665335e-16 -2
[3,] 0 0.000000e+00 4 0 0.000000e+00 0
[4,] 0 0.000000e+00 0 4 0.000000e+00 0
[5,] 0 8.881784e-17 0 0 5.000000e+00 -2
[6,] 0 -2.000000e+00 0 0 -2.000000e+00 5
AplikGenZ Příklad Animal Modelu LS2014
TGU Stránka 7
> Zc <- matrix(c(0,0,0,0),4,1)
> LS1 <- cbind(XX,cbind(XZ,Zc))
> LS1
[,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
> LS2 <- cbind(ZX0,ZZAK)
> LS2
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 0 1 0 4 0.000000e+00 0 0 0.000000e+00 0
[2,] 1 0 1 0 0 5.000000e+00 0 0 1.665335e-16 -2
[3,] 1 0 0 1 0 0.000000e+00 4 0 0.000000e+00 0
[4,] 0 1 0 1 0 0.000000e+00 0 4 0.000000e+00 0
[5,] 0 1 1 0 0 8.881784e-17 0 0 5.000000e+00 -2
[6,] 0 0 0 0 0 -2.000000e+00 0 0 -2.000000e+00 5
> LS <- rbind(LS1,LS2)
> LS
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 3 0 2 1 1 1.000000e+00 1 0 0.000000e+00 0
[2,] 0 2 1 1 0 0.000000e+00 0 1 1.000000e+00 0
[3,] 2 1 3 0 1 1.000000e+00 0 0 1.000000e+00 0
[4,] 1 1 0 2 0 0.000000e+00 1 1 0.000000e+00 0
[5,] 1 0 1 0 4 0.000000e+00 0 0 0.000000e+00 0
[6,] 1 0 1 0 0 5.000000e+00 0 0 1.665335e-16 -2
[7,] 1 0 0 1 0 0.000000e+00 4 0 0.000000e+00 0
[8,] 0 1 0 1 0 0.000000e+00 0 4 0.000000e+00 0
[9,] 0 1 1 0 0 8.881784e-17 0 0 5.000000e+00 -2
[10,] 0 0 0 0 0 -2.000000e+00 0 0 -2.000000e+00 5
AplikGenZ Příklad Animal Modelu LS2014
TGU Stránka 8
> Xy <- t(X)%*%y
> Xy
[,1]
[1,] 16000
[2,] 15000
[3,] 16500
[4,] 14500
> Zy <- t(Z)%*%y
> Zy
[,1]
[1,] 4500
[2,] 5000
[3,] 6500
[4,] 8000
[5,] 7000
> PS <- rbind(Xy,Zy)
> PS
[,1]
[1,] 16000
[2,] 15000
[3,] 16500
[4,] 14500
[5,] 4500
[6,] 5000
[7,] 6500
[8,] 8000
[9,] 7000
/* rovněž u matice pravé strany musíme > PS přidat nulu, aby vznikl vektor o 10 řádcích */
> PS0 <- rbind(PS,0)
AplikGenZ Příklad Animal Modelu LS2014
TGU Stránka 9
> PS0
[,1]
[1,] 16000
[2,] 15000
[3,] 16500
[4,] 14500
[5,] 4500
[6,] 5000
[7,] 6500
[8,] 8000
[9,] 7000
[10,] 0
> LSinv <- solve(LS)
Error in solve.default(LS) :
system is computationally singular: reciprocal condition number = 1.33628e-17
> det <- round(det(LS))
> det
[1] 0
/* 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*/
AplikGenZ Příklad Animal Modelu LS2014
TGU Stránka 10
> bu <- ginv(LS)%*%PS0
> bu
[,1]
[1,] 2302.23138
[2,] 4229.74702
[3,] 2547.96760
[4,] 3984.01080
[5,] -87.54974
[6,] 47.47015
[7,] 53.43945
[8,] -53.43945
[9,] 61.96703
[10,] 43.77487
Tomáš Urban
UMFGZ MENDELU
urban@mendelu.cz