Nothing
R/fit_model.R
In qtlpoly: Random-Effect Multiple QTL Mapping in Autopolyploids
Defines functions
mmer_adapted
summary.qtlpoly.fitted
fit_model
Documented in
fit_model
mmer_adapted
summary.qtlpoly.fitted
#' Fits multiple QTL models
#'
#' Fits alternative multiple QTL models by performing variance component estimation using REML.
#'
#' @param data an object of class \code{qtlpoly.data}.
#'
#' @param model an object of class \code{qtlpoly.profile} or \code{qtlpoly.remim}.
#'
#' @param probs a character string indicating if either \code{"joint"} (genotypes) or \code{"marginal"} (parental gametes) conditional probabilities should be used.
#'
#' @param polygenes a character string indicating if either \code{"none"}, \code{"most"} or \code{"all"} QTL should be used as polygenes.
#'
#' @param keep if \code{TRUE} (default), stores all matrices and estimates from fitted model; if \code{FALSE}, nothing is stored.
#'
#' @param verbose if \code{TRUE} (default), current progress is shown; if
#' \code{FALSE}, no output is produced.
#'
#' @param object an object of class \code{qtlpoly.fitted} to be summarized.
#'
#' @param pheno.col a numeric vector with the phenotype column numbers to be summarized; if \code{NULL} (default), all phenotypes from \code{'data'} will be included.
#'
#' @param ... currently ignored
#'
#' @return An object of class \code{qtlpoly.fitted} which contains a list of \code{results} for each trait with the following components:
#'
#' \item{pheno.col}{a phenotype column number.}
#' \item{fitted}{a \pkg{sommer} object of class \code{mmer}.}
#' \item{qtls}{a data frame with information from the mapped QTL.}
#'
#' @seealso \code{\link[qtlpoly]{read_data}}, \code{\link[qtlpoly]{remim}}
#'
#' @examples
#' \donttest{
#' # Estimate conditional probabilities using mappoly package
#' library(mappoly)
#' library(qtlpoly)
#' genoprob4x = lapply(maps4x[c(5)], calc_genoprob)
#' data = read_data(ploidy = 4, geno.prob = genoprob4x, pheno = pheno4x, step = 1)
#'
#' # Search for QTL
#' remim.mod = remim(data = data, pheno.col = 1, w.size = 15, sig.fwd = 0.0011493379,
#' sig.bwd = 0.0002284465, d.sint = 1.5, n.clusters = 1)
#'
#' # Fit model
#' fitted.mod = fit_model(data=data, model=remim.mod, probs="joint", polygenes="none")
#' }
#'
#' @author Guilherme da Silva Pereira, \email{gdasilv@@ncsu.edu}
#'
#' @references
#' Covarrubias-Pazaran G (2016) Genome-assisted prediction of quantitative traits using the R package sommer. \emph{PLoS ONE} 11 (6): 1–15. \doi{10.1371/journal.pone.0156744}.
#'
#' Pereira GS, Gemenet DC, Mollinari M, Olukolu BA, Wood JC, Mosquera V, Gruneberg WJ, Khan A, Buell CR, Yencho GC, Zeng ZB (2020) Multiple QTL mapping in autopolyploids: a random-effect model approach with application in a hexaploid sweetpotato full-sib population, \emph{Genetics} 215 (3): 579-595. \doi{10.1534/genetics.120.303080}.
#'
#' @useDynLib qtlpoly, .registration = TRUE
#' @importFrom Rcpp evalCpp
#' @export fit_model
fit_model <- function(data, model, probs="joint", polygenes="none", keep=TRUE, verbose=TRUE, pheno.col = NULL) {
results <- vector("list", length(model$results))
names(results) <- names(model$results)
if(probs=="marginal") {
ploidy <- data$ploidy
Palleles <- letters[1:ploidy]
Pgametes <- lapply(combn(Palleles, ploidy/2, simplify = FALSE), paste, collapse="")
Qalleles <- letters[(ploidy+1):(2*ploidy)]
Qgametes <- lapply(combn(Qalleles, ploidy/2, simplify = FALSE), paste, collapse="")
sibs <- sapply( Pgametes, FUN=function(x) paste(x, Pgametes, sep="") )
Pia <- matrix(data = NA, nrow = length(Pgametes), ncol = length(Pgametes))
for(i in 1:ncol(sibs)) {
for(j in 1:nrow(sibs)) {
Pia[i,j] <- (ploidy-length(unique(strsplit(sibs[i,j], "")[[1]])))/(ploidy/2)
}
}
Pib <- Pia # Pi1 = Pi2
colnames(Pia) <- rownames(Pia) <- unlist(Pgametes)
colnames(Pib) <- rownames(Pib) <- unlist(Qgametes)
indnames <- dimnames(data$Z)[[3]]
mrknames <- dimnames(data$Z)[[2]]
nmrk = length(mrknames)
nind = length(indnames)
Za <- array(data = NA, dim = c(length(Pgametes), nmrk, nind), dimnames = list(c(Pgametes), c(mrknames), c(indnames))) # Za = Z1
Zb <- array(data = NA, dim = c(length(Qgametes), nmrk, nind), dimnames = list(c(Qgametes), c(mrknames), c(indnames))) # Zb = Z2
for(m in 1:nmrk) {
for(i in 1:nind) {
Za[,m,i] <- colSums(matrix(data$Z[,m,i], nrow=length(Pgametes)))
Zb[,m,i] <- rowSums(matrix(data$Z[,m,i], nrow=length(Qgametes)))
}
}
}
for(p in 1:length(model$results)) {
pheno.col <- model$results[[p]]$pheno.col
qtl.mrk <- model$results[[p]]$qtl[,"Nmrk"]
qtl.lgr <- model$results[[p]]$qtl[,"LG"]
qtl.pos <- model$results[[p]]$qtl[,"Pos"]
if(verbose) {
if(length(qtl.mrk) == 1 && verbose) cat("There is ", length(qtl.mrk), " QTL in the model for trait ", pheno.col, " ", sQuote(colnames(data$pheno)[pheno.col]), ". Fitting model... ", sep="")
if(length(qtl.mrk) >= 2 && verbose) cat("There are ", length(qtl.mrk), " QTL in the model for trait ", pheno.col, " ", sQuote(colnames(data$pheno)[pheno.col]), ". Fitting model... ", sep="")
}
if(!is.null(model$results[[p]]$qtls)) {
nqtl <- dim(model$results[[p]]$qtls)[1]
ind <- rownames(data$pheno)[which(!is.na(data$pheno[,pheno.col]))]
Y <- data$pheno[ind,pheno.col]
Zstar <- diag(length(ind))
ETA <- NULL
if(nqtl > 1) {
if(probs=="joint" && polygenes=="none") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
for(q in 1:nqtl) {
eta <- list(list(Z=t(data$Z[,markers[q],ind]),K=data$Pi))
names(eta) <- paste("g", q, sep="")
ETA <- c(ETA, eta)
}
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
qtl <- c()
for(q in 1:nqtl) {
H2 <- sum(c(fitted$sigma)[q])/sum(c(fitted$sigma))
qtl <- c(qtl, c(model$results[[p]]$qtls[q,c(1:3)], NA, c(fitted$sigma)[q], NA, H2))
}
H2 <- sum(c(fitted$sigma)[1:nqtl])/sum(c(fitted$sigma))
qtl <- c(qtl, NA, NA, NA, c(fitted$Beta), NA, c(fitted$sigma)[nqtl+1], H2)
qtls <- as.data.frame(matrix(qtl, ncol=7, byrow=TRUE), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(g)", "Var(e)", "h2")
}
if(probs=="joint" && polygenes=="most") {
fitted <- vector("list", nqtl)
qtl <- c()
for(q in 1:nqtl) {
markers <- model$results[[p]]$qtls[-q,"Nmrk"]
Gstar <- apply(data$G[ind,ind,markers], MARGIN = c(1,2), sum)/length(markers); Gstar[1:5,1:5]
m <- model$results[[p]]$qtls[q,"Nmrk"]
ETA <- list(g=list(Z=t(data$Z[,m,ind]),K=data$Pi),gstar=list(Z=Zstar,K=Gstar))
fitted[[q]] <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- c(fitted[[q]]$sigma)[1]/sum(c(fitted[[q]]$sigma))
qtl <- c(qtl, c(model$results[[p]]$qtls[q,c(1:3)], c(fitted[[q]]$Beta), c(fitted[[q]]$sigma), H2))
}
qtls <- as.data.frame(matrix(qtl, ncol=8, byrow=TRUE), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(g)", "Var(g*)", "Var(e)", "h2")
}
if(probs=="joint" && polygenes=="all") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
Gstar <- apply(data$G[ind,ind,markers], MARGIN = c(1,2), sum)/length(markers); Gstar[1:5,1:5]
ETA <- list(gstar=list(Z=Zstar,K=Gstar))
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- c(fitted$sigma)[1]/sum(c(fitted$sigma))
qtl <- c(rep(NA, 4*nqtl), c(fitted$Beta), c(fitted$sigma), H2)
qtls <- as.data.frame(cbind(rbind(model$results[[p]]$qtls[,c(1:3)], rep(NA, 3)), matrix(qtl, ncol=4, byrow=TRUE), stringsAsFactors=FALSE))
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(g*)", "Var(e)", "h2")
}
if(probs=="marginal" && polygenes=="none") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
for(q in 1:nqtl) {
eta <- list(list(Z=t(Za[,markers[q],ind]),K=Pia),list(Z=t(Zb[,markers[q],ind]),K=Pib),list(Z=t(data$Z[,markers[q],ind]),K=data$Pi))
names(eta) <- paste(c("a","b","c"), rep(q,3), sep="")
ETA <- c(ETA, eta)
}
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
qtl <- c()
for(q in 1:nqtl) {
H2 <- sum(c(fitted$sigma)[(1:3)+((q-1)*3)])/sum(c(fitted$sigma))
qtl <- c(qtl, c(model$results[[p]]$qtls[q,c(1:3)], NA, c(fitted$sigma)[(1:3)+((q-1)*3)], NA, H2))
}
H2 <- sum(c(fitted$sigma)[1:(nqtl*3)])/sum(c(fitted$sigma))
qtl <- c(qtl, NA, NA, NA, c(fitted$Beta), rep(NA, 3), c(fitted$sigma)[(nqtl*3)+1], H2)
qtls <- as.data.frame(matrix(qtl, ncol=9, byrow=TRUE), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(a)", "Var(b)", "Var(d)", "Var(e)", "h2")
}
if(probs=="marginal" && polygenes=="most") {
fitted <- vector("list", nqtl)
qtl <- c()
for(q in 1:nqtl) {
markers <- model$results[[p]]$qtls[-q,"Nmrk"]
Gstar <- apply(data$G[ind,ind,markers], MARGIN = c(1,2), sum)/length(markers); Gstar[1:5,1:5]
m <- model$results[[p]]$qtls[q,"Nmrk"]
ETA <- list(a=list(Z=t(Za[,m,ind]),K=Pia),b=list(Z=t(Zb[,m,ind]),K=Pib),d=list(Z=t(data$Z[,m,ind]),K=data$Pi),gstar=list(Z=Zstar,K=Gstar))
fitted[[q]] <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- sum(c(fitted[[q]]$sigma)[(1:3)])/sum(c(fitted[[q]]$sigma))
qtl <- c(qtl, c(model$results[[p]]$qtls[q,c(1:3)], c(fitted[[q]]$Beta), c(fitted[[q]]$sigma), H2))
}
qtls <- as.data.frame(matrix(qtl, ncol=10, byrow=TRUE), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(a)", "Var(b)", "Var(d)", "Var(g*)", "Var(e)", "h2")
}
if(probs=="marginal" && polygenes=="all") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
A <- array(data = NA, dim = c(data$nind, data$nind, nqtl))
B <- array(data = NA, dim = c(data$nind, data$nind, nqtl))
for(q in 1:nqtl) {
A[,,q] <- t(Za[,markers[q],])%*%Pia%*%Za[,markers[q],]
B[,,q] <- t(Zb[,markers[q],])%*%Pib%*%Zb[,markers[q],]
}
Astar <- apply(A[ind,ind,], MARGIN = c(1,2), sum)/length(markers); Astar[1:5,1:5]
Bstar <- apply(B[ind,ind,], MARGIN = c(1,2), sum)/length(markers); Bstar[1:5,1:5]
Dstar <- apply(data$G[ind,ind,markers], MARGIN = c(1,2), sum)/length(markers); Dstar[1:5,1:5]
ETA <- list(astar=list(Z=Zstar,K=Astar),bstar=list(Z=Zstar,K=Bstar),dstar=list(Z=Zstar,K=Dstar))
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- sum(c(fitted$sigma)[1:3])/sum(c(fitted$sigma))
qtl <- c(rep(NA, 6*nqtl), c(fitted$Beta), c(fitted$sigma), H2)
qtls <- as.data.frame(cbind(rbind(model$results[[p]]$qtls[,c(1:3)], rep(NA, 3)), matrix(qtl, ncol=6, byrow=TRUE), stringsAsFactors=FALSE))
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(a*)", "Var(b*)", "Var(d*)", "Var(e)", "h2")
}
}
if(nqtl == 1) {
if(probs=="joint") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
ETA <- list(g=list(Z=t(data$Z[,markers,ind]),K=data$Pi))
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- sum(c(fitted$sigma)[1])/sum(c(fitted$sigma))
qtl <- c(c(fitted$Beta), c(fitted$sigma), H2)
qtls <- as.data.frame(c(model$results[[p]]$qtls[,c(1:3)], matrix(qtl, ncol=4, byrow=TRUE)), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(g)", "Var(e)", "h2")
}
if(probs=="marginal") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
ETA <- list(a=list(Z=t(Za[,markers,ind]),K=Pia),b=list(Z=t(Zb[,markers,ind]),K=Pib),d=list(Z=t(data$Z[,markers,ind]),K=data$Pi))
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- sum(c(fitted$sigma)[1:3])/sum(c(fitted$sigma))
qtl <- c(c(fitted$Beta), c(fitted$sigma), H2)
qtls <- as.data.frame(c(model$results[[p]]$qtls[,c(1:3)], matrix(qtl, ncol=6, byrow=TRUE)), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(a)", "Var(b)", "Var(d)", "Var(e)", "h2")
}
}
if(verbose) cat("Done! \n\n", sep="")
}
if(is.null(model$results[[p]]$qtls)) {
if(verbose) cat("There are no QTL in the model for trait ", pheno.col, " ", sQuote(colnames(data$pheno)[pheno.col]), ". No model has been fitted! \n\n", sep="")
fitted <- NULL
qtls <- NULL
}
if(!keep) fitted <- NULL
results[[p]] <- list(
pheno.col=pheno.col,
fitted=fitted,
qtls=qtls)
}
structure(list(data=deparse(substitute(data)),
model=deparse(substitute(model)),
pheno.col=model$pheno.col,
probs=probs,
polygenes=polygenes,
results=results),
class="qtlpoly.fitted")
}
#' @rdname fit_model
#' @export
summary.qtlpoly.fitted <- function(object, pheno.col=NULL, ...) {
if(any(class(object) == "qtlpoly.fitted")) cat("This is an object of class 'qtlpoly.fitted'\n")
if(is.null(pheno.col)) {
pheno.col <- 1:length(object$results)
} else {
pheno.col <- which(object$pheno.col %in% pheno.col)
}
for(p in pheno.col) {
cat("\n* Trait", object$results[[p]]$pheno.col, sQuote(names(object$results)[[p]]), "\n")
if(!is.null(object$results[[p]]$qtls)) print(object$results[[p]]$qtls)
else cat("There are no QTL in the model \n")
}
}
#' Adapt to sommer MNR core function
#'
#' Adapts genomic incidence and relationship (varcov) matrices to run using sommer's C++ core function (v. 4.0 or higher)
#' Function adapted from sommer v. 3.6 (Author: Giovanny Covarrubias-Pazaran)
#'
#' @keywords internal
#'
#' @author Gabriel de Siqueira Gesteira, \email{gdesiqu@@ncsu.edu}
#'
#' @references
#' Covarrubias-Pazaran G (2016) Genome-assisted prediction of quantitative traits using the R package sommer. \emph{PLoS ONE} 11 (6): 1–15. \doi{10.1371/journal.pone.0156744}.
#'
#' Pereira GS, Gemenet DC, Mollinari M, Olukolu BA, Wood JC, Mosquera V, Gruneberg WJ, Khan A, Buell CR, Yencho GC, Zeng ZB (2020) Multiple QTL mapping in autopolyploids: a random-effect model approach with application in a hexaploid sweetpotato full-sib population, \emph{Genetics} 215 (3): 579-595. \doi{10.1534/genetics.120.303080}.
mmer_adapted <- function(Y,X=NULL,Z=NULL,R=NULL,W=NULL,method="NR",init=NULL,iters=20,tolpar=1e-3,
tolparinv=1e-6,draw=FALSE,silent=FALSE, constraint=TRUE,
EIGEND=FALSE, forced=NULL, IMP=FALSE, complete=TRUE,
check.model=TRUE, restrained=NULL, REML=TRUE, init.equal=TRUE,
date.warning=TRUE, verbose=FALSE){
#R=NULL
if (inherits(Y, "formula")){
stop("\nYou have to use mmer2 function for formula-based models.\n", call. = FALSE)
}
diso <- dim(as.data.frame(Y))[2]
## control for 2-level list structure
if(!is.list(Z)){
stop("Please provide the Z parameter as a 2 level list structure.\nFor example for 2 random effects 'A' and 'B' do:\n ETA <- list( A=list( Z=myZ1, K=myK1 ) , B=list( Z=myZ2, K=myK2 ) )\n mod <- mmer(Y=y, Z=ETA)\nwhere Z's and K's are the incidence and var-covar matrices respectively.\nIf any Z or K is not provided, an identity matrix will be assumed. ",call. = FALSE)
}else{
if(!is.list(Z[[1]])){
stop("Please provide the Z parameter as a 2 level list structure.\nFor example for 2 random effects 'A' and 'B' do:\n ETA <- list( A=list( Z=myZ1, K=myK1 ) , B=list( Z=myZ2, K=myK2 ) )\n mod <- mmer(Y=y, Z=ETA)\nwhere Z's and K's are the incidence and var-covar matrices respectively.\nIf any Z or K is not provided, an identity matrix will be assumed. ",call. = FALSE)
}
}
## control for Z-K names
zzzkkk <- unlist(lapply(Z,function(x){length(names(x))}))
badRE <- which(zzzkkk==0) # BAD RE WITH NO NAMES
if(length(badRE)>0){
stop("Please when specifying a random effect use the names; \n'Z' for incidence and 'K' for variance-covariance matrices.\nFor example for 1 random effect (i.e. named 'A') model do:\n ETA <- list( A=list( Z=M1, K=M2) )\n mod <- mmer(Y=y, Z=ETA)\nSpecifying at least one; Z or K. You need to specify if is a 'Z' or 'K' \nsince this is the only way the program distinguishes between matrices.",call. = FALSE)
}
#########*****************************
## make sure user don't provide the same names for random effects
his.names <- names(Z)
if(!is.null(his.names)){
badnames <- which(duplicated(his.names))
if(length(badnames)>0){
his.names[badnames] <- paste(his.names[badnames],1:(length(his.names[badnames])), sep=".")
names(Z) <- his.names
}
}
dZ <- unlist(lapply(Z,function(x){dim(x$Z)[1]}))
if(is.null(dZ)){ #sometimes user don't specify the Z matrices
dZ <- unlist(lapply(Z,function(x){dim(x$K)[1]}))
}
if(!is.null(X)){
dZ <- c(dZ,dim(X)[1])
}
dall <- unlist(c(dZ,dim(as.matrix(Y))[1]))
if(length(which(!duplicated(dall))) > 1){
if(is.null(X)){
stop("Matrices Y and Z's should have the same number of individuals. \nPlease check the dimensions of your matrices.", call. = FALSE)
}else{
stop("Matrices Y, X and Z's should have the same number of individuals. \nPlease check the dimensions of your matrices.", call. = FALSE)
}
}
#########*****************************
for(bb in 1:length(Z)){
ss1 <- colnames(Z[[bb]]$Z) == colnames(Z[[bb]]$K)
if(length(which(!ss1))>0){
print(paste("Names of columns in matrices Z and K for the",bb,"th random effect do not match."))
print("This can lead to incorrect estimation of variance components. Double check.")
}
}
#########*****************************
# if(!is.null(W)){
# cat("Response is imputed for estimation of variance components in GWAS models.\n")
# }
#########*****************************
if(!is.null(X)){
if(is.list(X)){
stop("Multivariate models only accept one incidence matrix for fixed effects (X). Please modifiy your X argument.",call. = FALSE)
}
}
#if(!silent){cat("Running multivariate model\n")}
if(check.model){
if(is.list(Z)){
if(is.list(Z[[1]])){ ### -- if is a 2 level list -- ##
provided <- lapply(Z, names)
for(s in 1:length(provided)){ #for each random effect =============================
provided2 <- names(Z[[s]])
if(length(provided2) ==1){ #----the 's' random effect has one matrix only----
if(provided2 == "K"){ #user only provided K
#zz <- diag(length(y))#model.matrix(~rownames(Z[[s]][[1]]))
zz <- diag(nrow(as.matrix(Y)))
#colnames(zz) <- rownames(Z[[s]][[1]])
Z[[s]] <- list(Z=zz, K=Z[[s]][[1]])
}
if(provided2 == "Z"){ # user only provided Z
#kk <- diag(dim(Z[[s]][[1]])[2])
kk <- diag(dim(Z[[s]][[1]])[2])
attributes(kk)$diagon <- TRUE
#rownames(kk) <- colnames(Z[[s]][[1]]); colnames(kk) <- rownames(kk)
Z[[s]] <- list(Z=Z[[s]][[1]],K=kk)
}
}else{ #----the 's' random effect has two matrices----
dido<-lapply(Z[[s]], dim) # dimensions of Z and K
condi<-(dido$Z[2] == dido$K[1] & dido$Z[2] == dido$K[2])
# condition, column size on Z matches with a square matrix K
if(!condi){
stop(paste("ERROR! In the",s,"th random effect you have provided or created an incidence \nmatrix with dimensions:",dido$Z[1],"rows and",dido$Z[2],"columns. Therefore the \nvariance-covariance matrix(K) for this random effect expected was a \nsquare matrix with dimensions",dido$Z[2],"x",dido$Z[2]),", but you provided a",dido$K[1],"x",dido$K[2]," matrix \nas a variance-covariance matrix. Please double check your matrices. Might be that you didn't drop all levels while subsetting your data, if so try: yourdata <- droplevels(yourdata).")
## stop()
}
}#---------------------------------------------------------------------------
} #for each random effect end =================================================
}else{ # if is a one-level list !!!!!!!!!!!!!
if(length(Z) == 1){ ## -- if the user only provided one matrix -- ##
provided <- names(Z)
if(provided == "K"){
#zz <- diag(length(y))
zz <- diag(nrow(as.matrix(Y)))
Z <- list(Z=zz, K=Z[[1]])
}
if(provided == "Z"){
#kk <- diag(dim(Z[[1]])[2])
kk <- diag(dim(Z[[1]])[2])
attributes(kk)$diagon <- TRUE
#rownames(kk) <- colnames(Z[[1]]); colnames(kk) <- rownames(kk)
Z <- list(Z=Z[[1]],K=kk)
}
}else{ # there's 2 matrices in Z
dido<-lapply(Z, dim) # dimensions of Z and K
condi<-(dido$Z[2] == dido$K[1] & dido$Z[2] == dido$K[2])
# condition, column size on Z matches with a square matrix K
if(!condi){
stop(paste("ERROR! In the",s,"th random effect you have provided or created an incidence \nmatrix with dimensions:",dido$Z[1],"rows and",dido$Z[2],"columns. Therefore the \nvariance-covariance matrix(K) for this random effect expected was a \nsquare matrix with dimensions",dido$Z[2],"x",dido$Z[2]),", but you provided a",dido$K[1],"x",dido$K[2]," matrix \nas a variance-covariance matrix. Please double check your matrices. Might be that you didn't drop all levels while subsetting your data, if so try: yourdata <- droplevels(yourdata).")
## stop()
}else{Z=list(Z=Z)}
}
}
}else{
if(is.null(Z)){ # the user is not using the random part
stop("Error. No random effects specified in the model. \nPlease use 'lm' or provide a diagonal matrix in Z\ni.e. Zu = list(A=list(Z=diag(length(y))))\n")
## stop()
}else{
#stop;
if(verbose) {
cat("\nThe parameter 'Z' needs to be provided in a 2-level list structure. \n\nPlease see help typing ?mmer and look at the 'Arguments' section\n")
cat("\nIf no random effects provided, the model will be fitted using the 'lm' function\n\n")
}
}
}
}
Y2 = matrix(Y)
rownames(Y2) = names(Y)
colnames(Y2) = "Yield"
Y = Y2
x = matrix(rep(1, nrow(Y))) # fixed effect (only intercept)
colnames(x) = "(Intercept)"
X = list(x)
Gx = list(matrix(1)) # list of constraints for fixed effects
ZETA = lapply(Z, function(x) return(x[[1]])) # list containing only incidence matrices of random effects (sparse form)
## names(ZETA) = NULL
for (i in 1:length(ZETA)){
attr(ZETA[[i]], "assign") = rep(1, ncol(ZETA[[i]]))
attr(ZETA[[i]], "contrasts") = list("contr.treatment")
names(attr(ZETA[[i]], "contrasts")) = paste0("QTL ", i)
}
ZETA <- lapply(ZETA,function(x){as(x, Class = "sparseMatrix")})
K = lapply(Z, function(x) return(x[[2]])) # list of varcov matrices of random effects (sparse form)
## names(K) = NULL
## for (i in 1:length(K)){
## dimnames(K[[i]]) = NULL
## }
R = list(diag(1, nrow = nrow(Y), ncol = nrow(Y))) # list of residual matrices (sparse form)
rownames(R[[1]]) = rownames(Y)
colnames(R[[1]]) = paste0("units",sort(as.character(seq(1, nrow(Y), 1)), decreasing = F))
attr(R[[1]], "assign") = rep(1, nrow(R[[1]]))
attr(R[[1]], "constrasts") = list(units = "contr.treatment")
R <- lapply(R,function(x){as(x, Class = "sparseMatrix")})
GES = rep(list(matrix(0.15)), length(ZETA)) # list of initial values for variance components
GES[[length(GES)+1]] = matrix(0.75)
GESI = rep(list(matrix(1)), length(ZETA)+1) # list of constraints applied to random effects
selected = FALSE
getPEV = TRUE
stepweight <- rep(0.9,iters); stepweight[1:2] <- c(0.5,0.7)
emupdate <- rep(0, iters)
ws = rep(1, nrow(Y))
## RES = MNR(Y, X,Gx,ZETA,K,R,GES,GESI, ws, iters, tolpar, tolparinv, selected,getPEV,verbose, FALSE, stepweight, emupdate)
RES = .Call("_qtlpoly_MNR", PACKAGE = "qtlpoly",Y, X,Gx,ZETA,K,R,GES,GESI, ws, iters, tolpar, tolparinv, selected,getPEV,verbose, FALSE, stepweight, emupdate)
RES$alleles = rownames(K[[1]])
## RES <- MNR(Y=Y,X=X,ZETA=Z,R=R,W=W,init=init,iters=iters,tolpar=tolpar,
## tolparinv = tolparinv,draw=draw,silent=silent,
## constraint = constraint,EIGEND = EIGEND,
## forced=forced,IMP=IMP,restrained=restrained, REML=REML,
## init.equal = init.equal)
## class(RES)<-c("MMERM")
## layout(matrix(1,1,1))
return(RES)
}
Try the qtlpoly package in your browser
Any scripts or data that you put into this service are public.
qtlpoly documentation built on May 29, 2024, 2:14 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions mmer_adapted summary.qtlpoly.fitted fit_model
Documented in fit_model mmer_adapted summary.qtlpoly.fitted
#' Fits multiple QTL models
#'
#' Fits alternative multiple QTL models by performing variance component estimation using REML.
#'
#' @param data an object of class \code{qtlpoly.data}.
#'
#' @param model an object of class \code{qtlpoly.profile} or \code{qtlpoly.remim}.
#'
#' @param probs a character string indicating if either \code{"joint"} (genotypes) or \code{"marginal"} (parental gametes) conditional probabilities should be used.
#'
#' @param polygenes a character string indicating if either \code{"none"}, \code{"most"} or \code{"all"} QTL should be used as polygenes.
#'
#' @param keep if \code{TRUE} (default), stores all matrices and estimates from fitted model; if \code{FALSE}, nothing is stored.
#'
#' @param verbose if \code{TRUE} (default), current progress is shown; if
#' \code{FALSE}, no output is produced.
#'
#' @param object an object of class \code{qtlpoly.fitted} to be summarized.
#'
#' @param pheno.col a numeric vector with the phenotype column numbers to be summarized; if \code{NULL} (default), all phenotypes from \code{'data'} will be included.
#'
#' @param ... currently ignored
#'
#' @return An object of class \code{qtlpoly.fitted} which contains a list of \code{results} for each trait with the following components:
#'
#' \item{pheno.col}{a phenotype column number.}
#' \item{fitted}{a \pkg{sommer} object of class \code{mmer}.}
#' \item{qtls}{a data frame with information from the mapped QTL.}
#'
#' @seealso \code{\link[qtlpoly]{read_data}}, \code{\link[qtlpoly]{remim}}
#'
#' @examples
#' \donttest{
#' # Estimate conditional probabilities using mappoly package
#' library(mappoly)
#' library(qtlpoly)
#' genoprob4x = lapply(maps4x[c(5)], calc_genoprob)
#' data = read_data(ploidy = 4, geno.prob = genoprob4x, pheno = pheno4x, step = 1)
#'
#' # Search for QTL
#' remim.mod = remim(data = data, pheno.col = 1, w.size = 15, sig.fwd = 0.0011493379,
#' sig.bwd = 0.0002284465, d.sint = 1.5, n.clusters = 1)
#'
#' # Fit model
#' fitted.mod = fit_model(data=data, model=remim.mod, probs="joint", polygenes="none")
#' }
#'
#' @author Guilherme da Silva Pereira, \email{gdasilv@@ncsu.edu}
#'
#' @references
#' Covarrubias-Pazaran G (2016) Genome-assisted prediction of quantitative traits using the R package sommer. \emph{PLoS ONE} 11 (6): 1–15. \doi{10.1371/journal.pone.0156744}.
#'
#' Pereira GS, Gemenet DC, Mollinari M, Olukolu BA, Wood JC, Mosquera V, Gruneberg WJ, Khan A, Buell CR, Yencho GC, Zeng ZB (2020) Multiple QTL mapping in autopolyploids: a random-effect model approach with application in a hexaploid sweetpotato full-sib population, \emph{Genetics} 215 (3): 579-595. \doi{10.1534/genetics.120.303080}.
#'
#' @useDynLib qtlpoly, .registration = TRUE
#' @importFrom Rcpp evalCpp
#' @export fit_model
fit_model <- function(data, model, probs="joint", polygenes="none", keep=TRUE, verbose=TRUE, pheno.col = NULL) {
results <- vector("list", length(model$results))
names(results) <- names(model$results)
if(probs=="marginal") {
ploidy <- data$ploidy
Palleles <- letters[1:ploidy]
Pgametes <- lapply(combn(Palleles, ploidy/2, simplify = FALSE), paste, collapse="")
Qalleles <- letters[(ploidy+1):(2*ploidy)]
Qgametes <- lapply(combn(Qalleles, ploidy/2, simplify = FALSE), paste, collapse="")
sibs <- sapply( Pgametes, FUN=function(x) paste(x, Pgametes, sep="") )
Pia <- matrix(data = NA, nrow = length(Pgametes), ncol = length(Pgametes))
for(i in 1:ncol(sibs)) {
for(j in 1:nrow(sibs)) {
Pia[i,j] <- (ploidy-length(unique(strsplit(sibs[i,j], "")[[1]])))/(ploidy/2)
}
}
Pib <- Pia # Pi1 = Pi2
colnames(Pia) <- rownames(Pia) <- unlist(Pgametes)
colnames(Pib) <- rownames(Pib) <- unlist(Qgametes)
indnames <- dimnames(data$Z)[[3]]
mrknames <- dimnames(data$Z)[[2]]
nmrk = length(mrknames)
nind = length(indnames)
Za <- array(data = NA, dim = c(length(Pgametes), nmrk, nind), dimnames = list(c(Pgametes), c(mrknames), c(indnames))) # Za = Z1
Zb <- array(data = NA, dim = c(length(Qgametes), nmrk, nind), dimnames = list(c(Qgametes), c(mrknames), c(indnames))) # Zb = Z2
for(m in 1:nmrk) {
for(i in 1:nind) {
Za[,m,i] <- colSums(matrix(data$Z[,m,i], nrow=length(Pgametes)))
Zb[,m,i] <- rowSums(matrix(data$Z[,m,i], nrow=length(Qgametes)))
}
}
}
for(p in 1:length(model$results)) {
pheno.col <- model$results[[p]]$pheno.col
qtl.mrk <- model$results[[p]]$qtl[,"Nmrk"]
qtl.lgr <- model$results[[p]]$qtl[,"LG"]
qtl.pos <- model$results[[p]]$qtl[,"Pos"]
if(verbose) {
if(length(qtl.mrk) == 1 && verbose) cat("There is ", length(qtl.mrk), " QTL in the model for trait ", pheno.col, " ", sQuote(colnames(data$pheno)[pheno.col]), ". Fitting model... ", sep="")
if(length(qtl.mrk) >= 2 && verbose) cat("There are ", length(qtl.mrk), " QTL in the model for trait ", pheno.col, " ", sQuote(colnames(data$pheno)[pheno.col]), ". Fitting model... ", sep="")
}
if(!is.null(model$results[[p]]$qtls)) {
nqtl <- dim(model$results[[p]]$qtls)[1]
ind <- rownames(data$pheno)[which(!is.na(data$pheno[,pheno.col]))]
Y <- data$pheno[ind,pheno.col]
Zstar <- diag(length(ind))
ETA <- NULL
if(nqtl > 1) {
if(probs=="joint" && polygenes=="none") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
for(q in 1:nqtl) {
eta <- list(list(Z=t(data$Z[,markers[q],ind]),K=data$Pi))
names(eta) <- paste("g", q, sep="")
ETA <- c(ETA, eta)
}
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
qtl <- c()
for(q in 1:nqtl) {
H2 <- sum(c(fitted$sigma)[q])/sum(c(fitted$sigma))
qtl <- c(qtl, c(model$results[[p]]$qtls[q,c(1:3)], NA, c(fitted$sigma)[q], NA, H2))
}
H2 <- sum(c(fitted$sigma)[1:nqtl])/sum(c(fitted$sigma))
qtl <- c(qtl, NA, NA, NA, c(fitted$Beta), NA, c(fitted$sigma)[nqtl+1], H2)
qtls <- as.data.frame(matrix(qtl, ncol=7, byrow=TRUE), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(g)", "Var(e)", "h2")
}
if(probs=="joint" && polygenes=="most") {
fitted <- vector("list", nqtl)
qtl <- c()
for(q in 1:nqtl) {
markers <- model$results[[p]]$qtls[-q,"Nmrk"]
Gstar <- apply(data$G[ind,ind,markers], MARGIN = c(1,2), sum)/length(markers); Gstar[1:5,1:5]
m <- model$results[[p]]$qtls[q,"Nmrk"]
ETA <- list(g=list(Z=t(data$Z[,m,ind]),K=data$Pi),gstar=list(Z=Zstar,K=Gstar))
fitted[[q]] <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- c(fitted[[q]]$sigma)[1]/sum(c(fitted[[q]]$sigma))
qtl <- c(qtl, c(model$results[[p]]$qtls[q,c(1:3)], c(fitted[[q]]$Beta), c(fitted[[q]]$sigma), H2))
}
qtls <- as.data.frame(matrix(qtl, ncol=8, byrow=TRUE), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(g)", "Var(g*)", "Var(e)", "h2")
}
if(probs=="joint" && polygenes=="all") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
Gstar <- apply(data$G[ind,ind,markers], MARGIN = c(1,2), sum)/length(markers); Gstar[1:5,1:5]
ETA <- list(gstar=list(Z=Zstar,K=Gstar))
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- c(fitted$sigma)[1]/sum(c(fitted$sigma))
qtl <- c(rep(NA, 4*nqtl), c(fitted$Beta), c(fitted$sigma), H2)
qtls <- as.data.frame(cbind(rbind(model$results[[p]]$qtls[,c(1:3)], rep(NA, 3)), matrix(qtl, ncol=4, byrow=TRUE), stringsAsFactors=FALSE))
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(g*)", "Var(e)", "h2")
}
if(probs=="marginal" && polygenes=="none") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
for(q in 1:nqtl) {
eta <- list(list(Z=t(Za[,markers[q],ind]),K=Pia),list(Z=t(Zb[,markers[q],ind]),K=Pib),list(Z=t(data$Z[,markers[q],ind]),K=data$Pi))
names(eta) <- paste(c("a","b","c"), rep(q,3), sep="")
ETA <- c(ETA, eta)
}
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
qtl <- c()
for(q in 1:nqtl) {
H2 <- sum(c(fitted$sigma)[(1:3)+((q-1)*3)])/sum(c(fitted$sigma))
qtl <- c(qtl, c(model$results[[p]]$qtls[q,c(1:3)], NA, c(fitted$sigma)[(1:3)+((q-1)*3)], NA, H2))
}
H2 <- sum(c(fitted$sigma)[1:(nqtl*3)])/sum(c(fitted$sigma))
qtl <- c(qtl, NA, NA, NA, c(fitted$Beta), rep(NA, 3), c(fitted$sigma)[(nqtl*3)+1], H2)
qtls <- as.data.frame(matrix(qtl, ncol=9, byrow=TRUE), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(a)", "Var(b)", "Var(d)", "Var(e)", "h2")
}
if(probs=="marginal" && polygenes=="most") {
fitted <- vector("list", nqtl)
qtl <- c()
for(q in 1:nqtl) {
markers <- model$results[[p]]$qtls[-q,"Nmrk"]
Gstar <- apply(data$G[ind,ind,markers], MARGIN = c(1,2), sum)/length(markers); Gstar[1:5,1:5]
m <- model$results[[p]]$qtls[q,"Nmrk"]
ETA <- list(a=list(Z=t(Za[,m,ind]),K=Pia),b=list(Z=t(Zb[,m,ind]),K=Pib),d=list(Z=t(data$Z[,m,ind]),K=data$Pi),gstar=list(Z=Zstar,K=Gstar))
fitted[[q]] <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- sum(c(fitted[[q]]$sigma)[(1:3)])/sum(c(fitted[[q]]$sigma))
qtl <- c(qtl, c(model$results[[p]]$qtls[q,c(1:3)], c(fitted[[q]]$Beta), c(fitted[[q]]$sigma), H2))
}
qtls <- as.data.frame(matrix(qtl, ncol=10, byrow=TRUE), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(a)", "Var(b)", "Var(d)", "Var(g*)", "Var(e)", "h2")
}
if(probs=="marginal" && polygenes=="all") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
A <- array(data = NA, dim = c(data$nind, data$nind, nqtl))
B <- array(data = NA, dim = c(data$nind, data$nind, nqtl))
for(q in 1:nqtl) {
A[,,q] <- t(Za[,markers[q],])%*%Pia%*%Za[,markers[q],]
B[,,q] <- t(Zb[,markers[q],])%*%Pib%*%Zb[,markers[q],]
}
Astar <- apply(A[ind,ind,], MARGIN = c(1,2), sum)/length(markers); Astar[1:5,1:5]
Bstar <- apply(B[ind,ind,], MARGIN = c(1,2), sum)/length(markers); Bstar[1:5,1:5]
Dstar <- apply(data$G[ind,ind,markers], MARGIN = c(1,2), sum)/length(markers); Dstar[1:5,1:5]
ETA <- list(astar=list(Z=Zstar,K=Astar),bstar=list(Z=Zstar,K=Bstar),dstar=list(Z=Zstar,K=Dstar))
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- sum(c(fitted$sigma)[1:3])/sum(c(fitted$sigma))
qtl <- c(rep(NA, 6*nqtl), c(fitted$Beta), c(fitted$sigma), H2)
qtls <- as.data.frame(cbind(rbind(model$results[[p]]$qtls[,c(1:3)], rep(NA, 3)), matrix(qtl, ncol=6, byrow=TRUE), stringsAsFactors=FALSE))
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(a*)", "Var(b*)", "Var(d*)", "Var(e)", "h2")
}
}
if(nqtl == 1) {
if(probs=="joint") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
ETA <- list(g=list(Z=t(data$Z[,markers,ind]),K=data$Pi))
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- sum(c(fitted$sigma)[1])/sum(c(fitted$sigma))
qtl <- c(c(fitted$Beta), c(fitted$sigma), H2)
qtls <- as.data.frame(c(model$results[[p]]$qtls[,c(1:3)], matrix(qtl, ncol=4, byrow=TRUE)), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(g)", "Var(e)", "h2")
}
if(probs=="marginal") {
markers <- model$results[[p]]$qtls[,"Nmrk"]
ETA <- list(a=list(Z=t(Za[,markers,ind]),K=Pia),b=list(Z=t(Zb[,markers,ind]),K=Pib),d=list(Z=t(data$Z[,markers,ind]),K=data$Pi))
fitted <- mmer_adapted(Y=Y, Z=ETA, silent=TRUE, date.warning=FALSE)
H2 <- sum(c(fitted$sigma)[1:3])/sum(c(fitted$sigma))
qtl <- c(c(fitted$Beta), c(fitted$sigma), H2)
qtls <- as.data.frame(c(model$results[[p]]$qtls[,c(1:3)], matrix(qtl, ncol=6, byrow=TRUE)), stringsAsFactors=FALSE)
colnames(qtls) <- c("LG", "Pos", "Nmrk", "Intercept", "Var(a)", "Var(b)", "Var(d)", "Var(e)", "h2")
}
}
if(verbose) cat("Done! \n\n", sep="")
}
if(is.null(model$results[[p]]$qtls)) {
if(verbose) cat("There are no QTL in the model for trait ", pheno.col, " ", sQuote(colnames(data$pheno)[pheno.col]), ". No model has been fitted! \n\n", sep="")
fitted <- NULL
qtls <- NULL
}
if(!keep) fitted <- NULL
results[[p]] <- list(
pheno.col=pheno.col,
fitted=fitted,
qtls=qtls)
}
structure(list(data=deparse(substitute(data)),
model=deparse(substitute(model)),
pheno.col=model$pheno.col,
probs=probs,
polygenes=polygenes,
results=results),
class="qtlpoly.fitted")
}
#' @rdname fit_model
#' @export
summary.qtlpoly.fitted <- function(object, pheno.col=NULL, ...) {
if(any(class(object) == "qtlpoly.fitted")) cat("This is an object of class 'qtlpoly.fitted'\n")
if(is.null(pheno.col)) {
pheno.col <- 1:length(object$results)
} else {
pheno.col <- which(object$pheno.col %in% pheno.col)
}
for(p in pheno.col) {
cat("\n* Trait", object$results[[p]]$pheno.col, sQuote(names(object$results)[[p]]), "\n")
if(!is.null(object$results[[p]]$qtls)) print(object$results[[p]]$qtls)
else cat("There are no QTL in the model \n")
}
}
#' Adapt to sommer MNR core function
#'
#' Adapts genomic incidence and relationship (varcov) matrices to run using sommer's C++ core function (v. 4.0 or higher)
#' Function adapted from sommer v. 3.6 (Author: Giovanny Covarrubias-Pazaran)
#'
#' @keywords internal
#'
#' @author Gabriel de Siqueira Gesteira, \email{gdesiqu@@ncsu.edu}
#'
#' @references
#' Covarrubias-Pazaran G (2016) Genome-assisted prediction of quantitative traits using the R package sommer. \emph{PLoS ONE} 11 (6): 1–15. \doi{10.1371/journal.pone.0156744}.
#'
#' Pereira GS, Gemenet DC, Mollinari M, Olukolu BA, Wood JC, Mosquera V, Gruneberg WJ, Khan A, Buell CR, Yencho GC, Zeng ZB (2020) Multiple QTL mapping in autopolyploids: a random-effect model approach with application in a hexaploid sweetpotato full-sib population, \emph{Genetics} 215 (3): 579-595. \doi{10.1534/genetics.120.303080}.
mmer_adapted <- function(Y,X=NULL,Z=NULL,R=NULL,W=NULL,method="NR",init=NULL,iters=20,tolpar=1e-3,
tolparinv=1e-6,draw=FALSE,silent=FALSE, constraint=TRUE,
EIGEND=FALSE, forced=NULL, IMP=FALSE, complete=TRUE,
check.model=TRUE, restrained=NULL, REML=TRUE, init.equal=TRUE,
date.warning=TRUE, verbose=FALSE){
#R=NULL
if (inherits(Y, "formula")){
stop("\nYou have to use mmer2 function for formula-based models.\n", call. = FALSE)
}
diso <- dim(as.data.frame(Y))[2]
## control for 2-level list structure
if(!is.list(Z)){
stop("Please provide the Z parameter as a 2 level list structure.\nFor example for 2 random effects 'A' and 'B' do:\n ETA <- list( A=list( Z=myZ1, K=myK1 ) , B=list( Z=myZ2, K=myK2 ) )\n mod <- mmer(Y=y, Z=ETA)\nwhere Z's and K's are the incidence and var-covar matrices respectively.\nIf any Z or K is not provided, an identity matrix will be assumed. ",call. = FALSE)
}else{
if(!is.list(Z[[1]])){
stop("Please provide the Z parameter as a 2 level list structure.\nFor example for 2 random effects 'A' and 'B' do:\n ETA <- list( A=list( Z=myZ1, K=myK1 ) , B=list( Z=myZ2, K=myK2 ) )\n mod <- mmer(Y=y, Z=ETA)\nwhere Z's and K's are the incidence and var-covar matrices respectively.\nIf any Z or K is not provided, an identity matrix will be assumed. ",call. = FALSE)
}
}
## control for Z-K names
zzzkkk <- unlist(lapply(Z,function(x){length(names(x))}))
badRE <- which(zzzkkk==0) # BAD RE WITH NO NAMES
if(length(badRE)>0){
stop("Please when specifying a random effect use the names; \n'Z' for incidence and 'K' for variance-covariance matrices.\nFor example for 1 random effect (i.e. named 'A') model do:\n ETA <- list( A=list( Z=M1, K=M2) )\n mod <- mmer(Y=y, Z=ETA)\nSpecifying at least one; Z or K. You need to specify if is a 'Z' or 'K' \nsince this is the only way the program distinguishes between matrices.",call. = FALSE)
}
#########*****************************
## make sure user don't provide the same names for random effects
his.names <- names(Z)
if(!is.null(his.names)){
badnames <- which(duplicated(his.names))
if(length(badnames)>0){
his.names[badnames] <- paste(his.names[badnames],1:(length(his.names[badnames])), sep=".")
names(Z) <- his.names
}
}
dZ <- unlist(lapply(Z,function(x){dim(x$Z)[1]}))
if(is.null(dZ)){ #sometimes user don't specify the Z matrices
dZ <- unlist(lapply(Z,function(x){dim(x$K)[1]}))
}
if(!is.null(X)){
dZ <- c(dZ,dim(X)[1])
}
dall <- unlist(c(dZ,dim(as.matrix(Y))[1]))
if(length(which(!duplicated(dall))) > 1){
if(is.null(X)){
stop("Matrices Y and Z's should have the same number of individuals. \nPlease check the dimensions of your matrices.", call. = FALSE)
}else{
stop("Matrices Y, X and Z's should have the same number of individuals. \nPlease check the dimensions of your matrices.", call. = FALSE)
}
}
#########*****************************
for(bb in 1:length(Z)){
ss1 <- colnames(Z[[bb]]$Z) == colnames(Z[[bb]]$K)
if(length(which(!ss1))>0){
print(paste("Names of columns in matrices Z and K for the",bb,"th random effect do not match."))
print("This can lead to incorrect estimation of variance components. Double check.")
}
}
#########*****************************
# if(!is.null(W)){
# cat("Response is imputed for estimation of variance components in GWAS models.\n")
# }
#########*****************************
if(!is.null(X)){
if(is.list(X)){
stop("Multivariate models only accept one incidence matrix for fixed effects (X). Please modifiy your X argument.",call. = FALSE)
}
}
#if(!silent){cat("Running multivariate model\n")}
if(check.model){
if(is.list(Z)){
if(is.list(Z[[1]])){ ### -- if is a 2 level list -- ##
provided <- lapply(Z, names)
for(s in 1:length(provided)){ #for each random effect =============================
provided2 <- names(Z[[s]])
if(length(provided2) ==1){ #----the 's' random effect has one matrix only----
if(provided2 == "K"){ #user only provided K
#zz <- diag(length(y))#model.matrix(~rownames(Z[[s]][[1]]))
zz <- diag(nrow(as.matrix(Y)))
#colnames(zz) <- rownames(Z[[s]][[1]])
Z[[s]] <- list(Z=zz, K=Z[[s]][[1]])
}
if(provided2 == "Z"){ # user only provided Z
#kk <- diag(dim(Z[[s]][[1]])[2])
kk <- diag(dim(Z[[s]][[1]])[2])
attributes(kk)$diagon <- TRUE
#rownames(kk) <- colnames(Z[[s]][[1]]); colnames(kk) <- rownames(kk)
Z[[s]] <- list(Z=Z[[s]][[1]],K=kk)
}
}else{ #----the 's' random effect has two matrices----
dido<-lapply(Z[[s]], dim) # dimensions of Z and K
condi<-(dido$Z[2] == dido$K[1] & dido$Z[2] == dido$K[2])
# condition, column size on Z matches with a square matrix K
if(!condi){
stop(paste("ERROR! In the",s,"th random effect you have provided or created an incidence \nmatrix with dimensions:",dido$Z[1],"rows and",dido$Z[2],"columns. Therefore the \nvariance-covariance matrix(K) for this random effect expected was a \nsquare matrix with dimensions",dido$Z[2],"x",dido$Z[2]),", but you provided a",dido$K[1],"x",dido$K[2]," matrix \nas a variance-covariance matrix. Please double check your matrices. Might be that you didn't drop all levels while subsetting your data, if so try: yourdata <- droplevels(yourdata).")
## stop()
}
}#---------------------------------------------------------------------------
} #for each random effect end =================================================
}else{ # if is a one-level list !!!!!!!!!!!!!
if(length(Z) == 1){ ## -- if the user only provided one matrix -- ##
provided <- names(Z)
if(provided == "K"){
#zz <- diag(length(y))
zz <- diag(nrow(as.matrix(Y)))
Z <- list(Z=zz, K=Z[[1]])
}
if(provided == "Z"){
#kk <- diag(dim(Z[[1]])[2])
kk <- diag(dim(Z[[1]])[2])
attributes(kk)$diagon <- TRUE
#rownames(kk) <- colnames(Z[[1]]); colnames(kk) <- rownames(kk)
Z <- list(Z=Z[[1]],K=kk)
}
}else{ # there's 2 matrices in Z
dido<-lapply(Z, dim) # dimensions of Z and K
condi<-(dido$Z[2] == dido$K[1] & dido$Z[2] == dido$K[2])
# condition, column size on Z matches with a square matrix K
if(!condi){
stop(paste("ERROR! In the",s,"th random effect you have provided or created an incidence \nmatrix with dimensions:",dido$Z[1],"rows and",dido$Z[2],"columns. Therefore the \nvariance-covariance matrix(K) for this random effect expected was a \nsquare matrix with dimensions",dido$Z[2],"x",dido$Z[2]),", but you provided a",dido$K[1],"x",dido$K[2]," matrix \nas a variance-covariance matrix. Please double check your matrices. Might be that you didn't drop all levels while subsetting your data, if so try: yourdata <- droplevels(yourdata).")
## stop()
}else{Z=list(Z=Z)}
}
}
}else{
if(is.null(Z)){ # the user is not using the random part
stop("Error. No random effects specified in the model. \nPlease use 'lm' or provide a diagonal matrix in Z\ni.e. Zu = list(A=list(Z=diag(length(y))))\n")
## stop()
}else{
#stop;
if(verbose) {
cat("\nThe parameter 'Z' needs to be provided in a 2-level list structure. \n\nPlease see help typing ?mmer and look at the 'Arguments' section\n")
cat("\nIf no random effects provided, the model will be fitted using the 'lm' function\n\n")
}
}
}
}
Y2 = matrix(Y)
rownames(Y2) = names(Y)
colnames(Y2) = "Yield"
Y = Y2
x = matrix(rep(1, nrow(Y))) # fixed effect (only intercept)
colnames(x) = "(Intercept)"
X = list(x)
Gx = list(matrix(1)) # list of constraints for fixed effects
ZETA = lapply(Z, function(x) return(x[[1]])) # list containing only incidence matrices of random effects (sparse form)
## names(ZETA) = NULL
for (i in 1:length(ZETA)){
attr(ZETA[[i]], "assign") = rep(1, ncol(ZETA[[i]]))
attr(ZETA[[i]], "contrasts") = list("contr.treatment")
names(attr(ZETA[[i]], "contrasts")) = paste0("QTL ", i)
}
ZETA <- lapply(ZETA,function(x){as(x, Class = "sparseMatrix")})
K = lapply(Z, function(x) return(x[[2]])) # list of varcov matrices of random effects (sparse form)
## names(K) = NULL
## for (i in 1:length(K)){
## dimnames(K[[i]]) = NULL
## }
R = list(diag(1, nrow = nrow(Y), ncol = nrow(Y))) # list of residual matrices (sparse form)
rownames(R[[1]]) = rownames(Y)
colnames(R[[1]]) = paste0("units",sort(as.character(seq(1, nrow(Y), 1)), decreasing = F))
attr(R[[1]], "assign") = rep(1, nrow(R[[1]]))
attr(R[[1]], "constrasts") = list(units = "contr.treatment")
R <- lapply(R,function(x){as(x, Class = "sparseMatrix")})
GES = rep(list(matrix(0.15)), length(ZETA)) # list of initial values for variance components
GES[[length(GES)+1]] = matrix(0.75)
GESI = rep(list(matrix(1)), length(ZETA)+1) # list of constraints applied to random effects
selected = FALSE
getPEV = TRUE
stepweight <- rep(0.9,iters); stepweight[1:2] <- c(0.5,0.7)
emupdate <- rep(0, iters)
ws = rep(1, nrow(Y))
## RES = MNR(Y, X,Gx,ZETA,K,R,GES,GESI, ws, iters, tolpar, tolparinv, selected,getPEV,verbose, FALSE, stepweight, emupdate)
RES = .Call("_qtlpoly_MNR", PACKAGE = "qtlpoly",Y, X,Gx,ZETA,K,R,GES,GESI, ws, iters, tolpar, tolparinv, selected,getPEV,verbose, FALSE, stepweight, emupdate)
RES$alleles = rownames(K[[1]])
## RES <- MNR(Y=Y,X=X,ZETA=Z,R=R,W=W,init=init,iters=iters,tolpar=tolpar,
## tolparinv = tolparinv,draw=draw,silent=silent,
## constraint = constraint,EIGEND = EIGEND,
## forced=forced,IMP=IMP,restrained=restrained, REML=REML,
## init.equal = init.equal)
## class(RES)<-c("MMERM")
## layout(matrix(1,1,1))
return(RES)
}
Try the qtlpoly package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.