R/Get_lambda1.R

In william-denault/ffw: A package to perform Fast functional association analysis based on wavelets

#'@title Compute lambda1 for BF null distribution
#'@description Compute the lambda parameter of the Bayes factors null distribution. We use the standard SVD R routine.
#'@param Y a vector of the variable of interest.
#'@param confounder the confounding matrix with the same number of lines as the length of Y. The intercept should not be included. If missing will generate an intercept matrix.
#'@param sigma_b value of the prior used in the Wavelet screaming.
#'@export
#'@references Quan Zhou and Yongtao Guan, On the Null Distribution of Bayes Factors in linear Regression, Journal of the American Statistical Association, 518, 2017.
#'@details Compute the lambda parameter of the Bayes factors distribution using its closed-form, which is provided in the paper of Quan Zhou and Yongtao Guan.
#'@seealso \code{\link{get_ncp}}
#'@examples \dontrun{
#'Y <- rnorm(n=1000)
#'sigma <-0.2
#'get_lambda1(Y=Y,sigma_b=sigma)
#'}


get_lambda1 <- function(Y, confounder, sigma_b)
{
  Y <- as.vector(Y)

  # INPUT CHECKS
  print("Input dimensions:")
  if(!is.numeric(Y) || length(Y)==0){
    stop("ERROR: Y is not a numeric vector")
  } else {
    print(sprintf("%i phenotypes detected", length(Y)))
    if(all(Y %in% c(0, 1, NA))){
      print("Binary phenotype detected")
    } else if(!is.vector(Y)){
      stop("ERROR: Y is not a vector. Multi-phenotype analysis not implemented yet.")
    } else {
      print("Continuous phenotype detected")
    }
  }
  if(missing(confounder)) {
    print("no covariates provided, using intercept only")
    confounder <- data.frame(confounding=rep(1,length(Y)) )
  } else if(nrow(confounder)!=length(Y)) {
    stop("ERROR: number of samples in Y and confounder does not match")
  } else {
    print(sprintf("%i covariates for %i samples detected", ncol(confounder), nrow(confounder)))
    confounder <- cbind(rep(1,length(Y)),confounder)
  }

  # Clean missing samples from all inputs
  keepY <- complete.cases(Y)
  keepC <- complete.cases(confounder)
  nonmissing_index <- which(keepY & keepC)
  if(length(nonmissing_index) != length(Y)){
    print(sprintf("Warning: %i individuals will be removed due to missingness",
                  length(Y) - length(nonmissing_index)))
  }

  Y <- Y[nonmissing_index]
  confounder <- confounder[nonmissing_index,]

  W <- as.matrix(confounder, ncol=ncol(confounder))
  L <- as.matrix(Y , ncol=ncol(Y)) #reversed regression
  n = nrow(W)
  q = ncol(W)
  L <- as.matrix(L,ncol=1)
  PW = diag(n) - W %*% solve(t(W) %*% W) %*% t(W)
  p=1 #1df test
  X= PW %*% L
  HB = X %*% solve(t(X) %*% X + diag(1/sigma_b/sigma_b,p)) %*% t(X)

  lambda <- svd(HB, nu=0, nv=0)$d[1]

  return(lambda)
}