#######################
# diallels by Griffing (1956)
######################

data <- readRDS("pheno2") 
head(data)
str(data)
dim(data)

# atribute identification for parents
data$female <- as.character(data$female)
data$male <- as.character(data$male)

# replacing NAs by self generations
data[data$type == "parent", 2] <- data[data$type == "parent", 1]
data[data$type == "parent", 3] <- data[data$type == "parent", 1]
data$female <- as.factor(data$female)
data$male <- as.factor(data$male)
head(data)

data$female <- as.numeric(as.character(data$female))
data$male <- as.numeric(as.character(data$male))

# test for presence of parents
(data$parent <- data$female == data$male)
data$female <- as.factor(data$female)
data$male <- as.factor(data$male)

# hybrids and their parents combination
data$combination <- data$gid


# install.packages("EMMREML")
require(EMMREML)

########################################
# Griffing Method 2 - parents and F1's
#######################################
dim(data)
str(data)

levels(data$female)
levels(data$male)
levels(data$combination)

# phenotypes
Y <- data[,"MPSadj"]

# mean + replicates + N
X <- model.matrix(~ rep + N, data = data)
head(X)
dim(X)

# Z for parents
Zp <- model.matrix(~ female -1, data = data) + model.matrix(~ male -1, data = data)
dim(Zp) # 210  x 14


# Z for hybrids and selfs
Zc <- model.matrix(~ combination -1, data = data)
dim(Zc) # 210 x 55


# let`s create the kinship   
M  <- readRDS("M")

library(snpReady)
G <- G.matrix(M, method = "VanRaden", format = "wide", plot = F)

# kinship between parents and hybrids
Ga <- G$Ga[1:14, 1:14]
dim(Ga) # 14 x 14
Ga[1:4,1:4]
Gd <- G$Gd
dim(Gd) # 55 x 55
Gd[1:4,1:4]


# Full Z
colnames(Zp) <- colnames(Ga)
colnames(Zc) <- colnames(Gd)
Z <- cbind(Zp, Zc)
dim(Z) # 210 x 69

# kinship for non-related
Ki <- diag(69)
colnames(Ki) <- c(colnames(Ga), colnames(Gd))


# single kernel with K = diag, # parents non-related
sol.single.I.2 <- emmreml(y = as.matrix(Y), X = X, Z = Z, K = Ki, varbetahat = T,varuhat = T, PEVuhat = T, test = T)
sol.single.I.2


# genomic kinships for parents and hybrids
sol.multi.G.2 <- emmremlMultiKernel(y = as.matrix(Y), X = X, Z = list(Zp, Zc), K = list(Ga, Gd), varbetahat = T, varuhat = T, PEVuhat = T, test = T)
sol.multi.G.2

M2 <- data.frame(Single.I = sol.single.I.2$uhat,
                   Multi.G = sol.multi.G.2$uhat)

round(cor(M2), 2)


######################################
# Griffing Method 4 - Just F1's
######################################    
head(data)
data <- data[data$parent == FALSE,]
dim(data) # 154 x 
data$combination <- droplevels(data$combination)
str(data)

#phenotypes
Y <- data[,"MPSadj"]

# mean + replicates + N
X <- model.matrix(~ rep + N, data = data)
head(X)
dim(X)

levels(data$female)
levels(data$male)
levels(data$combination)


# single kernel
Z <- cbind((model.matrix(~ female -1, data = data) + model.matrix(~ male -1, data = data)), model.matrix(~combination -1, data = data))
dim(Z) # 154 + 55

#multikernel
p <- length(unique(c(data$female, data$male)))
c <- ncol(Z) - p
Zp <- as.matrix(Z[, 1:p])
dim(Zp)
Zc <- Z[,(p + 1):ncol(Z)]
dim(Zc)

colnames(Zp) <- colnames(Ga)
colnames(Zc) <- colnames(Gd[(p + 1):ncol(Z), (p + 1):ncol(Z)])

Ki <- diag(55)
colnames(Ki) <- colnames(Gd)

Gd.sc <- Gd[(p + 1):ncol(Z), (p + 1):ncol(Z)]
dim(Gd.sc)  

# single kernel with K = diag
sol_single.I.4 <- emmreml(y = as.matrix(Y), X = X, Z = Z, K = Ki, varbetahat = T,varuhat = T, PEVuhat = T, test = T)
sol_single.I.4

# multi kernel
sol_multi.G.4 <- emmremlMultiKernel(y = as.matrix(Y), X = X, Z = list(Zp, Zc), K = list(Ga, Gd.sc), varbetahat = T, varuhat = T, PEVuhat = T, test = T)
sol_multi.G.4

M4 <- data.frame(Single.I = sol_single.I.4$uhat,
                   Multi.G = sol_multi.G.4$uhat)

round(cor(M4), 2)

############