v1.11

jendelman · Jul 24, 2024 · ca61ea1 · ca61ea1
1 parent 3c1e618
commit ca61ea1
Show file tree

Hide file tree

Showing 12 changed files with 111 additions and 110 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.10
+Version: 1.11
 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.11
+* dominance function is now scaled to population, not just var component
+
 Changes in 1.10
 * Added support for multi-population/multi-ploidy analysis
 * Fixed error in Stage1 related to multiple experiments per env

diff --git a/R/Stage2.R b/R/Stage2.R
@@ -118,8 +118,8 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL,
                   vars=new(Class="class_var",
                            geno1=geno1, geno2=geno2, resid=resid.vc, B=as.matrix(B),
                            model=result[[1]]$vars@model,
-                           diagG=mean(sapply(result,function(x){x$vars@diagG})),
-                           diagD=mean(sapply(result,function(x){x$vars@diagD})),
+                           #diagG=mean(sapply(result,function(x){x$vars@diagG})),
+                           #diagD=mean(sapply(result,function(x){x$vars@diagD})),
                            vars=array(NA,dim=c(0,0,0)),
                            fix.eff.marker=result[[1]]$vars@fix.eff.marker)))
     }  
@@ -134,7 +134,7 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL,
     data <- data[!(data$env.id %in% data$env.id[missing]),]
   }
 
-  diagG <- diagD <- numeric(0)
+  #diagG <- diagD <- numeric(0)
   dom <- NULL
   if (!is.null(geno)) {
     stopifnot(is(geno,"class_geno"))
@@ -147,7 +147,7 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL,
     id.weights <- table(factor(data$id,levels=id))
     .GlobalEnv$asremlG <- geno@G[id,id]
     meanG <- as.numeric(pvar(V=.GlobalEnv$asremlG,weights=id.weights))
-    diagG <- mean(diag(.GlobalEnv$asremlG))
+    #diagG <- mean(diag(.GlobalEnv$asremlG))
 
     if (non.add=="dom") {
       .GlobalEnv$asremlD <- geno@D[id,id]
@@ -158,7 +158,7 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL,
       names(dom.covariate) <- id
       data$dom <- dom.covariate[data$id]
 
-      diagD <- mean(diag(.GlobalEnv$asremlD))
+      #diagD <- mean(diag(.GlobalEnv$asremlD))
       dom <- "dom"
     } 
   } else {
@@ -696,8 +696,8 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL,
   }
 
   out$vars <- new(Class="class_var",geno1=geno1.vc,geno2=geno2.vc,model=model,
-                  resid=resid.vc,diagG=diagG,diagD=diagD,
-                  vars=vars,B=B,fix.eff.marker=fix.eff.marker)
+                  resid=resid.vc,vars=vars,
+                  B=B,fix.eff.marker=fix.eff.marker)
 
   #random effects
   if (n.trait==1) {
@@ -783,12 +783,18 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL,
     }
   }
 
-  if (!is.null(geno))
+  if (!is.null(geno)) {
+    out$params$scale <- list(add=as.numeric(pvar(V=.GlobalEnv$asremlG)))
     rm("asremlG",envir = .GlobalEnv)
+  }
   if (!is.null(vcov))
     rm("asremlOmega",envir = .GlobalEnv)
-  if (non.add=="dom")
+  if (non.add=="dom") {
+    out$params$scale <- c(out$params$scale,
+                          list(dom=as.numeric(pvar(V=.GlobalEnv$asremlD)),
+                             heterosis=as.numeric(pvar(mu=dom.covariate))))
     rm("asremlD",envir = .GlobalEnv)
+  }
 
   return(out)
 }
diff --git a/R/blup_prep.R b/R/blup_prep.R
@@ -32,6 +32,8 @@ blup_prep <- function(data,vcov=NULL,geno=NULL,vars,mask=NULL,method=NULL) {
     stopifnot(method %in% c("MME","VINV"))
   }
 
+  if (vars@model > 0)
+    stopifnot(inherits(geno,"class_geno"))
   if (vars@model==3L)
     stopifnot(is(geno,"class_genoD"))
 
@@ -43,9 +45,9 @@ blup_prep <- function(data,vcov=NULL,geno=NULL,vars,mask=NULL,method=NULL) {
     stopifnot(n.trait > 1)
   } 
 
-  if (length(vars@diagG)>0) {
-    stopifnot(!is.null(geno))
-  }
+  #if (length(vars@diagG)>0) {
+  #  stopifnot(!is.null(geno))
+  #}
 
   data <- data[!is.na(data$BLUE),]
 
@@ -65,8 +67,9 @@ blup_prep <- function(data,vcov=NULL,geno=NULL,vars,mask=NULL,method=NULL) {
   }
 
   if (!is.null(geno)) {
-    stopifnot(is(geno,"class_geno"))
-    stopifnot(length(vars@diagG)>0)
+    stopifnot(inherits(geno,"class_geno"))
+    #stopifnot(length(vars@diagG)>0)
+    stopifnot(vars@model > 0)
     id <- intersect(data$id,rownames(geno@G))
     data <- data[data$id %in% id,]
     id <- rownames(geno@G)

diff --git a/R/class_var.R b/R/class_var.R
@@ -4,8 +4,6 @@
 #' @slot geno2 second genetic effect
 #' @slot model 0=no markers, 1=add, 2=add+g.resid, 3=add+dom
 #' @slot resid residual
-#' @slot diagG average diagonal element of the G matrix
-#' @slot diagD average diagonal element of the D matrix
 #' @slot vars variances for reporting
 #' @slot B var-cov matrix of fixed effects for gain
 #' @slot fix.eff.marker names of fixed effect markers
@@ -17,8 +15,6 @@ class_var <- setClass("class_var",slots=c(geno1="Matrix",
                                           geno2="Matrix",
                                           model="integer",
                                           resid="Matrix",
-                                          diagG="numeric",
-                                          diagD="numeric",
                                           vars="array",
                                           B="matrix",
                                           fix.eff.marker="character"))
diff --git a/R/dominance.R b/R/dominance.R
@@ -2,22 +2,22 @@
 #'
 #' Report dominance parameters
 #' 
-#' 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/loc 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}}). 
+#' The dominance variance (Vd) and baseline heterosis (b) are quantified relative to additive variance (Va) and std. dev. (SDa), respectively. As of v1.11, the variances are scaled to the population (previously, it was just the variance components). For a multi-trait/loc model, \code{index.coeff} specifies the coefficients of the standardized true values (see also \code{\link{blup}}), with \code{gamma} indicating the relative weight of non-additive to additive genetic merit for the standardization (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
+#' @return data frame with estimates 
 #'
-#' @references Rice JA (2007) Mathematical Statistics and Data Analysis, 3rd ed. Duxbury, Pacific Grove.
 #' @export
 
-dominance <- function(params, digits=2, index.coeff=NULL, gamma=0) {
+dominance <- function(params, index.coeff=NULL, gamma=0) {
   if (!is.element("heterosis",names(params)))
     stop("Dominance model was not used")
 
+  stopifnot("scale" %in% names(params))
+
   id <- grep("dominance",params$vc$name,fixed=T)
   ia <- grep("additive",params$vc$name,fixed=T)
   vc <- params$vc[,-1]
@@ -28,7 +28,6 @@ dominance <- function(params, digits=2, index.coeff=NULL, gamma=0) {
     stopifnot(!is.null(index.coeff))
     if (colnames(params$heterosis)[1]=="trait") {
       traits <- params$heterosis$trait
-      #Vd <- Va <- matrix(0,nrow=nhet,ncol=nhet,dimnames=list(traits,traits))
       Va <- f.cov.trait(vc[ia,],traits,us=(nhet>2))
       Vd <- f.cov.trait(vc[id,],traits,us=(nhet>2))
     } else {
@@ -49,31 +48,16 @@ dominance <- function(params, digits=2, index.coeff=NULL, gamma=0) {
     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 {
-
     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)
   }
-}
+
+  Vd <- Vd.mean*params$scale$dom + b.mean^2*params$scale$heterosis
+  Va <- Va.mean*params$scale$add
+
+  ans <- data.frame(ratio=c("Vd/Va","b/SDa"),
+                      estimate=c(Vd/Va, b.mean/sqrt(Va)))
+  return(ans)
+} 
diff --git a/docs/Vignette1.html b/docs/Vignette1.html