weimea: Weighted Mean Analysis

#' Testing the relationship between weighted-mean of species attributes with sample attributes
#' 
#' Function calculating relationship between weighted-mean of species attributes and sample attributes and performing standard (row based), modified (column based), or max test of significance.
#' 
#' @name test_cwm
#' @param cwm An object of the class \code{cwm}. 
#' @param env Vector or matrix with variables. See details.
#' @param method Statistical method used to analyse the relationship between cwm (of class \code{cwm}) and env (sample attributes); partial match to \code{'cor'}, \code{'lm'} and \code{'aov'}. 
#' @param wcor Logical; should the correlation be weighted by rowsums of \code{com} (extracted from \code{cwm}?). Default \code{FALSE}.
#' @param wstand Logical; should the variables in correlation be first weighted-standardized? Default \code{FALSE}. If \code{wstand = TRUE}, both \code{env} and \code{traits} are weighted standardized prior to calcuation; weights are derived from \code{com} matrix (extracted from \code{cwm}); \code{env} is weighted by rowsums of \code{com}, while \code{traits} are weighted by colsums of \code{com}.
#' @param wreg Logical; should weights be used in the regression (\code{method = 'lm'})? If \code{TRUE}, weighted least squares are calculated instead of ordinary least squares.
#' @param wsamp Either \code{NULL} (default) or a numeric vector; if \code{wreg = TRUE}, weights provided here are used to calculate weighted least square regression; if \code{wsamp = NULL}, \code{rowSums (com)} are used as weights. NOT IMPLEMENTED YET.
#' @param dependence Should \code{cwm} be dependent variable and \code{env} independent (\code{'cwm ~ env'}), or opposite? Applicable only for \code{method = 'lm'}. Partial match to \code{'cwm ~ env'} and \code{'env ~ cwm'}.
#' @param perm Number of permutations.
#' @param test Vector of character values. Which test should be conducted? Partial match to \code{'standard'} or \code{'rowbased'} for standard (row based) permutation test, \code{'modified'} or \code{'colbased'} for modified (column based) permutation test, \code{'max'} for max test (selecting the higher from \code{rowbased} and \code{colbased} result), and \code{'all'} including all three tests. See \code{Details}.
#' @param adjustP Logical, default FALSE. Should be the P-values adjusted? If \code{adjustP = TRUE}, the last column in the results is adjusted by method selected in \code{p.adjust.method}.
#' @param p.adjust.method A string indicating the method of P-value adjustement, see \code{\link{p.adjust.methods}} for possible choices.
#' @param parallel NULL (default) or integer number. Number of cores for parallel calculation of modified permutation test. Maximum number of cores should correspond to number of available cores on the processor.
#' @param x,object object of the class \code{"cwm"}, generated by function \code{cwm}.
#' @param digits number of digits reported by \code{print} method on object of \code{"cwm"} class (default is 3).
#' @param missing.summary Logical; should be the summary of values missing in \code{env}, \code{cwm} and \code{traits} be printed along to the result of \code{test_cwm} analysis? Default is \code{TRUE}.
#' @param eps.Pvalue Values of P below this threshold will be printed as \code{< [eps]} in the output.
#' @param signif.stars Logical; if TRUE, P-values are additionally encoded visually as 'significance stars' in order to help scanning of long coefficient tables. It defaults to the show.signif.stars slot of \code{\link{options}}.
#' @param alpha,line,cex.lab,par.mar,box.col,box.lwd Graphical parameters for \code{plot} function.
#' @param ... Other arguments for \code{print}, \code{summary}, \code{coef} or \code{plot} functions (some not implemented yet).

#' @export
#' @examples 
#' data (vltava)
#' 
#' # Traits vs environment (tested by max test)
#' CWM_traits <- cwm (com = vltava$herbs$spe, traits = vltava$herbs$traits)
#' re_traits <- test_cwm (cwm = CWM_traits, env = vltava$env[,c('pH', 'COVERE32')],
#'  method = 'lm', adjustP = TRUE)
#' re_traits
#' plot (re_traits)
#'
#' # Ellenberg indicator values vs assignment of plots into groups (by cluster analysis)
#' # (tested by modified (column-based) permutation test)
#' CWM_ell <- cwm (com = vltava$spe, traits = vltava$civ[,1:5]) 
#' re_ell <- test_cwm (cwm = CWM_ell, env = as.factor (vltava$env$GROUP), test = 'modif',
#'  method = 'aov', adjustP = TRUE)
#' re_ell
#' plot (re_ell)
#' @details
#' Currently implemented statistical methods: \code{'cor'}, \code{'lm'} and \code{'aov'}. For fourth corner use \code{\link{test_fourth}}.
#'
#' Argument \code{env} can be vector or matrix with one column. Only in the case of linear regression (\code{method = 'lm'}) it is possible to use matrix with several variables, which will all be used as independent variables in the model. For ANOVA and Kruskal-Wallis test, make sure that 'env' is \code{factor} (warning will be returned if this is not the case, but the calculation will be conducted). 
#' 
#' Specific issue related to weighted mean is the case of missing species attributes. In current implementation, species with missing species attributes are removed from sample x species matrix prior to permutation of species attributes among species. 
#' 
#' The \code{plot} functions will plot pairwise relationships between CWM and environmental variables, as a scatterplots in case of \code{method = 'cor'} and \code{'lm'} and as boxplot in case of \code{method = 'aov'}. Significant relationships are highlighted in the figure by colorful border; the P-values used for this highlighting are the last one listed in the summary output. If you don't like this behaviour, limit the analysis to a single test only - in that case this test will be used to highlight the significant results.
#' @return  Function \code{cwm} returns list of the class \code{"cwm"} (with \code{print} and \code{summary} methods), which contains the following components:
#' \itemize{
#'  \item \code{call} Call to the function.
#'  \item \code{out} Matrix with analysis results (coefficients, statistics, P-values).
#'  \item \code{miss} Matrix with counts of missing values in \code{env}, \code{cwm} and \code{traits}.
#'  \item \code{param} List with the setting of the function parameters (arguments).
#' }
#' @seealso \code{\link{cwm}}, \code{\link{snc}}
#' @importFrom graphics box par plot title axis boxplot layout lines plot.window
#' @importFrom stats lm predict
#' @importFrom stats model.matrix
#' @export
test_cwm <- function (cwm, env, method = c('cor'), wcor = FALSE, wstand = FALSE, wreg = FALSE, wsamp = NULL, dependence = "cwm ~ env", perm = 499, test = "max", parallel = NULL, p.adjust.method = 'holm', adjustP = FALSE)
{
  CALL <- match.call ()
  METHOD <- c('cor', 'lm', 'aov')
  TEST <- c('none', 'parametric', 'standard', 'rowbased', 'modified', 'colbased', 'max', 'all')
  DEPENDENCE <- c("cwm ~ env", "env ~ cwm")
  p.adjust.method <- match.arg (p.adjust.method, p.adjust.methods)
  method <- match.arg (method, METHOD)
  test <- match.arg (test, TEST, several.ok = T)
  if ('none' %in% test) test <- 'none'
  if ('all' %in% test) test <- c('parametric', 'rowbased', 'colbased', 'max')
  dependence <- match.arg (dependence, DEPENDENCE)
  
  if (method == 'cor' & any (unlist (lapply (env, is.factor)))) stop ("Object env contains factors; to calculate correlation of factor and CWM is not meaningful. Consider using method = 'aov' instead.")
  if (method == 'lm' & any (unlist (lapply (env, is.factor)))) stop ("Object env contains factors; to calculate correlation of factor and CWM is not meaningful. Consider using method = 'aov' instead.")
  if (method == 'aov' & any (unlist (lapply (env, is.numeric)))) stop ("Object env contains one or more variables which are not factors; to calculate ANOVA is not meaningful.")
  
  if (!is.cwm (cwm) & ("modified" %in% test || "colbased" %in% test || "max" %in% test || "all" %in% test)) stop ("Object cwm must be of 'cwm' class")
  env <- as.data.frame (env)
  # if ((method == 'lm') & ('max' %in% test || 'colbased' %in% test || 'modified' %in% test) & (dependence == 'cwm ~ env') & (ncol (env) > 1)) stop ("Column-based (modified) and max permutation test is not available for multiple linear regression with more than one predictor variable (env)")  # This can be fixed and removed
  com <- attr (cwm, 'com')
  traits <- attr (cwm, 'traits')

# correlation - function 'cor' ----
  if (method == 'cor')
  {
    vars.names <- expand.grid (traits = colnames (traits), env = colnames (env), stringsAsFactors = FALSE)

    if (is.null (parallel)) {
      res <- apply (vars.names, 1, FUN = function (var12) test_cwm_cor (e = env[, var12['env'], drop = TRUE], L = as.matrix (com), t = traits[, var12['traits'], drop = TRUE], test = test, perm = perm, wcor = wcor, wstand = wstand))
    }
       else 
      {
        cl <- parallel::makeCluster (parallel)
        parallel::clusterExport (cl, varlist = c('vars.names', 'com', 'traits', 'env', 'test', 'perm'), envir = environment ())
        parallel::clusterEvalQ (cl, eval (call ('library', 'weimea')))
        res <- parallel::parApply (cl, vars.names, 1, FUN = function (var12) test_cwm_cor (e = as.matrix (env[, var12['env'], drop = F]), L = as.matrix (com), t = as.matrix (traits[, var12['traits'], drop = F]), test = test, perm = perm, wcor = wcor, wstand = wstand))
      }
    names (res) <- apply (vars.names, 1, FUN = function (x) paste (x[1], x[2], sep = ' / '))
   }

  # linear regression - function 'lm' ----
  if (method == 'lm')
  {
    vars.names <- expand.grid (traits = colnames (traits), env = colnames (env), stringsAsFactors = FALSE)
    
    if (is.null (parallel)) 
      res <- apply (vars.names, 1, FUN = function (var12) test_cwm_lm (e = as.matrix (env[, var12[2], drop = FALSE]), L = as.matrix (com), t = as.matrix (traits[, var12[1], drop = FALSE]), test = test, dependence = dependence, perm = perm, wstand = wstand))
    else
    {
      cl <- parallel::makeCluster (parallel)
      parallel::clusterExport (cl, varlist = c('vars.names', 'com', 'traits', 'env', 'test', 'perm'), envir = environment ())
      parallel::clusterEvalQ (cl, eval (call ('library', 'weimea')))
      res <- parallel::parApply (cl, vars.names, 1, FUN = function (var12) test_cwm_lm (e = as.matrix (env[, var12[2], drop = FALSE]), L = as.matrix (com), t = as.matrix (traits[, var12[1], drop = FALSE]), test = test, dependence = dependence, perm = perm, wstand = wstand))
    }
    names (res) <- apply (if (dependence == 'cwm ~ env') vars.names[,1:2] else vars.names[,2:1], 1, FUN = function (x) paste (x[1], x[2], sep = ' ~ '))
  }

  dummy <- function(df) {  
    as.dummy <- function (var) 
    {
      res <- model.matrix (~ as.matrix  (var) - 1)
      lev <- levels (var)
      colnames (res) <- lev
      return (res)
    }
    temp_res <- lapply (df, FUN = function (column) if (is.factor (column)) as.dummy (column) else column)
    do.call (cbind.data.frame, temp_res)
  }
  
  # analysis of variance - function 'aov' ----
  if (method == 'aov')
  {
    vars.names <- expand.grid (traits = colnames (traits), env = colnames (env), stringsAsFactors = FALSE)
    if (is.null (parallel)) 
      res <- apply (vars.names, 1, FUN = function (var12) test_cwm_aov (e = as.matrix (dummy (env[, as.character (var12[2]), drop = FALSE])), L = as.matrix (com), t = as.matrix (traits[, var12[1], drop = FALSE]), test = test, perm = perm, wstand = wstand))
    else
    {
      cl <- parallel::makeCluster (parallel)
      parallel::clusterExport (cl, varlist = c('vars.names', 'com', 'traits', 'env', 'test', 'perm'), envir = environment ())
      parallel::clusterEvalQ (cl, eval (call ('library', 'weimea')))
      res <- parallel::parApply (cl, vars.names, 1, FUN = function (var12) test_cwm_aov (e = as.matrix (env[, var12[2], drop = FALSE]), L = as.matrix (com), t = as.matrix (traits[, var12[1], drop = FALSE]), test = test,perm = perm, wstand = wstand))
    }
    names (res) <- apply (vars.names[,1:2], 1, FUN = function (x) paste (x[1], x[2], sep = ' ~ '))
  }
  
  res.miss <- lapply (res, FUN = function (x) x[(length (x)-2):length (x)])
  res.coef <- lapply (res, FUN = function (x) x[1:(length (x)-3)])
  if (method == 'lm') res.coef <- lapply (res.coef, FUN = function (x) c(interc = x$coef[1], interc_se = x$stderr[1], coef = x$coef[2], coef_se = x$stderr[2], x[-1:-2]) )
  res.coef.df <- do.call (rbind.data.frame, res.coef)
  test.choice <- c('parametric', 'standard', 'rowbased', 'modified', 'colbased', 'max')
  P.test <- c ('P_par', 'P_row', 'P_row', 'P_col', 'P_col', 'P_max')[test.choice %in% test]
  res.out <- cbind (res.coef.df [,!substr (colnames (res.coef.df), 1, 2) %in% 'P_'], res.coef.df[, P.test, drop = F])  # deletes all P values and adds only those specificed in "test" argument
  if (adjustP && test != 'none') {
    res.out <- cbind (res.out, p.adjust (res.out[, ncol (res.out)], method = p.adjust.method))
    colnames (res.out)[ncol (res.out)] <- paste (colnames (res.out)[ncol (res.out)-1], 'adj', sep = '_')
  }
  
  result <- list (call = CALL, out = res.out, miss = res.miss, param = list (method = method, wcor = wcor, wstand = wstand, dependence = dependence, perm = perm, test = test, P.test = P.test, adjustP = adjustP, p.adjust.method = p.adjust.method), var.names = list (env = colnames (env), traits = colnames (traits)), orig.data = list (cwm = cwm, env = env))
  class (result) <- 'testCWM'
  return (result)
}

#' @rdname test_cwm
#' @export
print.testCWM <- function (x, digits = max(3, getOption("digits") - 3), missing.summary = FALSE, eps.Pvalue = 0.001, signif.stars = getOption("show.signif.stars"), ...)
{
  cs.ind <- switch (x$param$method, 
                    'cor' = 1,
                    'lm' = 1:4,
                    'aov' = 1:2,
                    'fourthcorner' = 1:2)
  tst.ind <- switch (x$param$method,
                     'cor' = 2,
                     'lm' = c(5:7, 10),
                     'aov' = c(5:7, 10),
                     'fourthcorner' = 0)
  df.ind <- switch (x$param$method,
                    'cor' = 3:4,
                    'lm'= 8:9,
                    'aov' = c(3:4, 8:9),
                    'fourthcorner' = 3:4)
  cat("\nCall:\n", paste(deparse(x$call), sep = "\n", collapse = "\n"), 
      "\n\n", sep = "")
  printCoefmat2 (x$out, digits = digits, cs.ind = cs.ind, tst.ind = tst.ind, df.ind = df.ind, eps.Pvalue = 0.001, signif.stars = signif.stars, ...)
  if (missing.summary) printMissum (x$miss)
  
}

printMissum <- function (x){
  miss <- lapply (x, FUN = function (a) lapply (a, length))
  miss <- as.matrix (do.call (rbind.data.frame, miss))
  colnames (miss) <- c('env', 'traits', 'samples')
  if (sum (miss) == 0) cat("\nThere are no missing env, traits or sample values.\n") else
  {
    cat("\nCount of missing values:\n")
    print.default (miss)
  }
  invisible (miss)
  }

#' @rdname test_cwm
#' @export
coef.testCWM <- function (object, ...)
  return (object$out)
 
join_mat <- function (mat1, mat2, incl1names = TRUE, incl2names = TRUE) {
  mat12 <- lapply (1: min (ncol (mat1), ncol (mat2)), FUN = function (i) {
    mat12_temp <- cbind (mat1[,i], mat2[,i])
    if (incl1names || incl2names) colnames (mat12_temp) <- c (if (incl1names) colnames (mat1)[i] else "", if (incl2names) colnames (mat2)[i] else "")
    mat12_temp})
  mat12 <- as.matrix (data.frame (mat12, check.names = FALSE, fix.empty.names = TRUE))
  colnames (mat12)[seq (2, ncol (mat12), by = 2)] <- ''  # this makes sure that P.values *** will not have colnames (a bit stupid solution here)
  return (mat12)
}

printCoefmat2 <- function (x, digits = max(3L, getOption("digits") - 2L),
                           signif.stars = getOption("show.signif.stars"), 
                           signif.legend = signif.stars,
                           dig.tst = max(1L, min(5L, digits - 1L)),
                           cs.ind = NULL, tst.ind = NULL,
                           df.ind = NULL, zap.ind = integer(),
                           eps.Pvalue = 0.001,
                           na.print = "NA", ...) 
{
  d <- dim(x)
  if (is.null(d) || length(d) != 2L) stop("'x' must be coefficient matrix/data frame")
  nc <- d[2L]
  xm <- data.matrix(x)
  Cf <- array("", dim = c(d[1], sum (cs.ind > 0) + sum (tst.ind > 0) + sum (df.ind > 0)), dimnames = list (dimnames(xm)[[1]], dimnames (xm)[[2]][1:(sum (cs.ind > 0) + sum (tst.ind > 0) + sum (df.ind > 0))]))
  ok <- !(ina <- is.na(xm))
  for (i in zap.ind) xm[, i] <- zapsmall(xm[, i], digits)
  if (length(cs.ind)) {
    acs <- abs(coef.se <- xm[, cs.ind, drop = FALSE])
    if (any(ia <- is.finite(acs))) {
      digmin <- 1 + if (length(acs <- acs[ia & acs != 
                                          0])) 
        floor(log10(range(acs[acs != 0], finite = TRUE)))
      else 0
      Cf[, cs.ind] <- format(round(coef.se, max(1L, digits - 
                                                  digmin)), digits = digits)
    }
  }
  if (length(tst.ind)) 
    Cf[, tst.ind] <- format(round(xm[, tst.ind], digits = dig.tst), 
                            digits = digits)
  ok[, tst.ind] <- FALSE
  if (length(df.ind))
    Cf[, df.ind] <- format (xm[, df.ind])
  
  dec <- getOption("OutDec")
  if (dec != ".") 
    x1 <- chartr(dec, ".", x1)
  if (any(ina)) 
    Cf[ina] <- na.print
  P.values.log <- substr (colnames (xm), 1, 2) == 'P_'
  if (any (P.values.log)) P.values <- colnames (xm) [P.values.log] else P.values <- NULL
  if (!is.null (P.values))
   {
    if (!is.logical(signif.stars) || is.na(signif.stars)) {
      warning("option \"show.signif.stars\" is invalid: assuming TRUE")
      signif.stars <- TRUE
    }
    pv <- xm[, P.values, drop = F]
    pv_f <- matrix (vector (mode = 'character', length = ncol (pv)*nrow (pv)), ncol = ncol (pv), nrow = nrow (pv), dimnames = dimnames (pv))
    for (co in 1:ncol (pv_f)) pv_f[,co] <- format.pval(pv[,co], digits = dig.tst, eps = eps.Pvalue)
    signif.stars <- signif.stars && any (pv < 0.1)
    if (is.na (signif.stars)) signif.stars <- FALSE  # in case there are no P-values at all (test = 'none')
    if (signif.stars) {
      Signif <- symnum(as.matrix (pv), corr = FALSE, na = FALSE, 
                       cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), 
                       symbols = c("***", "**", "*", ".", " "))
      Cf <- cbind(Cf, join_mat (pv_f, format (Signif), incl2names = FALSE))
    } else Cf <- cbind (Cf, pv_f)
  }
  else signif.stars <- FALSE
  
  Cf <- cbind (rownames (Cf), Cf)
  rownames (Cf) <- 1:nrow (Cf)
  
  print.default(Cf, quote = FALSE, right = TRUE, na.print = na.print, 
                ...)
  if (signif.stars && signif.legend) {
    if ((w <- getOption("width")) < nchar(sleg <- attr(Signif, 
                                                       "legend"))) 
      sleg <- strwrap(sleg, width = w - 2, prefix = "  ")
    cat("---\nSignif. codes:  ", sleg, sep = "", fill = w + 
          4 + max(nchar(sleg, "bytes") - nchar(sleg)))
  }
  invisible(x)
}

#' @rdname test_cwm
#' @export
plot.testCWM <- function (x, alpha = 0.05, line = NA, cex.lab = 1.5, par.mar = c(0.5, 0.5, 0.5, 0.5), box.col = c('blue', 'red'), box.lwd = 2, ...){
  
  abline_lm <- function (x, y) {
    lm0 <- lm (y ~ x)
    lines (y = predict (lm0, newdata = data.frame (x = c(min (x), max (x)))), x = c(min (x), max (x)), xpd = F)
  }
  
  opar <- par (c('mfcol', 'mar', 'mgp', 'oma', 'xpd'))  # include all graphical par set globally
  
  trait.names <- x$var.names$traits
  env.names <- x$var.names$env
  CWM <- x$orig.data$cwm
  ENV <- x$orig.data$env
  
  no_traits <- length (trait.names)
  no_env <- length (env.names)
  
  par (oma = c(3, 3, 1, 1))  # outer margin
  par (mgp = c(2, 0.3, 0))  # where title, axis labels and axis line should be plotted
  layout_matrix <- matrix (seq (1, length (trait.names)*length (env.names)), nrow = if (x$param$method == 'lm' & x$param$dependence == 'env ~ cwm') no_env else no_traits, byrow = if (x$param$method == 'lm' & x$param$dependence == 'env ~ cwm') TRUE else FALSE)
  layout (layout_matrix)
  par (mar = par.mar, xpd = NA)
  
  test <- x$out[,max (which (substr (colnames (x$out), 1, 2) %in% 'P_'))]
  
  traits.env <- expand.grid (traits = trait.names, env = env.names)
  
  if (x$param$method == 'lm' & x$param$dependence == 'env ~ cwm'){
    left_label <- vector ('logical', length = no_traits*no_env)
    left_label[seq (1, no_traits*no_env, by = ncol(layout_matrix))] <- TRUE
    
    bottom_label <- vector ('logical', length = no_traits*no_env)
    bottom_label [(no_traits*(no_env-1)+1):(no_traits*no_env)] <- TRUE
  } else {
    left_label <- vector ('logical', length = no_traits*no_env)
    left_label[1:nrow(layout_matrix)] <- TRUE
    
    bottom_label <- vector ('logical', length = no_traits*no_env)
    bottom_label [seq (nrow(layout_matrix), no_traits*no_env, by = nrow(layout_matrix))] <- TRUE
  }
  
  
  traits.env <- cbind (traits.env, left_label = left_label, bottom_label = bottom_label)
  if (!is.null (test))
    if (x$param$method == 'cor') traits.env <- cbind (traits.env, test = test, r = x$out$r) else
      if (x$param$method == 'lm') traits.env <- cbind (traits.env, test = test, r = x$out$coef) else
        if (x$param$method == 'aov') traits.env <- cbind (traits.env, test = test, r = x$out$F)
  
  
  for (i in seq (1, nrow (traits.env)))
  {
    te <- traits.env[i,]
    if (x$param$method == 'lm' & x$param$dependence == 'env ~ cwm') {
      h <- CWM[[te$traits]]  # horizontal variable
      v <- ENV[,te$env]  # vertical variable
    } else {
      v <- CWM[[te$traits]] 
      h <- ENV[,te$env]  
    }
      
    if (is.factor (ENV[,te$env]))  # boxplot for factors, scatter for quantitative
      boxplot (v ~ h, pch = 16, ann = F, axes = F) else
        plot (v ~ h, pch = 16, ann = F, axes = F, asp = if (x$param$method == 'cor') diff (range (h))/diff (range (v)) else NA)
    if (is.factor (h))  # make sure x-axis in boxplot will be correctly plotted
      axis (1, at = seq_along (levels (h)), labels = if (te$bottom_label) levels (h) else
        rep ('', length (levels (h))), tck = 0.05) else axis (1, labels = ifelse (te$bottom_label, TRUE, FALSE), tck = 0.05)
    axis (2, labels = ifelse (te$left_label, TRUE, FALSE), tck = 0.05)
    box (col = 'grey')
    if (x$param$method == 'lm' & x$param$dependence == 'env ~ cwm') 
      title (xlab = ifelse (te$bottom_label, as.character (te$traits), ''), ylab = ifelse (te$left_label, as.character (te$env), ''), line = line, cex.lab = cex.lab) else
        title (xlab = ifelse (te$bottom_label, as.character (te$env), ''), ylab = ifelse (te$left_label, as.character (te$traits), ''), line = line, cex.lab = cex.lab)
    if (!is.null (test)) {
      bc <- ifelse (as.numeric (te$r) < 0, box.col[1], box.col[2])
      if (as.numeric (te$test) <= alpha){
        box (col = bc, lwd = box.lwd)
        if (x$param$method == 'lm') abline_lm (x = h, y = v)  # for lm draw reg line
      } 
    }
  }
  par (opar)
  invisible ()
}