####=========================================####
#### For CRAN time limitations most lines in the
#### examples are silenced with one '#' mark,
#### remove them and run the examples
####=========================================####
data(Technow_data)
A.flint <- Technow_data$AF # Additive relationship matrix Flint
A.dent <- Technow_data$AD # Additive relationship matrix Dent
M.flint <- Technow_data$MF # Marker matrix Flint
M.dent <- Technow_data$MD # Marker matrix Dent
pheno <- Technow_data$pheno # phenotypes for 1254 single cross hybrids
pheno$hy <- paste(pheno$dent, pheno$flint, sep=":");head(pheno);dim(pheno)
####=========================================####
#### CREATE A DATA FRAME WITH ALL POSSIBLE HYBRIDS
####=========================================####
#DD <- kronecker(A.dent,A.flint,make.dimnames=TRUE)
#hybs <- data.frame(sca=rownames(DD),yield=NA,matter=NA,gcad=NA, gcaf=NA)
#hybs$yield[match(pheno$hy, hybs$sca)] <- pheno$GY
#hybs$matter[match(pheno$hy, hybs$sca)] <- pheno$GM
#hybs$gcad <- as.factor(gsub(":.*","",hybs$sca))
#hybs$gcaf <- as.factor(gsub(".*:","",hybs$sca))
#head(hybs)
#dim(hybs)
####=========================================####
#### CREATE INCIDENCE MATRICES
####=========================================####
#Z1 <- model.matrix(~gcad-1, data=hybs); dim(Z1)
#Z2 <- model.matrix(~gcaf-1, data=hybs); dim(Z2)
#Z3 <- model.matrix(~sca-1, data=hybs); dim(Z3) #only full model
####=========================================####
#### SORT INCIDENCE MATRICES ACCORDING TO RELATIONSHIP MATRICES, REAL ORDERS
####=========================================####
#real1 <- match( colnames(A.dent), gsub("gcad","",colnames(Z1)))
#real2 <- match( colnames(A.flint), gsub("gcaf","",colnames(Z2)))
#real3 <- match( colnames(DD), gsub("sca","",colnames(Z3)))
#Z1 <- Z1[,real1]
#Z2 <- Z2[,real2]
#Z3 <- Z3[,real3]
####=========================================####
#### RUN THE PREDICTION MODEL
####=========================================####
#ETA2 <- list(GCA1=list(Z=Z1, K=A.dent), GCA2=list(Z=Z2, K=A.flint))
#anss2 <- mmer(y=hybs$yield, Z=ETA2, method="AI")
#summary(anss2)
#### try EM algorithm as well if only 2 variance components
#### you can try adding SCA effects by adding SCA=list(Z=Z3, K=DD)
#### to your Z parameter
Run the code above in your browser using DataLab