gwaR: Functions for performing GWA from GBLUP

#' A print method for gblup objects
#' @title print \code{\link{gblup}} object
#' @param gblup an object of the class gblup generated by \code{\link{gblup}}
#' @return Null. This function only prints to the screen:
#' @seealso \code{\link{gblup}} 
#' @export
#'  
print.gblup <- function(gblup) {
    nms <- all.vars(gblup$model$Vformula)
    fe <- gblup$coefm[!rownames(gblup$coef) %in% nms, , drop = F]
    tv <- fe[, 1]/fe[, 2]
    pv <- 2 * (1 - pnorm(abs(tv)))
    fe <- data.frame(fe, test.st = tv, p.value = pv)
    re <- gblup$coefm[rownames(gblup$coef) %in% nms, ]
    h2 <- re[, 1]/sum(re[, 1])
    re <- data.frame(re, prop.var = h2)
    cat("gblup analysis of trait:", gblup$name, "\n\n")
    
    cat("fixed effects equation:\n")
    print(gblup$model$formula, showEnv = F)
    cat("\nrandom effects equation:\n")
    print(gblup$model$Vformula, showEnv = F)
    
    cat("\nlog-likelihood:", gblup$llik, "converged in:", gblup$cycle, "iterations \n")
    
    cat("\nestimated fixed effects:\n")
    print(fe)
    cat("\nestimated variance components:\n")
    print(re)
}



#' A summary method for gblup objects
#' @title summary for \code{\link{gblup}} object
#' @param gblup an object of the class gblup generated by \code{\link{gblup}}
#' @param ehat a logical value to indicate if residuals should be reported
#' @param sigma a logical value to indicate if variance component matrix should be included
#' @param fe a logical value to indicate if a fixed effects table should be included
#' @return an object of the class summary.gblup with:
#' \itemize{ 
#'    \item {\code{trait}} {name of the trait} 
#'    \item {\code{vformula and formula}} {equations for random and fixed effects} 
#'    \item {\code{llik}} {log-likelihood at convergence} 
#'    \item {\code{cycle}} {number of iterations} 
#'    \item {\code{uhat}} {blup of individuals} 
#'    \item{\code{sigma}} {estimated variance components and h2}
#'    \item{\code{coefm}} {estimated fixed effects, s.e. and test statistics}
#'}
#' @seealso \code{\link{gblup}} 
#' @export
#'  
summary.gblup <- function(gblup, ehat = FALSE, sigma = FALSE, fe = FALSE) {
    # add trait name!!!
    
    uhat <- gblup$sigma[["G"]] * gblup$model$G %*% (gblup$Vinv %*% gblup$ehat)
    
    result <- list(uhat = uhat)
    result$name <- gblup$name
    if (sigma) 
        result$sigma <- gblup$sigma
    if (ehat) 
        result$ehat <- gblup$ehat
    if (fe) 
        result$coefm <- gblup$coefm
    result$llik <- gblup$llik
    result$cycle <- gblup$cycle
    result$formula <- gblup$model$formula
    result$Vformula <- gblup$model$Vformula
    class(result) <- "summary.gblup"
    return(result)
}

#' A print method for summary.gblup objects
#' @title print \code{\link{summary.gblup}} object
#' @param gblup an object of the class symmary.gblup generated by \code{\link{summary.gblup}}
#' @return Null. This function only prints to the screen:
#' @seealso \code{\link{summary.gblup}} 
#' @export
print.summary.gblup <- function(gblup) {
    nms <- all.vars(gblup$Vformula)
    fe <- gblup$coefm[!rownames(gblup$coefm) %in% nms, ]
    tv <- fe[, 1]/fe[, 2]
    pv <- 2 * (1 - pnorm(abs(tv)))
    fe <- data.frame(fe, test.st = tv, p.value = pv)
    re <- gblup$coefm[rownames(gblup$coefm) %in% nms, ]
    h2 <- re[, 1]/sum(re[, 1])
    re <- data.frame(re, prop.var = h2)
    cat("summary gblup analysis of trait:", gblup$name, "\n\n")
    
    cat("fixed effects equation:\n")
    print(gblup$formula, showEnv = F)
    cat("\nrandom effects equation:\n")
    print(gblup$Vformula, showEnv = F)
    
    cat("\nlog-likelihood:", gblup$llik, "converged in:", gblup$cycle, "iterations \n")
    
    cat("\nestimated fixed effects:\n")
    print(fe)
    cat("\nestimated variance components:\n")
    print(re)
    
    if (!is.null(gblup$ehat)) 
        cat("\n", length(gblup$ehat), "residual values are included\n")
    
    cat("\nbreeding value predicted for", length(gblup$uhat), "first 5 animals shown:\n")
    print(head(gblup$uhat))
    
}

#' A prints a table of estimated fixed effects and tests for gblup or summary.gblup objects
#' @title print \code{\link{gblup}} object
#' @param gblup an object of the class gblup or summary.gblup generated by \code{\link{gblup}}
#' @return a dataset with one row per fixed effect and columns for estimated value, S.E., test statistic and p-value
#' @seealso \code{\link{gblup}} 
#' @export
#'  

anova.gblup <- function(gblup) {
    nms <- all.vars(gblup$model$Vformula)
    fe <- gblup$coefm[!rownames(gblup$coef) %in% nms, ]
    tv <- fe[, 1]/fe[, 2]
    pv <- 2 * (1 - pnorm(abs(tv)))
    fe <- data.frame(fe, test.st = tv, p.value = pv)
    return(fe)
}


#' A prints a table of estimated fixed effects and tests for gblup or summary.gblup objects
#' @title print \code{\link{gblup}} object
#' @param gblup an object of the class gblup or summary.gblup generated by \code{\link{gblup}}
#' @return a dataset with one row per random effect (residuals are labeled In) and columns for estimated value, S.E. and proportion of total variance
#' @seealso \code{\link{gblup}} 
#' @export
#'  

varcomp <- function(gblup) {
      if (class(gblup) != "gblup") 
            stop("object has to be of class gblup")
      nms <- all.vars(gblup$model$Vformula)
      re <- gblup$coefm[rownames(gblup$coef) %in% nms, ]
      h2 <- re[, 1]/sum(re[, 1])
      se<-sevr(gblup)
      re <- data.frame(re, prop.var = h2,se=se)
      return(re)
  } 

#standard error of variance ratios

sevr<-function(obj){
      propvar<-obj$sigma/sum(obj$sigma)
      vrs<-diag(obj$V.CV)
      tv<-sum(vrs)
      cvs<-rowSums(obj$V.CV)
      T1<-propvar^2
      T2<-vrs/obj$sigma^2
      T3<-tv/sum(obj$sigma)^2
      T4<-2*cvs/(obj$sigma*sum(obj$sigma))
      res<-sqrt(T1*(T2+T3-T4))
      return(res)
}