Find an interval constraining the rho parameter for a non linear kernel

Usage

findRhoInterval(
  tZ,
  rho_init = seq(0.01, 20, length = 300) * nrow(tZ),
  kernel = c("gaussian", "poly"),
  d = NA,
  rate_range = c(1.5, 4),
  pca_thres = 0.9,
  warning_suppress = TRUE
)

Arguments

tZ: a P x N matrix of genomic covariates (i.e., the usual data array Z transposed)
rho_init: an initial large range of possible rhos, which will be considered to see if they are reasonable tuning parameters for the kernel. Default is seq(0.01, 20, length=300)*P. See Details.
kernel: character string specifying a nonlinear kernel. Currently supported options are: "gaussian" or "poly"
d: if kernel is "poly", the polynomial power (e.g. d=2 for quadratic kernel). Default is NA.
rate_range: a vector of length 2 indicating the range of alpha in the paper. Default is c(1.5,4).
pca_thres: a number between 0 and 1 giving the threshold to be used for PCA. Default is 0.9. If NULL, no PCA is performed.
warning_suppress: logical flag. Indicating whether the warnings should be suppress during the linear model fitting step. Default is TRUE. See details.

Value

an upper and lower bound to look for rho

Details

This function will print rho_init range and the range of valid tuning parameters. If that range butts up against either the upper or lower bound of rho_init, you can rerun this function with a bigger rho_init.

Finding the right tuning parameters includes a step of fitting a linear model which can fail because some tuning parameters yield only one eigenvector. We want to eliminate those tuning parameters, so this is OK. However, in case one want to suppress (numerous) annoying warning messages, use the warning_suppress argument.

Examples


## First generate some Data
feat_m_fun <- function(X){
 sin(X[,1]+X[,2]^2)-1
}
feat_d_fun <-  function(X){
 (X[,4]-X[,5])^2/8
}
mydata <- sim_SCR_data(data_size = 400, ncol_gene_mat = 20, feat_m = feat_m_fun,
                      feat_d = feat_d_fun, mu_cen = 30, cov=0.5)

#initial range
ind_gene <- c(7:ncol(mydata))
my_rho_init <- seq(0.01, 20, length=300)*length(ind_gene)
range(my_rho_init)
#> [1]   0.2 400.0

if(interactive()){
# compute the interval for rho
rho_set <- findRhoInterval(tZ=t(mydata[,ind_gene]), rho_init = my_rho_init, kernel="gaussian")
rho_set
range(my_rho_init) # good to check that the interval produced here is strictly contained in rho_init
# otherwise, expand rho.init and rerun

#rhos <- exp(seq(log(rho_set[1]),log(rho_set[2]), length=50))
}