R/suffDimReduct.R

In abnormally-distributed/cvreg: Cross Validation and Robust Estimation Utilities

#' sdr: methods for sufficient dimension reduction
#'
#' @description This defines the methods for the class "sdr", which is utilized for sufficient dimension reduction.
#' Sufficient dimension reduction is an approach for reducing the dimensionality of a set of covariates based on
#' the statistical principle of sufficiency. A sufficient statistic is a number that captures all of the information
#' about a data set. Likewise, a set of sufficient predictors capture all of the information about an outcome Y contained
#' in the set of covariates X. Sufficient predictors can be thought of as similar to principal components. \cr
#' \cr
#' Several methods are offered in this package. Many codes here are based upon Li (2018), with modifications
#' and improvements for more efficient calculation, and the addition of new features. Also included
#' are methods for partial sufficient dimension reduction, which estimates sufficient dimension reduction
#' subspaces integrated over a 'nuisance variable', typically a categorical predictor. In addition,
#' envelope models, regularized SIR and SAVE, structured linear regression (which allows predictors to be
#' grouped), and principle fitted components are also offered. \cr
#' \cr
#' Though not of the sdr class, the estimators in the functions \code{\link{gpcr}} and \code{\link{gpls}} are geared
#' towards the same goal as the sufficient dimension reduction approach. \cr
#'
#' @title sdr: methods for sufficient dimension reduction
#' @param x argument
#' @param ... other arguments
#' @rdname sdr
#' @export sdr
#'
#' @references Li, B. (2018). Sufficient dimension reduction: methods and applications with R. Boca Raton: CRC Press, Taylor & Francis Group.
#'
sdr <- function(x, ...){
  UseMethod("sdr")
}

#' Sliced Inverse Regression
#'
#' @description Sliced inverse regression (SIR) is a method used to compute the effective dimension reduction subspace
#' by finding a set of sufficient predictors that contain all of the information in the model matrix X
#' about the outcome Y (Li, 1991). The term inverse regression refers to the fact that the inverse regressions of E[X | Y]
#' is calculated, and the term sliced refers to the fact that Y is divided into smaller non-overlapping
#' intervals (in a manner similar to the \link[stats]{loess} algorithm).
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 5. for categorical response variables the maximum allowed is the number of response levels minus one. if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#'
#' @references Li, K-C. (1991) Sliced inverse regression for dimension reduction. Journal of the
#' American Statistical Association, 86(414), 316-327. doi: 10.1080/01621459.1991.10475035
#'
#'
#' @return an sdr object
#' @export
#'
SIR <- function(formula, data, rank = "all", slices = 5, ytype = c("numeric", "categorical")){

  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y
  ytype <- match.arg(ytype)
  if (is.factor(y) || is.character(y)){
    ytype <- "categorical"
    y.obs <- as.numeric(as.factor(y.obs))
    y <- as.numeric(as.factor(y))
  }
  if (ytype == "categorical") {
    y.obs <- as.numeric(as.factor(y.obs))
    if (min(y.obs) == 0) y.obs <- y.obs + 1
    slices <- min(slices, length(unique(y.obs)))
  }

  if (rank == "all"){
    Rank <- "all"
    rank <- ncol(X)
  }
  else {
    Rank <- rank
  }


  sdr <- .sir(X, y.obs, slices, ytype)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)

  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }

  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results = structure(list(basis = basis, sdrmat = sdrmat, predictors = predictors, evalues = evalues, evectors = evectors, y.obs = y.obs, fitted = y.pred, formula = formula, call = call, mf = model.frame(formula, data), slices = slices), class = "sdr", model = "SIR")
  return(results)
}



#' Sliced Inverse Regression 2
#'
#' @description Sliced inverse regression 2 (SIRII) is a variant of the original SIR (Li, 1991a) that
#' incorporates features of SAVE in order to improve on both (Cook & Weisberg, 1991; Li, 1991b).
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 5. for categorical response variables the maximum allowed is the number of response levels minus one. if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#'
#' @references Li, K-C. (1991a) Sliced Inverse Regression for Dimension Reduction. Journal of the
#' American Statistical Association, 86(414), 316-327. doi: 10.1080/01621459.1991.10475035 \cr
#' \cr
#' Cook, R.D., Weisberg, S. (1991) Sliced Inverse Regression for Dimension Reduction: Comment. Journal of the American Statistical Association, 86(414), 328-332 \cr
#' \cr
#' Li, K.-C. (1991b). Sliced Inverse Regression for Dimension Reduction: Rejoinder. Journal of the American Statistical Association, 86(414), 337–342. doi:10.1080/01621459.1991.10475040
#'
#'
#' @return an sdr object
#' @export
#'
SIRII =function(formula, data, slices=6, rank="all", ytype=c("numeric", "categorical")){

  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y
  ytype <- match.arg(ytype)
  if (is.factor(y) || is.character(y)){
    ytype <- "categorical"
    y.obs <- as.numeric(as.factor(y.obs))
    y <- as.numeric(as.factor(y))
  }
  if (ytype == "categorical") {
    y.obs <- as.numeric(as.factor(y.obs))
    if (min(y.obs) == 0) y.obs <- y.obs + 1
    slices <- min(slices, length(unique(y.obs)))
  }

  if (rank == "all"){
    Rank <- "all"
    rank <- ncol(X)
  }
  else {
    Rank <- rank
  }


  sdr <- .sir2(X, y.obs, slices, ytype)
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)


  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }

  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results = structure(list(basis = basis, sdrmat = sdrmat, predictors = predictors, y.obs = y.obs, fitted = y.pred, formula = formula, call = call, mf = model.frame(formula, data),  evalues = evalues, evectors = evectors, slices = slices), class = "sdr", model = "SIRII")
  return(results)
}


#' Localized Sliced Inverse Regression
#'
#' @description Localized sliced inverse regression (LSIR) is a variant of \code{\link{SIR}} (Li, 1991) that is
#' capable of sufficient dimension reduction on both linear and nonlinear subspaces. LSIR also avoids the issue
#' of degeneracy in SIR by taking into account of the local structure of the input variables
#' (Wu, Liang, & Mukherjee, 2010).
#'
#' @param formula a model formula
#' @param data a data frame
#' @param k the local neighbordhood size
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 5. for categorical response variables the maximum allowed is the number of response levels minus one. if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#'
#' @references
#' Li, K-C. (1991). Sliced inverse regression for dimension reduction. Journal of the American Statistical Association, 86(414):316-327 \cr
#' \cr
#' Wu, Q., Liang, F., & Mukherjee, S. (2010). Localized Sliced Inverse Regression. Journal of Computational and Graphical Statistics, 19(4), 843-860.
#' @return an sdr object
#' @export
#'
LSIR <- function(formula, data, k = 6, rank = "all", slices = 5, ytype = c("numeric", "categorical")){

  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]

  y.obs <- y
  ytype <- match.arg(ytype)
  if (is.factor(y) || is.character(y)){
    ytype <- "categorical"
    y.obs <- as.numeric(as.factor(y.obs))
    y <- as.numeric(as.factor(y))
  }
  if (ytype == "categorical") y.obs <- as.numeric(as.factor(y.obs))
  if (ytype == "categorical") slices <- min(slices, length(unique(y.obs))-1)

  if (rank == "all"){
    Rank <- "all"
    rank <- ncol(X)
  }
  else {
    Rank <- rank
  }

  sdr <- .lsir(X, y.obs, k, slices, ytype)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)

  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }

  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results = structure(list(basis = basis, sdrmat = sdrmat, predictors = predictors, y.obs = y.obs, fitted = y.pred, formula = formula, call = call, mf = model.frame(formula, data),evalues = evalues, evectors = evectors , slices = slices), class = "sdr", model = "LSIR")
  return(results)
}


#' Sliced Inverse Median Regression
#'
#' @description Sliced inverse regression (\code{\link{SIR}}) is a method used to compute the effective
#' dimension reduction subspace by finding a set of sufficient predictors that encapsulate all of the
#' information in the model matrix X about the outcome Y (Li, 1991). Sliced inverse median regression
#' (SIMR) is a simple modification to SIR that uses the conditional medians instead.
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 5.
#'  for categorical response variables the maximum allowed is the number of response levels minus one.
#'  if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#'
#' @references Li, K-C. (1991) Sliced inverse regression for dimension reduction. Journal of the
#' American Statistical Association, 86(414), 316-327. doi: 10.1080/01621459.1991.10475035 \cr
#' \cr
#' Christou, E. (2018). Robust dimension reduction using sliced inverse median regression.
#' Statistical Papers. doi:10.1007/s00362-018-1007-z
#' @return an sdr object
#' @export
#'
SIMR <- function(formula, data, rank = "all", slices = 5, med = c("med", "spa"), ytype = c("numeric", "categorical")){
  med <- match.arg(med)
  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y
  ytype <- match.arg(ytype)
  if (is.factor(y) || is.character(y)){
    ytype <- "categorical"
    y.obs <- as.numeric(as.factor(y.obs))
    y <- as.numeric(as.factor(y))
  }
  if (ytype == "categorical") {
    y.obs <- as.numeric(as.factor(y.obs))
    if (min(y.obs) == 0) y.obs <- y.obs + 1
    slices <- min(slices, length(unique(y.obs)))
  }

  if (rank == "all"){
    Rank <- "all"
    rank <- ncol(X)
  }
  else {
    Rank <- rank
  }

  sdr <- .simr(X, y.obs, slices, ytype, med)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)
  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }
  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results = structure(list(basis = basis, sdrmat = sdrmat, predictors = predictors, evalues = evalues, evectors = evectors, y.obs = y.obs, fitted = y.pred, formula = formula, call = call, mf = model.frame(formula, data),  slices = slices), class = "sdr", model = "SIR")
  return(results)
}



#' Regularized Sliced Inverse Regression
#'
#' @description Sliced inverse regression (SIR) is a method used to compute the effective dimension reduction subspace
#' by finding a set of sufficient predictors that contain all of the information in the model matrix X about the outcome
#' Y (Li, 1991). However, in high dimensional scenarios the solution of the SIR algorithm may be difficult and unreliable
#' to compute. Applying regularization has been shown to improve the performance and computational stability of the model
#' (Zhong et al., 2005; Bernard-Michel, et al., 2009a). It has been successfully applied to difficult problems such as
#' analyzing data from Mars (Bernard-Michel, et al., 2009b).
#'
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 5. for categorical response variables the maximum allowed is the number of response levels minus one. if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#' @param target one of "ident", "const", or "pca". "ident" uses as the shrinkage target the identity matrix,
#' which yields the ridge penalty, \eqn{(x'x + I\lambda)^{-1}}. "const" uses as the shrinakge target
#' a matrix with the geometric mean of the variances along the diagonal, and the mean covariance in each off-diagonal cell.
#' "pca"
#' @param lambda the regularization parameter. note that the value selected will have quite different impacts depending on the shrinkage target.
#'
#' @references Li, K-C. (1991) Sliced inverse regression for dimension reduction. Journal of the
#' American Statistical Association, 86(414), 316-327. doi: 10.1080/01621459.1991.10475035 \cr
#' \cr
#' Zhong, W., Zeng, P., Ma, P., Liu, J.S., & Zhu, M.Y. (2005). RSIR: regularized sliced inverse regression
#' for motif discovery. Bioinformatics, 21 22, 4169-75. \cr
#' \cr
#' Bernard-Michel, C., Gardes, L., & Girard, S. (2009a). Gaussian Regularized Sliced Inverse Regression.
#' Statistics and Computing, 19, 85-98. \cr
#' \cr
#' Bernard‐Michel, C., Douté, S., Fauvel, M., Gardes, L., & Girard, S. (2009b). Retrieval of Mars surface
#' physical properties from OMEGA hyperspectral images using regularized sliced inverse regression, J. Geophys.
#' Res., 114, E06005, doi:10.1029/2008JE003171.
#'
#' @return an sdr object
#' @export
#'
RSIR <- function(formula, data, rank = "all", slices = 5, ytype = c("numeric", "categorical"), target=c("ident", "const", "pca"), lambda = 1){

  call <- match.call()
  target <- match.arg(target)
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y
  ytype <- match.arg(ytype)
  if (is.factor(y) || is.character(y)){
    ytype <- "categorical"
    y.obs <- as.numeric(as.factor(y.obs))
    y <- as.numeric(as.factor(y))
  }
  if (ytype == "categorical") {
    y.obs <- as.numeric(as.factor(y.obs))
    if (min(y.obs) == 0) y.obs <- y.obs + 1
    slices <- min(slices, length(unique(y.obs)))
  }

  if (rank == "all"){Rank <- "all"; rank <- ncol(X)}else{Rank <- rank}

  sdr <- .rsir(X, y.obs, slices, ytype, target, lambda, rank)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)

  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
  }

  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results = structure(list(basis = basis, sdrmat = sdrmat, predictors = predictors, evalues = evalues, evectors = evectors, y.obs = y.obs, fitted = y.pred, formula = formula, call = call, mf = model.frame(formula, data),  slices = slices), class = "sdr", model = "SIR")
  return(results)
}


#' Parametric Inverse Regression
#'
#' @description Like sliced inverse regression, parametric inverse regression (PIR)
#' is a method used to compute the effective dimension reduction subspace
#' by finding a set of sufficient predictors that contain all of the information in
#' the model matrix X about the outcome Y (Bura, 1997; Bura & Cook, 2001). PIR was introduced to address the
#' occassional problem with \code{\link{SIR}} where choosing a different number of slices can yield a
#' different number of sufficient predictors, and hence can be an unreliable
#' estimator of the rank of the effective dimension reduction subspace. Instead of slicing
#' the vector Y into intervals, PIR fits inverse regressions of the form \emph{E[X | f(y)]} , where \emph{f(y)}
#' is an m-degree polynomial representation of y. This implementation offers the usual
#' orthogonal polynomials, as well as natural splines and basis splines.
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param df the tuning parameter used for the parametric inverse regression. defaults to 3.
#' @param sm the smoother function to be used. one of  cubic regression splines ("cr", the default),
#' basis splines ("bs"), or natural splines ("ns"), or orthogonal polynomials ("poly").
#' @return an sdr object
#' @export
#'
#' @references
#' Bura, E. (1997) Dimension reduction via parametric inverse regression. L1-statistical procedures and related topics: Papers from the 3rd International Conference on L1-Norm and Related Methods held in Neuchâtel, 215-228. doi:10.1214/lnms/1215454139 \cr \cr
#' Bura, E., Cook, R.D. (2001) Estimating the structural dimension of regressions via parametric inverse. Journal of the Royal Statistical Society, Series B, 63:393-410. \cr
#'
PIR <- function(formula, data, rank = "all", df = 3, sm = c("cr", "ns", "bs", "poly")){
  sm <- match.arg(sm)
  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y

  if (rank == "all"){
    Rank <- "all"; rank <- ncol(X)
  }
  else {
    Rank <- rank
  }

  sdr <- .pir(X, y.obs, df, sm)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)

  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }

  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results = structure(list(basis = basis, sdrmat = sdrmat, predictors = predictors, y.obs = y.obs, fitted = y.pred, evalues = evalues, evectors = evectors, formula = formula, call = call, mf = model.frame(formula, data)), class = "sdr", model = "PIR")
  return(results)
}



#' Sliced Average Variance Estimator
#'
#' @description A problem that \code{\link{SIR}}, \code{\link{KIR}}, and \code{\link{PIR}} face is that they cannot recover a vector in the
#' central subspace if the regression function is symmetric about 0. To improve the situation, methods
#' based on second-order conditional moments have been developed to improve the difficulty. The
#' sliced average variance estimator (SAVE) estimates the conditional expectation of the transposed crossproduct
#' of X given Y, E(XX'| Y). Like SIR, SAVE also splits Y into bins and estimates the inverse regression as a piecewise
#' function (Cook & Weisburg, 1991). Since SAVE is sensitive to directions that first order moment methods such as SIR
#' cannot detect, the central subspace estimated by SAVE compared to them will tend to be larger to or equal in size.
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 4. for categorical response variables the maximum allowed is the number of response levels minus one. if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#'
#' @return an sdr object
#' @export
#'
#' @references Cook, R.D., Weisberg, S. (1991) Sliced Inverse Regression for Dimension Reduction: Comment. Journal of the American Statistical Association, 86(414):328-332
#'
SAVE <- function(formula, data, rank = "all", slices = 4, ytype = c("numeric", "categorical")){

  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y
  ytype <- match.arg(ytype)
  if (is.factor(y) || is.character(y)){
    ytype <- "categorical"
    y.obs <- as.numeric(as.factor(y.obs))
    y <- as.numeric(as.factor(y))
  }
  if (ytype == "categorical") y.obs <- as.numeric(as.factor(y.obs))
  if (ytype == "categorical") slices <- min(slices, length(unique(y.obs))-1)

  if (rank == "all"){
    Rank <- "all"
    rank <- ncol(X)
  }
  else {
    Rank <- rank
  }

  sdr <- .save(X, y.obs, slices,ytype)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)

  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }

  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results = structure(list(basis = basis, sdrmat = sdrmat, predictors = predictors, y.obs = y.obs, fitted = y.pred,  evalues = evalues, evectors = evectors, formula = formula, call = call, mf = model.frame(formula, data)), class = "sdr", model = "SAVE")
  return(results)

}


#' Iterative Hessian Transformation Regression
#'
#' @description Iterative Hessian transformation (IHT) regression was developed as a means of relaxing the constant conditional variance
#' assumption required by \code{\link{SIR}} and \code{\link{SAVE}}. IHT only requires the linear conditional mean assumption, but like SAVE it is capable
#' of recovering more directions in the central subspace than SIR. To be precise, IHT recovers directions in the central mean
#' subspace (CMS)
#'
#' The algorithm underlying IHT is fairly simple: \cr
#'
#' 1. Calculate the covariance of X and Y, \eqn{Cov_{xy} = cov(X, Y)} and the covariance of the transposed crossproduct of X
#' and Y, \eqn{Cov_{xxy} = E[XX'Y]}. \cr
#'
#' 2. Next, calculate a matrix where the first column is Cov_{xy}, and each following column is the product of Cov_{xxy} raised to the
#' 1st to p-1th power in sequential order : \cr
#' \deqn{M = (Cov_{xy}, Cov_{xxy}Cov_{xy}, ... , Cov_{xxy}^{p-1}Cov_{xy})}
#'
#' 3. Compute \eqn{\Lambda = MM'} and take the first r eigenvectors of \eqn{\Lambda}, where r is the desired rank of the set of
#' sufficient predictors. The set of r eigenvectors are the sufficient predictors on the unit scale. \cr
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#'.
#' @return an sdr object
#' @export
#'
#'
IHT <- function(formula, data, rank = "all"){
  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y

  if (rank == "all"){
    Rank <- "all"
    rank <- ncol(X)
  }
  else {
    Rank <- rank
  }


  sdr <- .iht(X, y.obs)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)


  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }

  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)

  results = structure(list(basis = basis,
                           sdrmat = sdrmat,
                           predictors = predictors,
                           y.obs = y.obs,
                           fitted = y.pred,
                           evalues = evalues,
                           evectors = evectors,
                           formula = formula,
                           call = call,
                           mf = model.frame(formula, data)),
                      class = "sdr", model = "IHT")
  return(results)
}


#' Kernel Inverse Regression
#'
#' @description Kernel inverse regression (KIR) is very similar to \code{\link{PIR}} (Zhu & Fang, 1996). Instead of using a polynomial
#' representation of Y, Y is represented through a kernel function, k. The inverse regression E[X|Y] is calculted
#' as E[Xk(Y-y)] / E[k(Y-y)]. The covariance of E[X|Y] is calculated as the average values of the
#' transposed crossproduct of E[X|Y]-E[X], ie, Sigma = E[E[X|Y]-E[X] * E[X|Y]-E[X]']. The effective
#' dimension reduction subspace is then estimated through the generalized eigenvalue problem GEV(Sigma, cov(X)). \cr
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param kern the kernel function to be used. must be one of  "rbf", (the default), "gt, "cauchy", "spline", "laplace", "bessel", or "anova".
#' @param sigma the scale to be used for the kernel. defaults to "auto" which estimates the optimal sigma value using kernlab::sigest.
#' @param order the order to be used for the bessel function.
#' @param df the degrees of freedom to be used for the bessel kernel, anova kernel, and generalized t kernel. defaults to 3.
#' @param lambda a tuning parameter for regularization. defaults to 1e-4.
#'
#' @references
#' Zhu, L., Fang, K. (1996) Asymptotics for kernel estimate of sliced inverse regression. Ann. Statist. 24, 3, 1053-1068. doi:10.1214/aos/1032526955
#'
#' @return an sdr object
#' @export
#'
KIR <- function(formula, data, rank = "all", kern = c("rbf", "gt", "cauchy", "spline", "laplace", "bessel", "anova"), sigma = "auto", order = 1, df = 3, lambda = 1e-4){

  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y

  if (rank == "all"){
    Rank <- "all"
    rank <- ncol(X)
  }
  else {
    Rank <- rank
  }

  kern <- match.arg(kern)

  if (sigma == "auto"){
    sigma <- kernlab::sigest(formula, data, frac = 1)
    if (kern == "rbf") sigma <- sigma[2]
    if (kern == "laplace") sigma <- sigma[2] * sqrt(2)
    if (kern == "cauchy") sigma <- sigma[2] * sqrt(0.1428571)
    if (kern == "gt") sigma <- sigma[2] * sqrt((df-2) / df)
  }

  kern <- switch(kern,
    "rbf" = kernlab::rbfdot(sigma),
    "laplace" = kernlab::laplacedot(sigma),
    "bessel" = kernlab::besseldot(sigma, order, df),
    "anova" = kernlab::anovadot(sigma, df),
    "spline" = kernlab::splinedot(),
    "gt" = cvreg::gtdot(df, sigma),
    "cauchy" = cvreg::cauchydot(sigma)
  )

  sdr <- .kir(X, y.obs, kern, lambda)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)


  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }

  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results = structure(list(basis = basis,
                           predictors = predictors,
                           y.obs = y.obs,
                           fitted = y.pred,
                           evalues = evalues,
                           evectors = evectors,
                           formula = formula,
                           call = call,
                           mf = model.frame(formula, data)),
                      class = "sdr", model = "KIR")
  return(results)
}


#' Outer Product of Gradients Regression
#'
#' @description Unlike most other sufficient dimension reduction methods in this package
#' the outer product of gradients (OPG) method does not utilize the inverse regression of
#' E[X | Y] to estimate the central subspace. Rather, it uses a kernel method to calculate
#' the central mean subspace of E[Y | X] directly. \cr \cr
#'
#' 1. The algorithm procedes by using a kernel function K to calculate two matrices, A and B:
#'
#' \deqn{\hat { A }  = \sum _ { j = 1 } ^ { n }\left[\begin{bmatrix}1 \\ x_{j}-X_{i}\end{bmatrix} \begin{bmatrix}1 \\ x_{j}-X_{i}\end{bmatrix}^{T} \kappa  \left( \frac{X_{j} - X_{i}}{h} \right) \right]}
#' \cr
#' \deqn{\hat { B }  = \sum _ { j = 1 } ^ { n } \left[ \begin{bmatrix}1 \\ x_{j}-X_{i}\end{bmatrix}^{T} Y_{j} \kappa  \left(\frac{x_{j} - X_{i}}{ h } \right) \right] }
#' \cr
#' 2. Next, calculate \eqn{\beta}, which consists of the 2nd, . . . , p+1-th entries of the vector
#' \eqn{A^{-1} \cdot B}. Then sum \eqn{\beta \beta^{T}}.
#' \cr \cr
#' 3. Perform the eigendecomposition of \eqn{\sum{\beta \beta^{T}}}. The eigenvectors
#' are the estimate of the central mean subspace.
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param kern the kernel function to be used. must be one of "rbf" (the default), "gt", "cauchy, "laplace", "bessel", or "anova".
#' @param order the order to be used for the bessel function.
#' @param df the degrees of freedom to be used for the bessel kernel and anova kernel. defaults to 3.
#'
#' @return an sdr object
#' @export
#'
#' @references Xia, Y., Tong, H., Li, W., & Zhu, L. (2002). An Adaptive Estimation of Dimension Reduction Space. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 64(3), 363-410.
#'
OPG <- function(formula, data, rank = "all", kern = c("rbf", "gt", "cauchy", "laplace", "bessel", "anova"), order = 1, df = 3){

  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y

  if (rank == "all"){
    Rank <- "all"
    rank <- ncol(X)
  }
  else {
    Rank <- rank
  }

  kern <- match.arg(kern)
  kern <- switch(kern,
                 "rbf" = function(h) kernlab::rbfdot(h),
                 "laplace" = function(h) kernlab::laplacedot(h * sqrt(2)),
                 "bessel" = function(h) kernlab::besseldot(h, order, df),
                 "anova" = function(h) kernlab::anovadot(h, df),
                 "gt" = function(h) cvreg::gtdot(df, h),
                 "cauchy" = function(h) cvreg::cauchydot(h)
  )

  sdr <- .opg(X, y.obs, kern)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)

  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
  }

  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results = structure(list(basis = basis, sdrmat = sdrmat, predictors = predictors, y.obs = y.obs, fitted = y.pred,  evalues = evalues, evectors = evectors, formula = formula, call = call, mf = model.frame(formula, data)), class = "sdr", model = "OPG")
  return(results)
}


#' Contour Regression
#'
#' @description The idea behind contour regression (CR) is that when considering a surface S the
#' central dimension reduction subspace represents the directions in which the height of S changes the most.
#' In other words, the central subspace is aligned with the gradients of S.  Gradient directions are orthogonal
#' to the contours, where the height of S remains constant. Hence, the central subspace can be approximated by finding the
#' orthogonal complement to the contour directions (Li, Zha, & Chiaromonte, 2005).
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param pct the span. should be a number greater than zero and less than one.
#'
#' @return an sdr object
#' @export
#'
#' @references Li, B.; Zha, H.; Chiaromonte, F. (2005) Contour regression: A general approach to dimension reduction. Ann. Statist. 33, 4, 1580-1616. doi:10.1214/009053605000000192.
#'
CR <- function(formula, data, rank = "all", pct = 0.075){

  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y

  if (rank == "all"){
    Rank <- "all"
    rank <- ncol(X)
  }
  else {
    Rank <- rank
  }

  sdr <- .cr(X, y.obs, pct)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- X %*% basis
  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }
  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results <- structure(list(basis = basis, sdrmat = sdrmat, predictors = predictors, y.obs = y.obs, fitted = y.pred, evalues = evalues, evectors = evectors, formula = formula, call = call, mf = model.frame(formula, data)), class = "sdr", model = "CR")
  return(results)
}


#' Directional Regression
#'
#' @description Directional regression (DR) is a second-order method somewhat similar to contour regression,
#' but is simpler to compute and efficient than Contour Regression. While contour regression
#' can be seen as a method for selecting empirical directions within the central subspace,
#' DR is a method for reorganizing empirical directions (Li & Wang, 2007).
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 5. for categorical response variables the maximum allowed is the number of response levels minus one. if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#'
#' @return
#' @export
#'
#' @references
#' Li, B.; Wang, S. (2007) On directional regression for dimension reduction. Journal of the American Statistical Association, 102, 997–1008
#'
DR <- function(formula, data, rank = "all", slices = 5, ytype = c("numeric", "categorical")){
  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y
  ytype <- match.arg(ytype)
  if (is.factor(y) || is.character(y)){
    ytype <- "categorical"
    y.obs <- as.numeric(as.factor(y.obs))
    y <- as.numeric(as.factor(y))
  }
  if (ytype == "categorical") {
    y.obs <- as.numeric(as.factor(y.obs))
    if (min(y.obs) == 0) y.obs <- y.obs + 1
    slices <- min(slices, length(unique(y.obs)))
  }

  if (rank == "all"){Rank <- "all";rank <- ncol(X)}else {Rank <- rank}

  sdr <- .dr(X, y.obs, slices, ytype)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- X %*% basis

  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }

  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results <- structure(list(basis = basis, sdrmat = sdrmat, predictors = predictors, y.obs = y.obs, fitted = y.pred,  evalues = evalues, evectors = evectors, formula = formula, call = call, mf = model.frame(formula, data)), class = "sdr", model = "CR")
  return(results)
}


#' Principal Fitted Components Regression
#'
#' @description Principal fitted components (PFC) regression was developed to address two
#' flaws with principal components regression (Cook & Forzani, 2008). First, in PCR the components are
#' computed from the predictors alone and do not make apparent use of the response
#' -- they are simply the principal components of the design matrix. Second, the principal components are not invariant
#' or equivariant under full rank linear transformation of the predictors.
#'
#' @param formula a model formula
#' @param data a data frame
#' @param slices the number of slices to use. technically in this function the slices are produced by a spline function with df = slices. the default is 3.
#' @param rank the desired number of sufficient predictors to return
#' @param ytype either numeric or categorical
#' @param sfun for numeric variables, the basis function to use. must be one of "gp" (the default), "ns", "bs", "cr", or "poly",
#' which respectively correspond to gaussian process, natural splines, basis splines, cubic regression splines, and orthogonal polynomials.
#' @return an sdr object
#' @export
#'
#' @references Cook, R.D., Forzani, L. (2008) Principal Fitted Components for Dimension Reduction in Regression. Statist. Sci. 23(4), 485-501. doi:10.1214/08-STS275.
#'
PFC <- function(formula, data, slices = 3, rank = "all", ytype = c("numeric", "categorical"), sfun = c("gp", "ns", "bs", "cr", "poly")){

  call <- match.call()
  ytype <- match.arg(ytype)
  sfun <- match.arg(sfun)
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]

  if (is.factor(y) || is.character(y)){
    ytype <- "categorical"
    y <- as.numeric(as.factor(y))
  }

  sdr <- .pfcfit(X, y, slices, rank, ytype, sfun)
  basis <- sdr$basis
  sdrmat <- sdr$sdrmat
  beta <- sdr$beta
  evalues = sdr$evalues
  evectors = sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)
  if (rank != ncol(X)){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }
  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% crossprod(MX, y))
  y.pred <- as.vector(MX %*% betas)
  results = structure(list(basis = basis, beta = beta, sdrmat = sdrmat,predictors = predictors, y.obs = y, fitted = y.pred,  evalues = evalues, evectors = evectors, formula = formula, call = call, mf = model.frame(formula, data), slices = slices), class = "sdr", model = "PFC")
  return(results)
}

#' Regularized Sliced Average Variance Estimator
#'
#' @description This implements the same method of regularization used in \code{\link{RSIR}} to improve the stability
#' of \code{\link{SAVE}} in high dimensional scenarios.
#'
#' @param formula a model formula
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 5. for categorical response variables the maximum allowed is the number of response levels minus one. if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#' @param target one of "ident", "const", or "pca". "ident" uses as the shrinkage target the identity matrix,
#' which yields the ridge penalty, \eqn{(x'x + I\lambda)^{-1}}. "const" uses as the shrinakge target
#' a matrix with the geometric mean of the variances along the diagonal, and the mean covariance in each off-diagonal cell.
#' "pca"
#' @param lambda the regularization parameter. note that the value selected will have quite different impacts depending on the shrinkage target.
#' @references Cook, R.D., Weisberg, S. (1991) Sliced Inverse Regression for Dimension Reduction: Comment.
#' Journal of the American Statistical Association, 86(414):328-332 \cr
#' \cr
#' Zhong, W., Zeng, P., Ma, P., Liu, J.S., & Zhu, M.Y. (2005). RSIR: regularized sliced inverse regression
#' for motif discovery. Bioinformatics, 21 22, 4169-75. \cr
#' \cr
#' Bernard-Michel, C., Gardes, L., & Girard, S. (2009a). Gaussian Regularized Sliced Inverse Regression.
#' Statistics and Computing, 19, 85-98. \cr
#' \cr
#' Bernard‐Michel, C., Douté, S., Fauvel, M., Gardes, L., & Girard, S. (2009b). Retrieval of Mars surface
#' physical properties from OMEGA hyperspectral images using regularized sliced inverse regression, J. Geophys.
#' Res., 114, E06005, doi:10.1029/2008JE003171.
#'
#' @return an sdr object
#' @export
#'
RSAVE <- function(formula,
                  data,
                  rank = "all",
                  slices = 4,
                  ytype = c("numeric", "categorical"),
                  target = c("ident", "const", "pca"),
                  lambda = 1){

  target <- match.arg(target)
  call <- match.call()
  X <- model.matrix(formula, data)[,-1]
  y <- model.frame(formula, data)[,1]
  y.obs <- y
  ytype <- match.arg(ytype)
  if (is.factor(y) || is.character(y)){
    ytype <- "categorical"
    y.obs <- as.numeric(as.factor(y.obs))
    y <- as.numeric(as.factor(y))
  }
  if (ytype == "categorical") y.obs <- as.numeric(as.factor(y.obs))
  if (ytype == "categorical") slices <- min(slices, length(unique(y.obs))-1)

  if (rank == "all"){
    Rank <- "all"
    rank <- ncol(X)
  }
  else {
    Rank <- rank
  }

  sdr <- .rsave(X, y.obs, slices, ytype, target, lambda, rank)
  basis <- sdr$basis
  evalues <- sdr$evalues
  evectors <- sdr$evectors
  rownames(basis) <- colnames(X)
  colnames(basis) <- paste0("SP", 1:ncol(basis))
  predictors <- as.matrix(X %*% basis)

  if (Rank != "all"){
    basis = basis[,1:rank]
    predictors = predictors[,1:rank]
    evalues = evalues[1:rank]
    evectors = evectors[1:rank]
  }

  colnames(basis) <- paste0("SP", 1:ncol(predictors))
  colnames(predictors) <- paste0("SP", 1:ncol(predictors))
  MX = model.matrix(~ 0 + ., as.data.frame(predictors))
  betas = as.vector(pseudoinverse(crossprod(MX), tol = 1e-6) %*% (crossprod(MX, y.obs)))
  y.pred <- as.vector(MX %*% betas)
  results = structure(list(basis = basis, predictors = predictors, y.obs = y.obs, fitted = y.pred,  evalues = evalues, evectors = evectors, formula = formula, call = call, mf = model.frame(formula, data)), class = "sdr", model = "SAVE")
  return(results)
}





.sir=function(x,y,slices,ytype){
  p=ncol(x);n=nrow(x)
  signrt=cvreg:::matpower(mpd(cov(x)),-1/2)
  xc=t(t(x)-apply(x,2,mean))
  xst=xc%*%signrt
  if (ytype=="numeric") {
    ydis=cvreg::discretize(y,slices)
  }
  else if(ytype=="categorical"){
    ydis=y
  }
  yless=ydis;ylabel=numeric()
  for(i in 1:n) {
    if(var(yless)!=0) {
      ylabel = c(ylabel,yless[1]);yless=yless[yless!=yless[1]]
    }
  }
  ylabel=c(ylabel,yless[1])
  prob=numeric();exy=numeric()
  for(i in 1:slices) prob=c(prob,length(ydis[ydis==ylabel[i]])/n)
  for(i in 1:slices) exy=rbind(exy,apply(xst[ydis==ylabel[i],],2,mean))
  exy <- as.matrix(exy)
  sirmat=mpd(crossprod(exy,diag(prob))%*%exy)
  eig=eigen(sirmat)
  return(list(basis = signrt%*%eig$vectors, sdrmat = sirmat, evalues = eig$values, evectors=eig$vectors))
}



.sir2 <- function(x, y, slices, ytype){
  p=ncol(x)
  n=nrow(x)
  means <- colMeans(x)
  signrt=cvreg:::matpower(crossprod(sweep(x,2,means))/n, -0.50)
  xst= crossprod(t(x) - means, signrt)
  if(ytype=="numeric") {ydis=cvreg::discretize(y,slices)}
  if(ytype=="categorical") {ydis=y}
  yless=ydis
  ylabel=numeric()
  for(i in 1:n) {
    if(var(yless)!=0) {
      ylabel=c(ylabel,yless[1]);
      yless=yless[yless!=yless[1]]
    }
  }
  ylabel=c(ylabel,yless[1])
  prob=numeric()
  exy=numeric()
  for(i in 1:slices) {
    prob=c(prob,length(ydis[ydis==ylabel[i]])/n)
  }
  for(i in 1:slices) {
    exy=rbind(exy,apply(xst[ydis==ylabel[i],],2,mean))
  }
  sirmat = mpd(crossprod(exy, diag(prob)) %*% exy)
  vxy = array(0,c(p,p,slices))
  for(i in 1:slices) {
    vxy[,,i] = var(xst[ydis==ylabel[i],])
  }
  savemat=0
  for(i in 1:slices){
    savemat=savemat+prob[i]*(vxy[,,i]-diag(p))%*%(vxy[,,i]-diag(p))
  }
  siriimat=mpd(tcrossprod(savemat-sirmat, sirmat))
  eig<-eigen(siriimat,T,F)
  return(list(basis = signrt%*%eig$vectors, sdrmat = siriimat, evalues = eig$values, evectors=eig$vectors))
}


.lsir <- function(x,y,k,slices,ytype){
  scatter <- function(x, mu){
    p <- ncol(x)
    n <- nrow(x)
    s <- matrix(0, p, p)
    rowvec <- rep(0, p)
    colvec <- rep(0, p)
    for (i in 1:n){
      rowvec <- x[i,]-mu
      colvec <- t(rowvec)
      s <- s+(t(colvec)%*%t(rowvec))
    }
    return(s)
  }

  nbd <- function(tvec, tmat, k){
    n <- nrow(tmat)
    dis <- rep(0,n)
    for (i in 1:n){
      cntrd <- tvec-as.vector(tmat[i,])
      dis[i] <- sum(cntrd^2)
    }
    dorder <- order(dis)
    porder <- dorder[1:max(min(n,k),1)]
    lvec <- rep(FALSE,n)
    lvec[porder] <- TRUE
    return(lvec)
  }
  p<-ncol(x);n<-nrow(x)
  signrt<-cvreg:::matpower(mpd(cov(x)),-1/2)
  xc<-t(t(x)-apply(x,2,mean))
  xst<-xc%*%signrt
  if (ytype=="numeric") {
    ydis<-cvreg::discretize(y,slices)
  }
  else if(ytype=="categorical"){
    ydis<-y
  }
  yless=ydis;ylabel=numeric()
  for(i in 1:n) {
    if(var(yless)!=0) {
      ylabel = c(ylabel,yless[1]);
      yless=yless[yless!=yless[1]]
    }
  }
  ylabel=c(ylabel,yless[1])
  prob=numeric();exy=numeric()
  sigma <- covShrink(xc, target = "adaptive")
  exy <- array(0,c(n,p))
  for (i in 1:n){
    tgt <- setdiff(which(ydis==ydis[i]), i)
    tvec <- as.vector(xst[i,])
    tmat <- xst[tgt,]
    pidx  <- nbd(tvec, tmat, k)
    plab <- tgt[pidx]
    exy[i,] <- as.vector(colMeans(xst[plab,]))
  }
  sirmat <- mpd(scatter(exy, colMeans(xst))/n)
  eig <- eigen(sirmat,T)
  return(list(basis = signrt%*%eig$vectors, sdrmat = sirmat, evalues = eig$values, evectors=eig$vectors))
}

.rsir<-function(x, y, slices, ytype, target, lambda, rank){
  p=ncol(x); n=nrow(x)
  x = sweep(x, 2, colMeans(x))
  Sigma = cov(x)
  signrt = cvreg:::matpower(mpd(Sigma), -0.50)
  x = x%*%signrt
  if (ytype=="numeric") {
    ydis=cvreg::discretize(y,slices)
  }
  else if(ytype=="categorical"){
    ydis=y
  }
  yless=ydis;ylabel=numeric()
  for(i in 1:n) {
    if(var(yless)!=0) {
      ylabel = c(ylabel,yless[1]);yless=yless[yless!=yless[1]]
    }
  }
  ylabel=c(ylabel,yless[1])
  prob=numeric();exy=numeric()
  for(i in 1:slices) prob=c(prob,length(ydis[ydis==ylabel[i]])/n)
  for(i in 1:slices) exy=rbind(exy,apply(x[ydis==ylabel[i],],2,mean))
  exy <- as.matrix(exy)
  sirmat <- crossprod(exy,diag(prob))%*%exy


  if (target=="ident"){
    Omega = lambda*diag(p)
  }
  if (target=="const"){
    avgcov = mean(Sigma[lower.tri(Sigma)])
    avgvar = exp(mean(log(diag(Sigma))))
    shrinktarget <- matrix(avgcov, p, p)
    diag(shrinktarget) <- avgvar
    Omega = lambda*shrinktarget
  }
  if (target=="pca") {
    pcarun = eigen(Sigma)
    pcavector = pcarun$vector[,1:rank]
    Omega = array(0,c(p,p))
    for (i in 1:ncol(pcavector)){
      evec  = as.vector(pcavector[,i])
      Omega = Omega + lambda*outer(evec,evec)
    }
  }

  L = Re((Omega%*%Sigma)+diag(p))
  R = Re((Omega%*%sirmat))
  costInv = mpd(Re(pseudoinverse(L)%*%R))
  eig = eigen(costInv)
  P = eig$vectors
  evalues = eig$values
  p = ncol(P)
  Pid = crossprod(P)
  dnormval = base::norm(Re(diag(p)-Pid), type = "F")
  if (dnormval > 1e-8){
    P = try(Re(gram.schmidt(P)$Q),silent=T)
    if (class(P)=="try.error"){P=Re(qb(P,rand=F)$Q)}
  } else {
    P = try(qb(Re(P),rand=F)$Q,silent=T)
    if (class(P)=="try.error"){P=Re(P)}
  }
  n = nrow(P); p = ncol(P)
  basis = array(0,c(n,p))
  for (i in 1:p){
    evec = as.vector(P[,i])
    basis[,i] = Re(evec/sqrt(sum(evec*evec)))
  }
  evectors <- basis
  basis <- signrt%*%basis
  rownames(evectors) <- rownames(basis) <- colnames(x)
  rownames(evectors) <- colnames(basis) <- paste0("SP", 1:ncol(basis))
  result <- list()
  result$predictors <- x%*%basis
  result$evalues <- evalues
  result$evectors <- evectors
  colnames(result$predictors) <- paste0("SP", 1:ncol(basis))
  result$basis <- basis
  result$Omega <- Omega
  result$Sigma <- Sigma
  result$sdrmat <- costInv
  return(result)
}


.pir<-function(x, y, df = 3, sm = c("cr", "ns", "bs", "poly")){

  sm <- match.arg(sm)
  x = sweep(x,2,colMeans(x))
  signrt=cvreg:::matpower(mpd(cov(x)),-1/2)
  x=x%*%signrt
  ystand=(y-mean(y))/sd(y)
  if (sm == "poly") {
    f = poly(ystand, degree = df, raw = FALSE)
  }
  else if (sm == "ns"){
    f = splines::ns(ystand, df = df)
  }
  else if (sm == "bs"){
    f = splines::bs(ystand, degree = df)
  }
  else if (sm == "cr"){
    f = crSpline(ystand, df = df)
  }
  sigxf=cov(x,f);
  sigff=var(f)
  cand= mpd(sigxf%*%tcrossprod(pseudoinverse(sigff), sigxf))
  eig <- eigen(cand,T)
  return(list(basis = signrt%*%eig$vectors, sdrmat = cand, evalues = eig$values, evectors=eig$vectors))
}


.save <- function(x,y,slices,ytype){
  p=ncol(x);n=nrow(x)
  signrt=cvreg:::matpower(mpd(cov(x)),-1/2)
  xc=t(t(x)-apply(x,2,mean))
  xst=xc%*%signrt
  if(ytype=="numeric") ydis=cvreg::discretize(y,slices)
  if(ytype=="categorical") ydis=y
  yless=ydis;ylabel=numeric()
  for(i in 1:n) {if(var(yless)!=0) {ylabel=c(ylabel,yless[1]);
  yless=yless[yless!=yless[1]]}}
  ylabel=c(ylabel,yless[1])
  prob=numeric()
  for(i in 1:slices) prob=c(prob,length(ydis[ydis==ylabel[i]])/n)
  vxy = array(0,c(p,p,slices))
  for(i in 1:slices) vxy[,,i] = var(xst[ydis==ylabel[i],])
  savemat=0
  for(i in 1:slices){
    savemat=savemat+prob[i]*(vxy[,,i]-diag(p))%*%(vxy[,,i]-diag(p))
  }
  savemat = mpd(savemat)
  eig = eigen(savemat)
  return(list(basis = signrt%*%eig$vectors, sdrmat = savemat, evalues = eig$values, evectors=eig$vectors))
}



.iht <- function(x,y){
  x=sweep(x,2,colMeans(x))
  signrt=cvreg::matpower(mpd(cov(x)),-1/2);
  z=standmat(x%*%signrt);rm(x)
  szy=cov(z,y);
  szz=cov(z);
  p=dim(z)[2];
  imat=szy
  for(i in 1:(p-1)) {
    imat=cbind(imat,cvreg::matpower(szz,i)%*%szy)
  }
  imat <- mpd(tcrossprod(imat))
  eig <- eigen(imat,T)
  return(list(basis = signrt%*%eig$vectors, sdrmat=imat, evalues = eig$values, evectors=eig$vectors))
}


.kir=function(x,y,kern,lambda){
  n = length(y)
  xc=t(t(x)-apply(x,2,mean))
  signrt=cvreg::matpower(mpd(cov(x)),-1/2)
  xst=xc%*%signrt
  yst=(y-mean(y))/sd(y)
  kern = kernlab::kernelMatrix(kern, yst)
  mu=mean(c(kern%*%rep(1,n)))
  den=apply(cbind(kern%*%rep(1,n),rep(lambda*mu,n)),1,max)
  scale=eigen(kern)$values[1]
  exy=(kern%*%xst)*(1/den)
  mat=mpd(crossprod(exy))
  eig <- eigen(mat)
  return(list(basis = signrt%*%eig$vectors, sdrmat = mat, evalues = eig$values, evectors=eig$vectors))
}



.opg=function(x,y,kern){
  wls=function(X,Y,w){
    out = tcrossprod(pseudoinverse(crossprod(X,(diag(w) %*% X)),tol = 1e-6),X)%*%(diag(w)%*% Y)
    return(list(a = out[1], b = out[-1]))
  }
  standvec=function(x) return((x-mean(x))/sd(x))
  p=dim(x)[2];n=dim(x)[1]
  c0=2.34;p0=max(p,3);rn=n^(-1/(2*(p0+6)));h=c0*n^(-(1/(p0+6)))
  sig=diag(var(x));
  x=apply(x,2,standvec)
  kfun = kern(h);
  kmat=kernlab::kernelMatrix(kfun, x);
  bmat=numeric()
  for(i in 1:dim(x)[1]){
    wi=kmat[,i];xi=cbind(1,t(t(x)-x[i,]))
    bmat=cbind(bmat,wls(xi,y,wi)$b)
  }
  bmat <- mpd(tcrossprod(bmat))
  eig<-eigen(bmat)
  return(list(basis = diag(sig^(-1/2))%*%eig$vectors, evalues = eig$values, sdrmat = bmat))
}



.cr <- function(x, y, percent){
  swapindex = function(k,n){
    j=ceiling(k/n);i=k-(j-1)*n;return(c(i,j))
  }
  mu=apply(x,2,mean);
  signrt=cvreg::matpower(mpd(cov(x)),-1/2)
  z=t(t(x)-mu)%*%signrt
  n=dim(x)[1];p = dim(x)[2]
  ymat=matrix(y,n,n)
  deltay=c(abs(ymat-t(ymat)))
  singleindex=(1:n^2)[deltay < percent*mean(deltay)]
  contourmat=matrix(0,p,p)
  for(k in singleindex){
    doubleindex = swapindex(k,n)
    deltaz = z[doubleindex[1],]-z[doubleindex[2],]
    contourmat = contourmat + deltaz %*% t(deltaz)
  }
  countourmat <- mpd(contourmat)
  eig <- eigen(contourmat)
  return(list(basis = signrt%*%eig$vectors, evalues = eig$values, sdrmat = contourmat))
}


.dr <- function(x,y,slices,ytype){
  p=ncol(x);n=nrow(x)
  signrt=cvreg:::matpower(mpd(cov(x)),-1/2)
  xc=scale(x, center = TRUE, scale = FALSE)
  xst=xc%*%signrt
  if(ytype=="numeric"){ydis=cvreg::discretize(y,slices)}else if(ytype=="categorical"){ydis=y;slices=length(unique(ydis))}
  yless=ydis;ylabel=numeric();ylabel=1:slices;prob=numeric()
  for(i in 1:slices) prob=c(prob,length(ydis[ydis==ylabel[i]])/n)
  vxy = array(0,c(p,p,slices));exy=numeric()
  for(i in 1:slices) {vxy[,,i]=var(xst[ydis==ylabel[i],]);exy=rbind(exy,apply(xst[ydis==ylabel[i],],2,mean))}
  mat1 = matrix(0,p,p);mat2 = matrix(0,p,p)
  for(i in 1:slices){
    mat1 = mat1+prob[i]*(vxy[,,i]+tcrossprod(exy[i,],exy[i,]))%*%(vxy[,,i]+tcrossprod(exy[i,],exy[i,]))
    mat2 = mat2+prob[i]*tcrossprod(exy[i,],exy[i,])
  }
  out = mpd(2*mat1+2*mat2%*%mat2+2*sum(diag(mat2))*mat2-2*diag(p))
  eig <- eigen(out)
  return(list(basis = signrt%*%eig$vectors, sdrmat = out, evalues = eig$values, evectors=eig$vectors))
}


.pfcfit <- function(X, y, slices, rank, ytype, sfun){

  if (rank == "all"){
    rank = ncol(X)
  }
  if (ytype == "categorical"){
    fy <- splines::bs(y, degree = length(unique(y)))
    fy <- as.matrix(fy)
  }
  else if (ytype == "numeric"){
    fy <- switch(sfun,
                 "bs" = splines::bs(y, df = slices),
                 "ns" = splines::ns(y, df = slices),
                 "gp" = gpSpline(y, df = slices),
                 "cr" = crSpline(y, df = slices),
                 "poly" = poly(y, degree = slices)
    )
    fy <- as.matrix(fy)
  }

  iso <- function(i, Sigma, Xc, fy) {
    vnames <- colnames(X);
    ev <- eigen(Sigma);
    ev.fit <- eigen(Sigma_fit);
    all_evalues <-ev.fit$values
    evalues <- all_evalues[1:i]
    sigma2 <- Re(sum(ev$values)/p);
    basis <- Re(matrix(ev.fit$vectors[,1:i], ncol=i));
    dimnames(basis) <- list(vnames, paste("SP", 1:i, sep=""))
    Beta <- Re(tcrossprod(t(basis),Xc)%*%fy%*%pseudoinverse(mpd(crossprod(fy))));
    sigma2 <- Re((sum(ev$values)-sum(evalues))/p);
    Delta <- sigma2*diag(1,p);
    dimnames(Delta)<-list(vnames, vnames)
    return(list(Beta=Beta,basis=basis,Delta=Delta,evalues=evalues));
  }

  r <- dim(fy)[2]
  op <- dim(X); n <- op[1]; p <- op[2]
  eff.rank <- rank
  Muhat <- apply(X, 2, mean)
  Xc <- base::scale(X, TRUE, FALSE)
  P_F <- fy%*%tcrossprod(pseudoinverse(crossprod(fy)), fy)
  Sigma <- crossprod(Xc)/n
  Sigma_fit <- cov(P_F%*%X)
  Sigma_res <- Sigma - Sigma_fit
  out <- iso(eff.rank, Sigma, Xc, as.matrix(fy));
  basis <- out$basis
  beta <- out$Beta
  return(list(basis = ans, evalues = out$evalues, beta = beta, sdrmat = Sigma_fit))
}


.rsave <-function(x, y, slices, ytype, target, lambda, rank){

  p=ncol(x);n=nrow(x)
  Sigma = cov(x);
  signrt=cvreg::matpower(mpd(Sigma),-0.50)
  xc=t(t(x)-apply(x,2,mean))
  xst=xc%*%signrt
  if(ytype=="numeric") ydis=cvreg::discretize(y,slices)
  if(ytype=="categorical") ydis=y
  yless=ydis;ylabel=numeric()
  for(i in 1:n) {if(var(yless)!=0) {ylabel=c(ylabel,yless[1]);
  yless=yless[yless!=yless[1]]}}
  ylabel=c(ylabel,yless[1])
  prob=numeric()
  for(i in 1:slices) prob=c(prob,length(ydis[ydis==ylabel[i]])/n)
  vxy = array(0,c(p,p,slices))
  for(i in 1:slices) vxy[,,i] = var(xst[ydis==ylabel[i],])
  savemat=0
  for(i in 1:slices){
    savemat=savemat+prob[i]*(vxy[,,i]-diag(p))%*%(vxy[,,i]-diag(p))
  }
  if (target=="ident"){
    Omega = lambda*diag(p)
  }
  if (target=="const"){
    avgcov = mean(Sigma[lower.tri(Sigma)])
    avgvar = exp(mean(log(diag(Sigma))))
    shrinktarget <- matrix(avgcov, p, p)
    diag(shrinktarget) <- avgvar
    Omega = lambda*shrinktarget
  }
  if (target=="pca") {
    pcarun = eigen(Sigma)
    pcavector = pcarun$vector[,1:rank]
    Omega = array(0,c(p,p))
    for (i in 1:ncol(pcavector)){
      evec  = as.vector(pcavector[,i])
      Omega = Omega + lambda*outer(evec,evec)
    }
  }

  L = Re((Omega%*%Sigma)+diag(p))
  R = Re((Omega%*%savemat))
  costInv = mpd(Re(pseudoinverse(L)%*%R))
  eig = eigen(costInv)
  evalues = eig$values
  P = eig$vectors
  p = ncol(P)
  Pid = crossprod(P)
  dnormval = base::norm(Re(diag(p)-Pid), type = "F")
  if (dnormval > 1e-8){
    P = try(Re(gram.schmidt(P)$Q),silent=T)
    if(class(P)=="try.error"){P=Re(qb(P)$Q)}
  } else {
    P = try(qb(Re(P))$Q,silent=T)
    if(class(P)=="try.error"){P=Re(P)}
  }
  n = nrow(P); p = ncol(P)
  basis = array(0,c(n,p))
  for (i in 1:p){
    evec = as.vector(P[,i])
    basis[,i] = Re(evec/sqrt(sum(evec*evec)))
  }
  evectors <- basis
  basis <- signrt%*%basis
  rownames(evectors) <- rownames(basis) <- colnames(x)
  colnames(basis) <- colnames(evectors) <- paste0("SP", 1:ncol(basis))
  result <- list()
  result$predictors <- x%*%basis
  result$evectors <- evectors
  result$evalues <- evalues
  colnames(result$predictors)<- paste0("SP", 1:ncol(basis))
  result$basis <- basis
  result$Omega <- Omega
  result$Sigma <- Sigma
  result$sdrmat <- costInv
  return(result)
}



.simr=function(x,y,slices,ytype,med){
  p=ncol(x);n=nrow(x)
  if(med=="med"){
    medians <- apply(x,2,hdmedian)
  }
  else if(med=="spa"){
    medians <- L1median(x)$center
  }
  signrt=cvreg::matpower(mpd(crossprod(sweep(x,2,medians))/n),-1/2)
  xc=t(t(x)-medians)
  xst=xc%*%signrt
  if (ytype=="numeric") {
    ydis=cvreg::discretize(y,slices)
  }
  else if(ytype=="categorical"){
    ydis=y
  }
  yless=ydis;ylabel=numeric()
  for(i in 1:n) {
    if(var(yless)!=0) {
      ylabel = c(ylabel,yless[1]);yless=yless[yless!=yless[1]]
    }
  }
  ylabel=c(ylabel,yless[1])
  prob=numeric();exy=numeric()
  for(i in 1:slices) prob=c(prob,length(ydis[ydis==ylabel[i]])/n)
  if(med=="med"){
    for(i in 1:slices) exy=rbind(exy,apply(xst[ydis==ylabel[i],],2,median))
  }
  else if (med=="spa"){
    for(i in 1:slices) exy=rbind(exy,L1median(xst[ydis==ylabel[i],])$center)
  }
  exy <- as.matrix(exy)
  sirmat= mpd(crossprod(exy,diag(prob))%*%exy)
  eig=eigen(sirmat)
  return(list(basis = signrt%*%eig$vectors, sdrmat = sirmat, evalues = eig$values, evectors=eig$vectors))
}