R/EBmvFR_workhorse.R

In stephenslab/susiF.alpha: Sum of Single Functions

#' @title Empirical Bayes multivariate functional regression
#'
#' @description Empirical Bayes multivariate functional regression
#'
#' @details Empirical Bayes multivariate functional regression
#'
#' @param obj an object of class EBmvFR
#' @param W a list in which element D contains matrix of wavelet d coefficients and
#'  element C contains the vector of scaling coefficients
#'
#' @param X matrix of size n by p contains the covariates
#' @param verbose If \code{verbose = TRUE}, the algorithm's progress,
#' and a summary of the optimization settings are printed to the
#' console.
#' @param lowc_wc list of wavelet coefficients that exhibit too little variance
#' @param indx_lst list generated by gen_wavelet_indx for the given level of resolution
#'
#' @param tol a small, non-negative number specifying the convergence
#' tolerance for the IBSS fitting procedure. The fitting procedure
#' will halt when the difference in the variational lower bound, or
#' \dQuote{ELBO} (the objective function to be maximized), is less
#' than \code{tol}.
#'
#' @param maxit Maximum number of IBSS iterations.
#'
#' @param init_pi0_w starting value of weight on null compoenent in mixsqp
#'  (between 0 and 1)
#' @param control_mixsqp list of parameter for mixsqp function see  mixsqp package
#' @param  cal_obj logical if set as TRUE compute ELBO for convergence monitoring
#' @param nullweight numeric value for penalizing likelihood at point mass 0
#' (useful in small sample size)
#' @param  max_step_EM see susiF function
#' @export

EBmvFR.workhorse <- function(obj,
                             W,
                             X,
                             tol,
                             lowc_wc,
                             init_pi0_w ,
                             control_mixsqp ,
                             indx_lst,
                             nullweight,
                             cal_obj ,
                             verbose,
                             maxit,
                             max_step_EM){


  Y_f      <-  cbind( W$D,W$C)
  # numerical value to check breaking condition of while
  check <- 3*tol

  if(verbose){
    print("Initialization done")
  }
  ##### Start While -----
  iter <- 1


  while( (check >tol & iter <maxit))
  {
    for( j in 1:ncol(X))
    {

      if(verbose){
        print(paste("Fitting effect ", j,", iter" ,  iter ))
      }
      obj   <-  fit_effect.EBmvFR (obj = obj,
                                          j          = j,
                                          X          = X,
                                          D          = W$D,
                                          C          = W$C,
                                          indx_lst   = indx_lst,
                                          lowc_wc    = lowc_wc)

    }#end for l in 1:L  -----



    sigma2    <- estimate_residual_variance(obj,Y=Y_f,X)
   # print(sigma2)
    obj <- update_residual_variance(obj, sigma2 = sigma2 )

    obj <- update_prior_EBmvFR( obj            = obj,
                                max_step       = max_step_EM,
                                espsilon       = 0.0001,
                                init_pi0_w     = init_pi0_w ,
                                control_mixsqp = control_mixsqp,
                                indx_lst       = indx_lst,
                                lowc_wc        = lowc_wc,
                                nullweight     = nullweight   )

    obj <- test_stop_cond( obj = obj,
                                 check      = check,
                                 cal_obj    = cal_obj,
                                 Y          = Y_f,
                                 X          = X,
                                 D          = W$D,
                                 C          = W$C,
                                 indx_lst    = indx_lst)

    #print(obj$alpha)
    #print(obj$ELBO)
    check <- obj$check



    iter <- iter +1


  }#end while

  obj$niter <- iter
  return(obj)
}