Animal Model - viz přednáška č. 11 Předpokládáme, že naměřená užitkovost krávy (y) 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). 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. jedince 1 2* 3 4 5* stádo 1 1 1 2 2 laktace užitkovost 4500 5000 6500 8000 7000 V populaci byl odhadnuta hodnota h2 = 0,25. modelová rovnice: maticový zápis: Sj + Lj + uk + eijkl Yijki y = Xb +Zu + e Odvozená soustava normálních rovnic smíšeného modelu XX X'Z Z'X Z'Z + A'/C IV x'y u _z'y_ K - \-h2 Řešení-matice aditivně genetické Matice A příbuznosti Jedince Otec 1 2* 3 4 5* 0 6 0 0 6 a„ = 1 + 0,5(aom) ay = 0,5(ajo + aim) onv 2 3 4 5 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 Programování v R/RStudiu Designová matice X XI <- matrix (c (1, 1, 1, 0, 0, 0, 0, 0, 1, 1) , 5) XI Pro pevný efekt stádo [,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 Pro pevný efekt pořadí laktace X <- cbind(Xl,X2) X [,1][,2][,3] [,4] [1,] 10 10 [2,] 10 10 [3,] 10 0 1 [4,] 0 10 1 [5,] 0 110 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] [1J 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,] 4500 [2,] 5000 [3,] 6500 z <- diag(l,5) [4,] 8000 z [5/]7000 L1][,2][,3][,4][,5] [1,] 1 0 0 0 0 h2 <- 0.25 [2,] 0 10 0 0 K <- (l-h2)/h2 [3J 0 0 10 0 *1]3 [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 */ ZO <- matrix(c (0,0,0,0,0)) ZZ1 <- cbind(ZZ,ZO) ZOl <- matrix(c(0,0,0,0,0,0),1,6) ZZ2 <- rbind(ZZ1,ZOl) 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) LSI <- cbind(XX,XZO) LS2 <- cbind(ZXO,ZZAK) LS <- rbind(LSI,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,] 101040000 0 [6,] 101005000 -2 [7,] 100100400 0 [8,] 010100040 0 [9,] 01 1000005 -2 [10,] 00000 -2 00 -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 */ [,i] [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 [110 > 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 Li] [1,] 2302.23138 [2,] 4229.74702 [3,] 2547.96760 [4,] 3984.01080 [5,] [6,] [7,] [8,] [9,] -87.54974 47.47015 53.43945 -53.43945 61.96703 [10,] 43.77487 odhadnutá odchylka stáda 1 odhadnutá odchylka stáda 2 (odchylka 1927) odhadnutá odchylka 1. laktace odhadnutá odchylka 2. laktace (odchylka 1436) ~ OPH krávy č. 1 ~ OPH krávy č. 2 ~ OPH krávy č. 3 ~ OPH krávy č. 4 ~ OPH krávy č. 5 ~ 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žitkovosti, 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ů. poř1adl T Z" 2 3 +53 6500 3 2 +47 5000 4 4 -53 8000 5 1 -88 4500