R/operation_on_prior.R

In stephenslab/susiF.alpha: Sum of Single Functions

######################## SuSiF operation on prior ###########################
#' @title Initialize the prior
#
#' @description generate list of object corresponding to the parameters of the prior set for analysis
#
#' @return an object of the class "normal", "mixture_normal" or "mixture_normal_per_scale"
#
#' @export
#' @keywords internal
init_prior <- function(  ...)
  UseMethod("init_prior")

#' @title Initialize the prior
#
#' @description generate list of object corresponding to the parameters of the prior
#
#' @param Y  functional phenotype, matrix of size N by size J. The underlying algorithm uses wavelet which assume that J is of the form J^2. If J not a power of 2, susif internally remaps the data into grid of length 2^J
#
#' @param X matrix of size n by p in
#
#' @param prior Three choice are available "normal", "mixture_normal", "mixture_normal_per_scale"
#
#' @param v1 a vector of ones of length equal to nrow(Y)
#' @param control_mixsqp list of parameter for mixsqp function see  mixsqp package
#' @param nullweight numeric value for penalizing likelihood at point mass 0 (should be between 0 and 1)
# (usefull in small sample size)
#' @param indx_lst list generated by \code{\link{gen_wavelet_indx}} for the given level of resolution, used only with class  mixture_normal_per_scale
#
#' @param lowc_wc wavelet coefficient with low count to be discarded
#' @param gridmult numeric used to control the number of component used in the mixture prior (see ashr package
#  for more details). From the ash function:  multiplier by which the default grid values for mixsd differ by one another.
#   (Smaller values produce finer grids.). Increasing this value may reduce computational time
#' @param ind_analysis, optional, specify index for the individual to be analysied, allow analyis data with different entry with NA
#' if a vector is provided, then we assume that the entry of Y have NA at the same place, if a list is provide
#
#' @return an object of the class "normal", "mixture_normal" or "mixture_normal_per_scale"
#
#' @importFrom ashr ash
#
#' @export
#' @keywords internal
init_prior.default <- function(Y,
                               X,
                               prior,
                               v1 ,
                               indx_lst,
                               lowc_wc,
                               control_mixsqp,
                               nullweight ,
                               gridmult=sqrt(2),
                               ind_analysis,

                               max_SNP_EM=100,
                               max_step_EM=1,
                               cor_small=FALSE,
                               tol_null_prior=0.001,... )
{


  if (cor_small){
    df = nrow(X)
  }else{
    df =NULL
  }


  if(missing(ind_analysis)){


    temp <- cal_Bhat_Shat(Y, X ,v1,
                          lowc_wc=lowc_wc
                         )   ## Speed Gain would be good to call directly cal_Bhat_Shat in the ash function

  }else{
    temp <- cal_Bhat_Shat(Y[ind_analysis,], X [ind_analysis,] ,v1,
                          lowc_wc=lowc_wc
                         )
  }
  if( prior == "mixture_normal")
  {


    G_prior <- list()
    if( !is.null(lowc_wc)){


      set.seed(1)
      betahat <- c(max(abs(temp$Bhat[,-lowc_wc])), sample(temp$Bhat[,-lowc_wc], size = min( prod(dim(temp$Bhat[,-lowc_wc])), 5000)) )
      set.seed(1)
        sdhat <-c(0.01 , sample(temp$Shat[,-lowc_wc], size = min( prod(dim(temp$Shat[,-lowc_wc])), 5000)) )

      t_ash <-  ashr::ash(betahat, sdhat,#ash can take quite some space
                          mixcompdist ="normal",
                          outputlevel=0,
                          gridmult =gridmult)

      t_ash$fitted_g$pi <- c(0.8 , rep( 0.2/  (length(t_ash$fitted_g$pi)-1),  (length(t_ash$fitted_g$pi)-1)))



      G_prior[[1]]  <-  t_ash


    }else {



      set.seed(1)
      betahat <- c(max(abs(temp$Bhat)), sample(temp$Bhat, size = min( prod(dim(temp$Bhat)), 5000)) )
      set.seed(1)
       sdhat <- c(0.01, sample(temp$Shat, size = min( prod(dim(temp$Shat)), 5000)) )

      t_ash <-  ashr::ash(betahat, sdhat,#ash can take quite some space
                          mixcompdist ="normal",
                          outputlevel=0,
                          gridmult =gridmult)

      t_ash$fitted_g$pi <- c(0.8 , rep( 0.2/  (length(t_ash$fitted_g$pi)-1),  (length(t_ash$fitted_g$pi)-1)))


      G_prior[[1]]   <-  t_ash
    }

    attr(G_prior, "class")  <- "mixture_normal"
  }

  if( prior == "mixture_normal_per_scale")
  {
    if( !is.null(lowc_wc)){
      set.seed(1)
      betahat <- c(max(abs(temp$Bhat[,-lowc_wc])), sample(temp$Bhat[,-lowc_wc], size = min( prod(dim(temp$Bhat[,-lowc_wc])), 50000)) )
      set.seed(1)
       sdhat <-c(0.01 , sample(temp$Shat[,-lowc_wc], size = min( prod(dim(temp$Shat[,-lowc_wc])), 50000)) )

      t_ash <-  ashr::ash(betahat, sdhat,#ash can take quite some space
                          mixcompdist ="normal",
                          outputlevel=0,
                          gridmult =gridmult)

      t_ash$fitted_g$pi <- c(0.8 , rep( 0.2/  (length(t_ash$fitted_g$pi)-1),  (length(t_ash$fitted_g$pi)-1)))



    }else {
      set.seed(1)
      betahat <- c(max(abs(temp$Bhat)), sample(temp$Bhat, size = min( prod(dim(temp$Bhat)), 50000)) )
      set.seed(1)
       sdhat <- c(0.01, sample(temp$Shat, size = min( prod(dim(temp$Shat)), 50000)) )

      t_ash <-  ashr::ash(betahat, sdhat,#ash can take quite some space
                          mixcompdist ="normal",
                          outputlevel=0,
                          gridmult =gridmult)

      t_ash$fitted_g$pi <- c(0.8 , rep( 0.2/  (length(t_ash$fitted_g$pi)-1),  (length(t_ash$fitted_g$pi)-1)))


    }

    G_prior  <- rep(list(t_ash),length(indx_lst))
    #first log2(Y_f)+1 element of G_prior   are ash prior fitted per level coefficient on var 1
    # element in  (log2(Y_f)+2):  2*( log2(Y_f)+1)   of G_prior   are ash prior fitted per level coefficient on var 2

    attr(G_prior, "class") <- "mixture_normal_per_scale"
  }

  tpi_k <-  EM_pi(G_prior=G_prior,
                  Bhat=temp$Bhat,
                  Shat= temp$Shat,
                  indx_lst=indx_lst,
                  max_step = max_step_EM,
                  espsilon = 0.0001,
                  init_pi0_w =1,
                  control_mixsqp=control_mixsqp,
                  lowc_wc=lowc_wc,
                  nullweight=nullweight,
                  max_SNP_EM=max_SNP_EM,
                  df=df,
                  tol_null_prior =tol_null_prior )$tpi_k

  G_prior <- update_prior(G_prior , tpi_k)
  return(list(G_prior=G_prior,
              tt=temp)
  )
}





#' @title Get mixture proportion for mixture prior
#'
#' @description Get mixture proportion for mixture   prior
#'
#' @param G_prior mixture normal prior
#' @param \dots Other arguments.
#' @return vector of mixture proportion
#'
#' @export
#' @keywords internal
get_pi_G_prior <- function(G_prior, ...)
  UseMethod("get_pi_G_prior")

#' @rdname get_pi_G_prior
#'
#' @method get_pi_G_prior mixture_normal
#'
#' @export get_pi_G_prior.mixture_normal
#'
#' @export
#' @keywords internal

get_pi_G_prior.mixture_normal <- function(G_prior, ...)
{
  out <- G_prior[[1]]$fitted_g$pi
  class(out)  <- "pi_mixture_normal"
  return(out)
}

#' @rdname get_pi_G_prior
#'
#' @method get_pi_G_prior mixture_normal_per_scale
#'
#' @export get_pi_G_prior.mixture_normal_per_scale
#'
#' @export
#' @keywords internal
get_pi_G_prior.mixture_normal_per_scale <- function(G_prior, ...)
{
  out <- lapply(G_prior, function(x) x$fitted_g$pi)
  class(out) <- "pi_mixture_normal_per_scale"
  return(out)
}




#' @title Get mixture proportion for mixture normal prior
#'
#' @description Add description here.
#'
#' @param tpi  object of class pi_mixture_normal
#'
#' @return numeric between 0 an 1
#'
#' @export
#' @keywords internal


get_pi0 <- function(tpi, ...)
  UseMethod("get_pi0")


#' @rdname get_pi0
#'
#' @method get_pi0 pi_mixture_normal
#'
#' @export get_pi0.pi_mixture_normal
#'
#' @export
#' @keywords internal
get_pi0.pi_mixture_normal  <- function(tpi, ...)
{
  out <- tpi[1]
  return(out)
}




#' @rdname get_pi0
#'
#' @method get_pi0 pi_mixture_normal_per_scale
#'
#' @export get_pi0.pi_mixture_normal_per_scale
#'
#' @export
#' @keywords internal
get_pi0.pi_mixture_normal_per_scale  <- function(tpi, ...)
{
  out <-    (unlist(lapply(tpi, function(y) y[[1]][1]) ))
  return(out)
}



#' @title Get mixture standard deviations for mixture normal prior
#'
#' @description Add description here.
#'
#' @param G_prior mixture normal prior
#'
#' @return vector of standard deviations
#'
#' @export
#' @keywords internal

get_sd_G_prior <- function(G_prior , ...)
  UseMethod("get_sd_G_prior")


#' @rdname get_sd_G_prior
#'
#' @method get_sd_G_prior mixture_normal
#'
#' @export get_sd_G_prior.mixture_normal
#'
#' @export
#' @keywords internal
get_sd_G_prior.mixture_normal <- function(G_prior, ...)
{
  out <- G_prior[[1]]$fitted_g$sd
  class(out) <- "sd_mixture_normal"
  return(out)
}



#' @rdname get_sd_G_prior
#'
#' @method get_sd_G_prior mixture_normal_per_scale
#'
#' @export get_sd_G_prior.mixture_normal_per_scale
#'
#' @export
#' @keywords internal
get_sd_G_prior.mixture_normal_per_scale <- function(G_prior, ...)
{
  out <- lapply(G_prior, function(x) x$fitted_g$sd)
  class(out) <- "sd_mixture_normal_per_scale"
  return(out)
}





#' @title Update mixture proportion for mixture normal prior
#'
#' @description Add description here.
#'
#' @param G_prior a prior of class "mixture_normal" or a prior of class "mixture_normal_per_scale"
#'
#' @param tpi a vector of proportion of class "pi_mixture_normal" resp "pi_mixture_normal_per_scale"
#'
#' @param \dots Other arguments.
#' @return a prior of class "mixture_normal" or a prior of class "mixture_normal_per_scale"
#'
#' @export
#' @keywords internal
#'

update_prior <- function(G_prior, tpi, ...)
  UseMethod("update_prior")


#' @rdname update_prior
#'
#' @method update_prior mixture_normal
#'
#' @export update_prior.mixture_normal
#'
#' @export
#' @keywords internal
update_prior.mixture_normal <- function(G_prior, tpi, ...)
{
  if( class(tpi)[1]=="pi_mixture_normal"){
   G_prior[[1]]$fitted_g$pi <- tpi
  }else{
    stop("Error: tpi is not of class pi_mixture_normal,\n please compute tpi using generic functions m_step or get_pi_G_prior")
  }

  return(G_prior)
}


#' @rdname update_prior
#'
#' @method update_prior mixture_normal_per_scale
#'
#' @export update_prior.mixture_normal_per_scale
#'
#' @export
#' @keywords internal

update_prior.mixture_normal_per_scale <- function(G_prior, tpi, ...)
{
  if( class(tpi)[1]=="pi_mixture_normal_per_scale"){
    out <- mapply(update_ash_pi ,G_prior, tpi, SIMPLIFY = FALSE)
    class(out ) <- "mixture_normal_per_scale"
  }else{
    stop("Error: tpi is not of class pi_mixture_normal_per_scale,\n please compute tpi using generic functions m_step or get_pi_G_prior")
  }

  return(out)
}

#' @title Update ash object mixture proportion
#'
#' @description  Update ash object mixture proportion
#'
#' @param G a ash object
#'
#' @param tpi a vector of proportion
#'
#' @return an ash object with updated mixture proportion
#'
#' @export
#' @keywords internal
update_ash_pi<- function(G , tpi)
{

  if( length(G$fitted_g$pi )==length(tpi))
  {
    G$fitted_g$pi <- tpi
  }else{
    stop("Error: when updating ash object length of new mixture proportion
    \nlonger than the mixture proportion in the ash object    ")

  }

  return(G)
}