library(tidyverse)
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Warning: package 'tibble' was built under R version 3.3.3
## Warning: package 'tidyr' was built under R version 3.3.3
## Warning: package 'readr' was built under R version 3.3.3
## Warning: package 'purrr' was built under R version 3.3.3
## Conflicts with tidy packages ----------------------------------------------
## filter(): dplyr, stats
## lag(): dplyr, stats
MIVQUE
df <- data.frame(cattle=1:9,
herd=c(1,1,2,2,2,3,3,3,3),
sire=c("A","D","B","D","D","C","C","D","D"),
# sire=c(1,4,2,4,4,3,3,4,4),
yield=c(110,100,110,100,100,110,110,100,100),
y=c(139.33,85.15,122.09,145.66,132.18,130.23,122.37,92.13,82.06))
df$herd <- factor(df$herd)
df$sire <- factor(df$sire)
I <- diag(x=1,nrow = 9)
X <- model.matrix(~herd-1, data = df)
Z <- model.matrix(~sire-1, data = df)
P <- I-X%*%solve(t(X)%*%X)%*%t(X)
# ped <- data.frame(id=df$cattle, sire=df$sire, dam=rep(0,9))
ped <- data.frame(id=unique(df$sire), sire=rep(NA,4), dam=rep(NA,4))
A <- 2*kinship2::kinship(ped$id,ped$sire,ped$dam)
A2 <- expm::sqrtm(A)
Vr <- P%*%Z
# %*%A2
H11 <- psych::tr(Vr%*%t(Vr)%*%Vr%*%t(Vr))
H12 <- H21 <- psych::tr(Vr%*%t(Vr)%*%P)
H22 <- psych::tr(P)
H1 <- cbind(H11,H12);H2 <- cbind(H21,H22)
H <- rbind(H1,H2)
Qr <- t(df$y)%*%Vr%*%t(Vr)%*%df$y
Qe <- t(df$y)%*%P%*%df$y
Q <- rbind(Qr,Qe)
para <- solve(H)%*%Q %>% as.data.frame
rownames(para) <- c("sigma_a^2","sigma_e^2")
para
## V1
## sigma_a^2 451.1045
## sigma_e^2 235.2104
BLUP
H1 <- cbind(t(X)%*%X, t(X)%*%Z)
H2 <- cbind(t(Z)%*%X, t(Z)%*%Z+solve(A)/(para$V1[1]/para$V1[2]))
H <- rbind(H1,H2)
Q <- rbind(t(X)%*%df$y,t(Z)%*%df$y)
D <- solve(H);D
## herd1 herd2 herd3 sireA sireB sireC
## herd1 1.1375988 0.4503074 0.4370848 -0.7477266 -0.2959803 -0.3466987
## herd2 0.4503074 0.9430378 0.5010258 -0.2959803 -0.6198446 -0.3974172
## herd3 0.4370848 0.5010258 0.9006352 -0.2872893 -0.3293167 -0.7143901
## sireA -0.7477266 -0.2959803 -0.2872893 1.1487545 0.1945434 0.2278799
## sireB -0.2959803 -0.6198446 -0.3293167 0.1945434 1.0646996 0.2612163
## sireC -0.3466987 -0.3974172 -0.7143901 0.2278799 0.2612163 0.9632627
## sireD -0.5274709 -0.6046344 -0.5868803 0.3466987 0.3974172 0.4655176
## sireD
## herd1 -0.5274709
## herd2 -0.6046344
## herd3 -0.5868803
## sireA 0.3466987
## sireB 0.3974172
## sireC 0.4655176
## sireD 0.7082431
effect <- D%*%Q;effect
## [,1]
## herd1 110.76312
## herd2 149.96605
## herd3 106.92453
## sireA 18.77658
## sireB -18.32251
## sireC 15.36875
## sireD -15.82282
ML
# lamta <-
# H <- Z%*%A%*%t(Z)%*%lamta+I
# beta <- solve(t(X)%*%solve(H)%*%X)%*%t(X)%*%solve(H)%*%y
# sigma_e <- 1/n*t(y-X*beta)%*%solve(H)%*%(y-X*beta)