Nothing
#' Kernel Partial Least Squares Fit
#'
#' This function computes the Partial Least Squares fit. This algorithm scales
#' mainly in the number of observations.
#'
#' We first standardize \code{X} to zero mean and unit variance.
#'
#' @param X matrix of predictor observations.
#' @param y vector of response observations. The length of \code{y} is the same
#' as the number of rows of \code{X}.
#' @param m maximal number of Partial Least Squares components. Default is
#' \code{m}=ncol(X).
#' @param compute.jacobian Should the first derivative of the regression
#' coefficients be computed as well? Default is \code{FALSE}
#' @param DoF.max upper bound on the Degrees of Freedom. Default is
#' \code{min(ncol(X)+1,nrow(X)-1)}.
#' @return \item{coefficients}{matrix of regression coefficients}
#' \item{intercept}{vector of regression intercepts} \item{DoF}{Degrees of
#' Freedom} \item{sigmahat}{vector of estimated model error} \item{Yhat}{matrix
#' of fitted values} \item{yhat}{vector of squared length of fitted values}
#' \item{RSS}{vector of residual sum of error} \item{covariance}{\code{NULL}
#' object.} \item{TT}{matrix of normalized PLS components}
#' @author Nicole Kraemer, Mikio L. Braun
#' @seealso \code{\link{linear.pls.fit}},
#' \code{\link{pls.cv}},\code{\link{pls.model}}, \code{\link{pls.ic}}
#' @references Kraemer, N., Sugiyama M. (2011). "The Degrees of Freedom of
#' Partial Least Squares Regression". Journal of the American Statistical
#' Association 106 (494)
#' \url{https://www.tandfonline.com/doi/abs/10.1198/jasa.2011.tm10107}
#'
#' Kraemer, N., Braun, M.L. (2007) "Kernelizing PLS, Degrees of Freedom, and
#' Efficient Model Selection", Proceedings of the 24th International Conference
#' on Machine Learning, Omni Press, 441 - 448
#' @keywords multivariate
#' @examples
#'
#' n<-50 # number of observations
#' p<-5 # number of variables
#' X<-matrix(rnorm(n*p),ncol=p)
#' y<-rnorm(n)
#'
#'
#' pls.object<-kernel.pls.fit(X,y,m=5,compute.jacobian=TRUE)
#'
#'
#'
#' @export kernel.pls.fit
kernel.pls.fit=function (X, y, m=ncol(X), compute.jacobian=FALSE,DoF.max=min(ncol(X)+1,nrow(X)-1)) {
p <- ncol(X)
n <- nrow(X)
m<-min(m,DoF.max)
#Beta <- matrix(, p, m) # matrix of regression coefficients
#W <- V <- Beta
#dW <- dBeta <- dV <- array(dim = c(m, p, n))
X0<-X
y0<-y
# scaling of the data
mean.y <- mean(y)
y <- scale(y, scale = FALSE)
mean.X <- apply(X, 2, mean)
sd.X <- apply(X, 2, sd)
sd.X[sd.X == 0] = 1
X <- X - rep(1, nrow(X)) %*% t(mean.X)
X <- X/(rep(1, nrow(X)) %*% t(sd.X))
K<-X%*%t(X)
#dcoefficients=NULL
yhat<-vector(length=m+1)
TT <- matrix(, n, m)
Yhat <- matrix(, n, m)
Alpha <- matrix(0, n, m)
Gamma <- matrix(, n, m)
DoF = NULL
dYhat = NULL
if (compute.jacobian == TRUE) {
dYhat = array(dim = c(m, n, n))
dtildeT <- array(dim = c(m, n, n))
dT <- dtildeT
DoF = vector(length = m)
}
for (i in 1:m) {
if (i == 1) {
ri <- y
gi <- ri
ti <- tildeti <- K %*% y
if (compute.jacobian == TRUE) {
dT[i, , ] <- dtildeT[i, , ] <- K
dT[i, , ] <- dnormalize(ti, dT[i, , ])
dtildeT[i, , ] <- dT[i, , ]
}
dummy <- normalize(ti, gi)
ti <- dummy$v
gi <- dummy$w
tildeti <- ti
TT[, i] <- ti
Gamma[, i] <- gi
Alpha[, i] <- Gamma[, i] * sum(ti * y)
Yhat[, i] = vvtz(ti, y)
if (compute.jacobian == TRUE) {
ti <- as.vector(ti)
dYhat[i, , ] = dvvtz(ti, y, dT[i, , ], diag(n))
DoF[i] <- sum(diag(dYhat[i, , ]))
}
}
if (i > 1) {
ri <- y - Yhat[, i - 1]
tildeti <- K %*% ri
if (compute.jacobian == TRUE) {
dtildeT[i, , ] <- K %*% (diag(n) - dYhat[i -
1, , ])
}
TTi <- TT[, 1:(i - 1), drop = FALSE]
ti <- tildeti - vvtz(TTi, tildeti)
dummy <- rep(0, n)
for (r in 1:(i - 1)) {
dummy <- dummy + sum(tildeti * TT[, r]) * Gamma[,r]
}
gi <- ri - dummy
if (compute.jacobian == TRUE) {
dT[i, , ] = dtildeT[i, , ] - dvvtz(TTi, tildeti,
dT[1:(i - 1), , , drop = FALSE], dtildeT[i,
, ])
dT[i, , ] = dnormalize(ti, dT[i, , ])
}
dummy <- normalize(ti, gi)
ti <- dummy$v
gi <- dummy$w
TT[, i] <- ti
Gamma[, i] <- gi
Yhat[, i] <- Yhat[, i - 1] + vvtz(ti, y)
#yhat[i]<-sum(Yhat[,i]^2)
Alpha[, i] <- Alpha[, i - 1] + Gamma[, i] * sum(TT[,
i] * y)
if (compute.jacobian == TRUE) {
ti <- as.vector(ti)
dYhat[i, , ] <- dYhat[i - 1, , ] + dvvtz(ti,
y, dT[i, , ], diag(n))
DoF[i] <- sum(diag(dYhat[i, , ]))
}
}
}
Yhat<-cbind(rep(0,n),Yhat)
Yhat<-Yhat+mean.y
yhat<-apply(Yhat^2,2,sum)
Alpha<-cbind(rep(0,n),Alpha)
if (compute.jacobian==TRUE){
DoF<-c(0,DoF)+1
}
if (compute.jacobian==FALSE){
DoF<-1:(m+1)
}
RSS<-vector(length=m+1)
for (i in 1:(m+1)) {
res <- Yhat[, i] - y0
RSS[i]<-sum(res^2)
}
if ((compute.jacobian) == TRUE) {
dummy<-array(1/n,dim=c(m+1,n,n))
dummy[2:(m+1),,]<-dummy[2:(m+1),,] + dYhat
dYhat<-dummy
denominator <- vector(length = m+1)
for (i in 1:(m+1)) {
denominator[i] <- sum(diag((diag(n) - dYhat[i, ,
]) %*% ((diag(n) - t(dYhat[i, , ])))))
}
sigmahat <- sqrt(RSS/denominator)
}
if (compute.jacobian==FALSE){
sigmahat<-sqrt(RSS/(n-DoF))
}
coefficients<-t(X)%*%Alpha
coefficients = coefficients/(sd.X %*% t(rep(1, m+1)))
intercept <- rep(mean.y, m+1) - t(coefficients) %*% mean.X
intercept<-as.vector(intercept)
DoF[DoF>DoF.max]=DoF.max
covariance=NULL
return(list(coefficients=coefficients, intercept=intercept,DoF = DoF, Yhat=Yhat,yhat=yhat,
RSS = RSS, sigmahat = sigmahat,covariance=covariance, TT = TT))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.