Příklad řešení soustavy rovnic > A <- matrix(c(3,1,4,6),2,2) > A [,1] [,2] [1,] 3 4 [2,] 1 6 > y <- matrix(c(4,2),2,1) > y [,1] [1,] 4 [2,] 2 > x <- solve(A) %*% y > x [,1] [1,] 1.1428571 [2,] 0.1428571 Otcovský (sire) model > X <- matrix(c(1,0,1,1,1,0, + 0,1,0,0,0,1),6,2) > X [,1] [,2] [1,] 1 0 [2,] 0 1 [3,] 1 0 [4,] 1 0 [5,] 1 0 [6,] 0 1 > Z <- matrix(c(1,1,0,0,0,0, + 0,0,1,1,0,0, + 0,0,0,0,1,1),6,3) > Z [,1] [,2] [,3] [1,] 1 0 0 [2,] 1 0 0 [3,] 0 1 0 [4,] 0 1 0 [5,] 0 0 1 [6,] 0 0 1 > R <- diag(c(1,1,1,1,1,1)*6 + R <- diag(c(1,1,1,1,1,1)*6> > R <- diag(c(1,1,1,1,1,1))*6 > R [,1] [,2] [,3] [,4] [,5] [,6] [1,] 6 0 0 0 0 0 [2,] 0 6 0 0 0 0 [3,] 0 0 6 0 0 0 [4,] 0 0 0 6 0 0 [5,] 0 0 0 0 6 0 [6,] 0 0 0 0 0 6 > Va <- 8 > Vo <- Va/4 > G <- diag(3)*Vo > G [,1] [,2] [,3] [1,] 2 0 0 [2,] 0 2 0 [3,] 0 0 2 > V <- Z %*% G%*%t(Z) + R > V [,1] [,2] [,3] [,4] [,5] [,6] [1,] 8 2 0 0 0 0 [2,] 2 8 0 0 0 0 [3,] 0 0 8 2 0 0 [4,] 0 0 2 8 0 0 [5,] 0 0 0 0 8 2 [6,] 0 0 0 0 2 8 > y <-matrix(c(9,12,11,6,7,14),6,1) > y [,1] [1,] 9 [2,] 12 [3,] 11 [4,] 6 [5,] 7 [6,] 14 > BLUE <- solve(t(X)%*%solve(V)%*%X)%*%(t(X)%*%solve(V)%*%y) > BLUE [,1] [1,] 8.222222 [2,] 13.055556 > BLUP <- G%*%t(Z)%*%solve(V) %*%(y-(X%*%BLUE)) > BLUP [,1] [1,] -0.05555556 [2,] 0.11111111 [3,] -0.05555556 > # metodou normálních rovnic smíšeného modelu > XRX <- t(X)%*%solve(R)%*%X > XRX [,1] [,2] [1,] 0.6666667 0.0000000 [2,] 0.0000000 0.3333333 > XRZ <- t(X)%*%solve(R)%*%Z > XRZ [,1] [,2] [,3] [1,] 0.1666667 0.3333333 0.1666667 [2,] 0.1666667 0.0000000 0.1666667 > ZRX <- t(Z)%*%solve(R)%*%X > ZRX [,1] [,2] [1,] 0.1666667 0.1666667 [2,] 0.3333333 0.0000000 [3,] 0.1666667 0.1666667 > ZRZG <- t(Z)%*%solve(R)%*%Z + solve(G) > ZRZG [,1] [,2] [,3] [1,] 0.8333333 0.0000000 0.0000000 [2,] 0.0000000 0.8333333 0.0000000 [3,] 0.0000000 0.0000000 0.8333333 > LS1 <-rbind(XRX,ZRX) > LS1 [,1] [,2] [1,] 0.6666667 0.0000000 [2,] 0.0000000 0.3333333 [3,] 0.1666667 0.1666667 [4,] 0.3333333 0.0000000 [5,] 0.1666667 0.1666667 > LS2 <- rbind(XRZ,ZRZG) > LS2 [,1] [,2] [,3] [1,] 0.1666667 0.3333333 0.1666667 [2,] 0.1666667 0.0000000 0.1666667 [3,] 0.8333333 0.0000000 0.0000000 [4,] 0.0000000 0.8333333 0.0000000 [5,] 0.0000000 0.0000000 0.8333333 > LS <- cbind(LS1,LS2) > LS [,1] [,2] [,3] [,4] [,5] [1,] 0.6666667 0.0000000 0.1666667 0.3333333 0.1666667 [2,] 0.0000000 0.3333333 0.1666667 0.0000000 0.1666667 [3,] 0.1666667 0.1666667 0.8333333 0.0000000 0.0000000 [4,] 0.3333333 0.0000000 0.0000000 0.8333333 0.0000000 [5,] 0.1666667 0.1666667 0.0000000 0.0000000 0.8333333 > XRy <- t(X)%*%solve(R)%*%y > XRy [,1] [1,] 5.500000 [2,] 4.333333 > ZRy <- t(Z)%*%solve(R)%*%y > ZRy [,1] [1,] 3.500000 [2,] 2.833333 [3,] 3.500000 > PS <- rbind(XRy,ZRy) > PS [,1] [1,] 5.500000 [2,] 4.333333 [3,] 3.500000 [4,] 2.833333 [5,] 3.500000 > bu <- solve(LS)%*%PS > bu [,1] [1,] 8.22222222 [2,] 13.05555556 [3,] -0.05555556 [4,] 0.11111111 [5,] -0.05555556 Animal model > h2 <- 0.25 > K <- (1-h2)/h2 > 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) > 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 > y <- matrix(c(4500,5000,6500,8000,7000),5,1) > y [,1] [1,] 4500 [2,] 5000 [3,] 6500 [4,] 8000 [5,] 7000 > 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) > 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 <-(tX)%*%X Error: object 'tX' not found > 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 > 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 > AK <- solve(A)*K > AK [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3 0 0 0 0.000000e+00 0 [2,] 0 4 0 0 1.665335e-16 -2 [3,] 0 0 3 0 0.000000e+00 0 [4,] 0 0 0 3 0.000000e+00 0 [5,] 0 0 0 0 4.000000e+00 -2 [6,] 0 -2 0 0 -2.000000e+00 5 > 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 > Z0 <- matrix(c(0,0,0,0,0),1,5) > Z0 [,1] [,2] [,3] [,4] [,5] [1,] 0 0 0 0 0 > Z00 <- matrix(c(0,0,0,0,0,0),6,1) > Z00 [,1] [1,] 0 [2,] 0 [3,] 0 [4,] 0 [5,] 0 [6,] 0 > ZZZ <- cbind(rbind(ZZ,Z0),Z00) > 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 > ZAK <- ZZZ+AK > ZAK [,1] [,2] [,3] [,4] [,5] [,6] [1,] 4 0 0 0 0.000000e+00 0 [2,] 0 5 0 0 1.665335e-16 -2 [3,] 0 0 4 0 0.000000e+00 0 [4,] 0 0 0 4 0.000000e+00 0 [5,] 0 0 0 0 5.000000e+00 -2 [6,] 0 -2 0 0 -2.000000e+00 5 > 01 <- matrix(c(0,0,0,0),4,1) Error in 1 <- matrix(c(0, 0, 0, 0), 4, 1) : invalid (do_set) left-hand side to assignment > X01 <- matrix(c(0,0,0,0),4,1) > XZ0 <- cbind(XZ,X01) > XZ0 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 1 1 0 0 0 [2,] 0 0 0 1 1 0 [3,] 1 1 0 0 1 0 [4,] 0 0 1 1 0 0 > Z01 <- matrix(c(0,0,0,0),1,4) > ZX0 <- rbind(ZX,Z01) > 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 > LS1 <- rbind(XX,ZX0) > LS1 <- rbind(XZ0,ZAK) > LS1 <- rbind(XX,ZX0) > LS2 <- rbind(XZ0,ZAK) > LS <- cbind(LS1,LS2) > LS [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 3 0 2 1 1 1 1 0 0.000000e+00 0 [2,] 0 2 1 1 0 0 0 1 1.000000e+00 0 [3,] 2 1 3 0 1 1 0 0 1.000000e+00 0 [4,] 1 1 0 2 0 0 1 1 0.000000e+00 0 [5,] 1 0 1 0 4 0 0 0 0.000000e+00 0 [6,] 1 0 1 0 0 5 0 0 1.665335e-16 -2 [7,] 1 0 0 1 0 0 4 0 0.000000e+00 0 [8,] 0 1 0 1 0 0 0 4 0.000000e+00 0 [9,] 0 1 1 0 0 0 0 0 5.000000e+00 -2 [10,] 0 0 0 0 0 -2 0 0 -2.000000e+00 5 > 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,0) > PS [,1] [1,] 16000 [2,] 15000 [3,] 16500 [4,] 14500 [5,] 4500 [6,] 5000 [7,] 6500 [8,] 8000 [9,] 7000 [10,] 0 > invLS <- solve(LS) Error in solve.default(LS) : system is computationally singular: reciprocal condition number = 1.23358e-17 Protože matice LS je singulární (má nulový determinant), nelze ji invertovat standardním způsobem. Řešením je úprava matice LS, nebo použít zobecněnou inverzi. Zde je potřeba načíst balíček MASS, aby bylo možné použít příkaz pro zobecněnou inverzi. > bu <- ginv(LS)%*%PS > 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 >