# lineární modely v genetice

# Př1: změřené hodnoty u tří jedinců 5000, 6000, 7000
# odhad průměrné hodnoty (mí) pomocí lineárního modelu

# y = mi + e
# y = Xb + e   - maticový zápis
# soustava normálních rovnic
# |X´X| b = X´y
#       b = |X´X|-1 X´y     - výpočet vektoru řešení b (průměrná hodnota)

# vektor fenotypových změřených hodnot
y <- matrix(c(5000,6000,7000),3,1)
y

# designová matice (matice plánu pokusu) pro pevný efekt
X <- matrix(c(1,1,1),3,1)
X

XX <- t(X) %*% X
XX
Xy <- t(X) %*% y
Xy

b <- solve(XX) %*% Xy
b

# Př. 2: prostřeďový efekt s dvěma úrovněmi A a B
# jedinec 1 byl v prostředí A, jedinci 2 a 3 v prostředí B, B
# y = AB + e
# y = Xb + e
# |X´X| b = X´y
#       b = |X´X|-1 X´y


y
X <- matrix(c(1,0,0,
              0,1,1),3,2)
X

XX <- t(X) %*% X
XX

Xy <- t(X) %*% y
Xy

b <- solve(XX) %*% Xy
b

det(XX)  # určení determinantu matice levé strany X´X

# Př. 3  Výpočet matice aditivně genetické příbuznosti z rodokmenových dat (soubor rodokment.pdf)
#### balíček pedigreemm
#### 
library(pedigreemm)
rodokmen <- pedigree(sire=c(NA,NA,1,1,4,5),
                     dam =c(NA,NA,2,NA,3,2), label=1:6)
rodokmen

A <- getA(rodokmen)
A

A <- matrix(A,6,6)
A

# odhad plemenné hodnoty metodou BLUP AM 
# (zadání v souboru Zadání Animal model v R.pdf)

y <- matrix(c(200,170,180),3,1)
y

X <- matrix(c(1,1,1),3,1)
X

Z0 <-  diag(0,3)
Z0
Z1 <- diag(3)
Z1
Z <- cbind(Z0,Z1)   # designová matice Z pro náhodný efekt jedinců
Z


rodokmen <- pedigree(sire=c(NA,NA,NA,1,3,3),
                     dam =c(NA,NA,NA,2,2,4), label=1:6)
rodokmen

A <- getA(rodokmen)
A

A <- matrix(A,6,6)
A

h2 <- 1/3
h2

k <- (1-h2)/h2
k


XX <- t(X) %*% X
XX


XZ <- t(X) %*% Z
XZ

ZX <- t(Z) %*% X
ZX

ZZ <- t(Z) %*%Z
ZZ
AAk <- solve(A)*k
AAk

ZZAk <- ZZ + AAk
ZZAk

LS <- rbind(cbind(XX,XZ),cbind(ZX,ZZAk))
LS

Xy <- t(X) %*% y
Xy

Zy <- t(Z) %*% y
Zy

PS <- rbind(Xy,Zy)
PS

det(LS)

bu <- solve(LS) %*%PS
bu