v1.08

jendelman · May 27, 2024 · 032e795 · 032e795
1 parent 088069a
commit 032e795
Show file tree

Hide file tree

Showing 8 changed files with 103 additions and 61 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: StageWise
 Title: Two-stage analysis of multi-environment trials for genomic selection and GWAS
-Version: 1.07
+Version: 1.08
 Author: Jeffrey B. Endelman
 Maintainer: Jeffrey Endelman <[email protected]>
 Description: Fully efficient, two stage analysis of multi-environment trials, including directional dominance and multi-trait genomic selection

diff --git a/NEWS b/NEWS
@@ -1,3 +1,6 @@
+Changes in 1.08
+* summary.var and dominance functions now work for multi-trait
+
 Changes in 1.07
 * coerce_dpo updated to account for changes to the Matrix package
 

diff --git a/R/Stage2.R b/R/Stage2.R
@@ -246,25 +246,11 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL,
       tmp <- gsub("+",":loc+",paste(c("env",fix.eff.marker,dom),collapse="+"),fixed=T)
       model <- sub("FIX",paste(tmp,"loc",sep=":"),model,fixed=T)
     }
-    #if (n.mark > 0 | (non.add=="dom") | !is.null(fix.eff)) {
-    #if (n.loc==1) {
-    #  model <- sub("FIX",paste(c("env",fix.eff,fix.eff.marker,dom),collapse="+"),model,fixed=T)
-    #} else {
-    #  model <- sub("FIX",paste(c("env",paste0(c(fix.eff.marker,dom),":loc")),collapse="+"),model,fixed=T)
-    #}
-    #} else {
-    #  model <- sub("FIX","env",model,fixed=T)
-    #}
   } else {
     model <- "asreml(data=data,fixed=BLUE~FIX-1,random=~RANDOM,residual=~id(env.id):us(Trait)"
     tmp <- gsub("+",":Trait+",paste(c("env",fix.eff.marker,dom),collapse="+"),fixed=T)
     model <- sub("FIX",paste(tmp,"Trait",sep=":"),model,fixed=T)
 
-    # if (n.mark > 0 | (non.add=="dom")) {
-    #   model <- sub("FIX",paste(paste0(c("env",fix.eff.marker,dom),":Trait"),collapse="+"),model,fixed=T)
-    # } else {
-    #   model <- sub("FIX","env:Trait",model,fixed=T)
-    # }
   }
 
   if (is.null(geno)) {

diff --git a/R/dominance.R b/R/dominance.R
@@ -2,49 +2,74 @@
 #'
 #' Report dominance parameters
 #' 
-#' The dominance variance (Vd) and baseline heterosis (b) are quantified relative to additive variance (Va) and std. dev. (SDa), respectively. The estimate and SE of the ratios are calculated based on the delta method (Rice 2006, p. 162-166). Only implemented for single trait/location models.
+#' The dominance variance (Vd) and baseline heterosis (b) are quantified relative to additive variance (Va) and std. dev. (SDa), respectively. For single traits, the estimate and SE of the ratios are calculated based on the delta method (Rice 2007, p. 162-166). For a multi-trait model, the result is just the ratio of the estimates, with \code{index.coeff} specifying the coefficients of the standardized true values (see also \code{\link{blup}}) and \code{gamma} specifying the relative weight for non-additive to additive genetic merit (see also \code{\link{gain}}). 
 #' 
 #' @param params list returned by \code{\link{Stage2}}
 #' @param digits number of digits for rounding
+#' @param index.coeff merit index coefficients
+#' @param gamma contribution of non-additive values for genetic merit
 #'
 #' @return data frame with estimates and SE
 #'
 #' @references Rice JA (2007) Mathematical Statistics and Data Analysis, 3rd ed. Duxbury, Pacific Grove.
 #' @export
 
-dominance <- function(params, digits=2) {
+dominance <- function(params, digits=2, index.coeff=NULL, gamma=0) {
   if (!is.element("heterosis",names(params)))
     stop("Dominance model was not used")
-  if (nrow(params$heterosis) > 1)
-    stop("Unavailable for multi-location/trait models")
-
+
   id <- grep("dominance",params$vc$name,fixed=T)
   ia <- grep("additive",params$vc$name,fixed=T)
-  Vd.mean <- params$vc$estimate[id]
-  Vd.var <- params$vc$SE[id]^2
-  Va.mean <- params$vc$estimate[ia]
-  Va.var <- params$vc$SE[ia]^2
-  SDa.mean <- sqrt(Va.mean)*(1-2*Va.var/Va.mean^2)
-  SDa.var <- Va.var/(4*Va.mean)
-  b.mean <- params$heterosis$estimate
-  b.var <- params$heterosis$SE^2
-
-  ratio1.mean <- Vd.mean/Va.mean*(1+Va.var/Va.mean^2)
-  ratio1.SE <- sqrt(Va.var*Vd.mean^2/Va.mean^2+Vd.var)/Va.mean
-  ratio2.mean <- b.mean/SDa.mean*(1+SDa.var/SDa.mean^2)
-  ratio2.SE <- sqrt(SDa.var*b.mean^2/SDa.mean^2+b.var)/SDa.mean
+
+  n.trait <- nrow(params$heterosis)
+  if (n.trait > 1) {
+    stopifnot(!is.null(index.coeff))
+    traits <- params$heterosis$trait
+    Vd <- Va <- matrix(0,nrow=n.trait,ncol=n.trait,dimnames=list(traits,traits))
+
+    tmp <- strsplit(params$vc$name[ia],split=":")  
+    ix <- cbind(sapply(tmp,"[[",2),sapply(tmp,"[[",3))
+    Va[ix] <- params$vc$estimate[ia]
+    Va[upper.tri(Va)] <- Va[lower.tri(Va)]
+
+    tmp <- strsplit(params$vc$name[id],split=":")  
+    ix <- cbind(sapply(tmp,"[[",2),sapply(tmp,"[[",3))
+    Vd[ix] <- params$vc$estimate[id]
+    Vd[upper.tri(Vd)] <- Vd[lower.tri(Vd)]
+
+    stopifnot(sort(names(index.coeff))==sort(traits))
+    index.coeff <- index.coeff[match(traits,names(index.coeff))]
+    trait.scale <- sqrt(diag(Va)+gamma^2*diag(Vd))
+    index.coeff <- index.coeff/trait.scale
+
+    Va.mean <- as.numeric(crossprod(index.coeff,Va%*%index.coeff))
+    Vd.mean <- as.numeric(crossprod(index.coeff,Vd%*%index.coeff))
+    b.mean <- sum(params$heterosis$estimate*index.coeff)
+
+    ans <- data.frame(ratio=c("Vd/Va","b/SDa"),
+                      estimate=round(c(Vd.mean/Va.mean,b.mean/sqrt(Va.mean)),digits))
+    return(ans)
+
+  } else {
 
-  # if (length(ia) > 1) {
-  #  ans <- list('Vd/Va'=data.frame(trait=params$heterosis$trait,
-  #                           estimate=round(ratio1.mean,digits),
-  #                           SE=round(ratio1.SE,digits)),
-  #              'b/SDa'=data.frame(trait=params$heterosis$trait,
-  #                            estimate=round(ratio2.mean,digits),
-  #                            SE=round(ratio2.SE,digits)))
-  # } else {
-  ans <- data.frame(ratio=c("Vd/Va","b/SDa"),
-                    estimate=round(c(ratio1.mean,ratio2.mean),digits),
-                    SE=round(c(ratio1.SE,ratio2.SE),digits))
+    Vd.mean <- params$vc$estimate[id]
+    Vd.var <- params$vc$SE[id]^2
+    Va.mean <- params$vc$estimate[ia]
+    Va.var <- params$vc$SE[ia]^2
+    SDa.mean <- sqrt(Va.mean)*(1-2*Va.var/Va.mean^2)
+    SDa.var <- Va.var/(4*Va.mean)
+    b.mean <- params$heterosis$estimate
+    b.var <- params$heterosis$SE^2
+
+    ratio1.mean <- Vd.mean/Va.mean*(1+Va.var/Va.mean^2)
+    ratio1.SE <- sqrt(Va.var*Vd.mean^2/Va.mean^2+Vd.var)/Va.mean
+    ratio2.mean <- b.mean/SDa.mean*(1+SDa.var/SDa.mean^2)
+    ratio2.SE <- sqrt(SDa.var*b.mean^2/SDa.mean^2+b.var)/SDa.mean
+
+    ans <- data.frame(ratio=c("Vd/Va","b/SDa"),
+                      estimate=round(c(ratio1.mean,ratio2.mean),digits),
+                      SE=round(c(ratio1.SE,ratio2.SE),digits))
 
-  return(ans)
+    return(ans)
+  }
 }
diff --git a/R/summary.var.R b/R/summary.var.R
@@ -2,10 +2,14 @@
 #'
 #' Displays variances and correlations
 #' 
-#' For a single trait, the 'var' output is a data frame with two columns of information for the various effects: the first is the variance and the second is the proportion of variance explained (PVE), excluding the environment effect. For multiple locations or traits, the 'cor' output is the correlation matrix for additive effects (does not include fixed effect markers). For multiple traits, the variance and PVE results are returned as separate data frames.
+#' For a single trait, the 'var' output is a data frame with two columns of information for the various effects: the first is the variance and the second is the proportion of variance explained (PVE), excluding the environment effect. For multiple locations or traits, the 'cor' output is the correlation matrix for additive effects (does not include fixed effect markers). For multiple traits, the variance and PVE results are returned as separate data frames, unless \code{index.coeff} is used to create an index. 
+#' 
+#' The \code{index.coeff} are the coefficients of the standardized true values (see also \code{\link{blup}}). The argument \code{gamma} is the relative weight for non-additive to additive genetic merit (see also \code{\link{gain}}). 
 #' 
 #' @param object object of \code{\link{class_var}}
 #' @param digits number of digits for rounding
+#' @param index.coeff merit index coefficients
+#' @param gamma contribution of non-additive values for genetic merit
 #'
 #' @return List output that varies depending on the situation (see Details)
 #'
@@ -16,42 +20,56 @@ NULL
 
 setGeneric("summary")
 setMethod("summary",c(object="class_var"),
-          definition=function(object,digits=3){
+          definition=function(object,digits=3,index.coeff=NULL,gamma=0){
 
             n.trait <- ncol(object@resid)
+            if (n.trait > 1 & !is.null(index.coeff)) {
+              traits <- rownames(object@geno1)
+              stopifnot(sort(names(index.coeff))==sort(traits))
+              index.coeff <- index.coeff[match(traits,names(index.coeff))]
+              trait.scale <- sqrt(diag(object@geno1)+gamma^2*diag(object@geno2))
+              index.coeff <- index.coeff/trait.scale
+            }
             if (dim(object@vars)[1]==0) {
               pvar <- FALSE
             } else {
               pvar <- TRUE
             }
-            if (pvar)
-              vars <- apply(object@vars,3,diag)
-
-            if (n.trait > 1) {
+            if (pvar) {
+              if (is.null(index.coeff)) {
+                vars <- apply(object@vars,3,diag)
+              } else {
+                vars <- apply(object@vars,3,function(V){
+                  as.numeric(crossprod(index.coeff,V%*%index.coeff))
+                })
+              }
+            }
+
+            if (n.trait > 1 & is.null(index.coeff)) {
               cor.mat <- cov_to_cor(object@geno1)
               if (pvar) {
                 vars <- t(vars)
                 prop.var <- t(t(vars[-1,])/apply(vars[-1,],2,sum,na.rm=T))
                 prop.var <- round(prop.var[!is.na(prop.var[,1]),],3)
-
+                
                 vars <- apply(vars[!is.na(vars[,1]),],2,sigdig,digits=digits)
                 return(list(var=data.frame(vars),
-                          PVE=data.frame(prop.var),
-                          cor.mat=round(cor.mat,3)))
+                            PVE=data.frame(prop.var),
+                            cor.mat=round(cor.mat,3)))
               } else {
                 return(round(cor.mat,3))
               }
             } else {
-              if (nrow(object@geno1) > 1) {
+              if (nrow(object@geno1) > 1 & n.trait==1) {
                 #multiple locations
                 cor.mat <- cov_to_cor(object@geno1)
                 x <- data.frame(Variance=sigdig(vars,digits),
-                           PVE=round(c(NA,vars[-1]/sum(vars[-1],na.rm=T)),3))
+                                PVE=round(c(NA,vars[-1]/sum(vars[-1],na.rm=T)),3))
                 return(list(var=x[!is.na(x[,1]),],cor=round(cor.mat,3)))
               } else {
                 x <- data.frame(Variance=sigdig(vars,digits),
-                           PVE=round(c(NA,NA,vars[-(1:2)]/sum(vars[-(1:2)],na.rm=T)),3))
+                                PVE=round(c(NA,NA,vars[-(1:2)]/sum(vars[-(1:2)],na.rm=T)),3))
                 return(x[!is.na(x[,1]),])
-                }
+              }
             }
           })
diff --git a/docs/manual.pdf b/docs/manual.pdf
diff --git a/man/dominance.Rd b/man/dominance.Rd
diff --git a/man/summary.var.Rd b/man/summary.var.Rd