Nothing
################################################################################
#
# iprior: Linear Regression using I-priors
# Copyright (C) 2018 Haziq Jamil
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
#' Load the kernel matrices for I-prior models
#'
#' @param y Vector of response variables
#' @param ... Only used when fitting using non-formula, enter the variables
#' (vectors or matrices) separated by commas.
#' @param formula The formula to fit when using formula interface.
#' @param data Data frame containing variables when using formula interface.
#' @param kernel Character vector indicating the type of kernel for the
#' variables. Available choices are: \itemize{ \item{\code{"linear"} -
#' (default) for the linear kernel} \item{\code{"canonical"} - alternative
#' name for \code{"linear"}} \item{\code{"fbm"}, \code{"fbm,0.5"} - for the
#' fBm kernel with Hurst coefficient 0.5 (default)} \item{\code{"se"},
#' \code{"se,1"} - for the SE kernel with lengthscale 1 (default)}
#' \item{\code{"poly"}, \code{"poly2"}, \code{"poly2,0"} - for the polynomial
#' kernel of degree 2 with offset 0 (default)} \item{\code{"pearson" - for the
#' Pearson kernel}}} The \code{kernel} argument can also be a vector of length
#' equal to the number of variables, therefore it is possible to specify
#' different kernels for each variables. Note that factor type variables are
#' assigned the Pearson kernel by default, and that non-factor types can be
#' forced to use the Pearson kernel (not recommended).
#' @param interactions Character vector to specify the interaction terms. When
#' using formulas, this is specified automatically, so is not required. Syntax
#' is \code{"a:b"} to indicate variable \code{a} interacts with variable
#' \code{b}.
#' @param est.lambda Logical. Estimate the scale parameters? Defaults to
#' \code{TRUE}.
#' @param est.hurst Logical. Estimate the Hurst coefficients for fBm kernels?
#' Defaults to \code{FALSE}.
#' @param est.lengthscale Logical. Estimate the lengthscales for SE kernels?
#' Defaults to \code{FALSE}.
#' @param est.offset Logical. Estimate the offsets for polynomial kernels?
#' Defaults to \code{FALSE}.
#' @param est.psi Logical. Estimate the error precision? Defaults to
#' \code{TRUE}.
#' @param fixed.hyp Logical. If \code{TRUE}, then no hyperparameters are
#' estimated, i.e. all of the above \code{est.x} are set to \code{FALSE}, and
#' vice versa. If \code{NULL} (default) then all of the \code{est.x} defaults
#' are respected.
#' @param lambda Initial/Default scale parameters. Relevant especially if
#' \code{est.lambda = FALSE}.
#' @param psi Initial/Default value for error precision. Relevant especially if
#' \code{est.psi = FALSE}.
#' @param nystrom Either logical or an integer indicating the number of Nystrom
#' samples to take. Defaults to \code{FALSE}. If \code{TRUE}, then
#' approximately 10\% of the sample size is used for the Nystrom
#' approximation.
#' @param nys.seed The random seed for the Nystrom sampling. Defaults to
#' \code{NULL}, which means the random seed is not fixed.
#' @param model DEPRECATED.
#' @param one.lam Logical. When using formula input, this is a convenient way of
#' letting the function know to treat all variables as a single variable (i.e.
#' shared scale parameter). Defaults to \code{FALSE}.
#' @param train.samp (Optional) A vector indicating which of the data points
#' should be used for training, and the remaining used for testing.
#' @param test.samp (Optional) Similar to \code{train.samp}, but on test samples
#' instead.
#' @param intercept (Optional) Intercept for response variables.
#'
#' @return An \code{ipriorKernel} object which contains the relevant material to
#' be passed to the \code{iprior} function for model fitting.
#'
#' @seealso \link[=iprior]{iprior}
#'
#' @examples
#'
#' str(ToothGrowth)
#' (mod <- kernL(y = ToothGrowth$len,
#' supp = ToothGrowth$supp,
#' dose = ToothGrowth$dose,
#' interactions = "1:2"))
#' kernL(len ~ supp * dose, data = ToothGrowth) # equivalent formula call
#'
#' # Choosing different kernels
#' str(stackloss)
#' kernL(stack.loss ~ ., stackloss, kernel = "fbm") # all fBm kernels
#' kernL(stack.loss ~ ., stackloss, kernel = "FBm") # cApS dOn't MatTeR
#' kernL(stack.loss ~ ., stackloss,
#' kernel = c("linear", "se", "poly3")) # different kernels
#'
#' # Sometimes the print output is too long, can use str() options here
#' print(mod, strict.width = "cut", width = 30)
#'
#' @export
kernL <- function(
y, ..., kernel = "linear", interactions = NULL, est.lambda = TRUE,
est.hurst = FALSE, est.lengthscale = FALSE, est.offset = FALSE,
est.psi = TRUE, fixed.hyp = NULL, lambda = 1, psi = 1, nystrom = FALSE,
nys.seed = NULL, model = list(), train.samp, test.samp, intercept
) UseMethod("kernL")
#' @export
kernL.default <- function(y, ..., kernel = "linear", interactions = NULL,
est.lambda = TRUE, est.hurst = FALSE,
est.lengthscale = FALSE, est.offset = FALSE,
est.psi = TRUE, fixed.hyp = NULL, lambda = 1,
psi = 1, nystrom = FALSE, nys.seed = NULL,
model = list(), train.samp, test.samp,
intercept) {
# Checks ---------------------------------------------------------------------
if (is.list(model) & length(model) > 0) {
stop("\'model\' option is deprecated. Use the arguments directly instead. See \'?kernL\' for details.", call. = FALSE)
}
kernel <- tolower(kernel)
Xl <- list(...)
# It is common to make the mistake and type kernels instead of kernel. This
# corrects it.
Xl.kernel.mistake <- match("kernels", names(Xl))
if ("kernels" %in% names(Xl)) {
kernel <- Xl[[Xl.kernel.mistake]]
Xl[[Xl.kernel.mistake]] <- NULL
}
# Check formula
Xl.formula <- match("Xl.formula", names(Xl))
formula.method <- FALSE
if ("Xl.formula" %in% names(Xl)) {
Xl <- Xl[[Xl.formula]]
formula.method <- TRUE
}
# Get names
xname <- names(Xl)
yname <- attr(y, "yname")
# Check for probit model
if (is.factor(y)) {
tmp <- get_y_and_levels(y)
y <- tmp$y
y.levels <- tmp$levels
} else {
y.levels <- NULL
}
one.lam <- attr(Xl, "one.lam")
# Were training samples specified? -------------------------------------------
train.check <- FALSE
if (!missing(train.samp) | !missing(test.samp)) {
if (!missing(train.samp) & !missing(test.samp)) {
stop("Use either train.samp or test.samp only.", call. = FALSE)
}
if (missing(train.samp)) train.samp <- seq_along(y)[-test.samp]
if (all(train.samp %in% seq_along(y))) {
train.check <- TRUE
test.samp <- seq_along(y)[-train.samp]
y.test <- y[test.samp]
Xl.test <- lapply(Xl, function(x) {
if (is.matrix(x) | is.data.frame(x)) return(x[test.samp, , drop = FALSE])
else return(x[test.samp])
})
y <- y[train.samp]
Xl <- lapply(Xl, function(x) {
if (is.matrix(x) | is.data.frame(x)) return(x[train.samp, , drop = FALSE])
else return(x[train.samp])
})
} else {
warning("Training samples incorrectly specified.", call. = FALSE)
}
if (!is.null(one.lam)) attr(Xl, "one.lam") <- attr(Xl.test, "one.lam") <-
one.lam
}
# Get intercept --------------------------------------------------------------
if (missing(intercept)) {
y <- scale(y, scale = FALSE) # centre variables
intercept <- attr(y, "scaled:center")
} else {
y <- scale(y, scale = FALSE, center = intercept)
}
# Meta -----------------------------------------------------------------------
n <- length(y)
p <- length(Xl)
if (is.null(xname)) {
xname <- paste0("X", seq_len(p))
} else {
where.blanks <- grepl("^$", xname)
xname[where.blanks] <- paste0("X", which(where.blanks))
}
if (is.null(yname)) yname <- "y"
# For Nystrom method: Reorder data and create Xl.Nys -------------------------
if (as.numeric(nystrom[1]) > 0 & as.numeric(nystrom[1]) != n) {
if (length(nystrom) > 1) {
# This is the samples to use
nys.samp <- c(nystrom, seq_along(y)[-nystrom])
nystrom <- length(nystrom)
} else {
if (as.numeric(nystrom) == 1) nystrom <- floor(0.1 * n)
if (!is.null(nys.seed)) set.seed(nys.seed)
nys.samp <- sample(seq_along(y))
}
y.tmp <- y[nys.samp]
mostattributes(y.tmp) <- attributes(y)
y <- y.tmp
tmp <- lapply(Xl, reorder_x, smp = nys.samp)
mostattributes(tmp) <- attributes(Xl)
Xl <- tmp
Xl.nys <- lapply(Xl, reorder_x, smp = seq_len(nystrom))
mostattributes(Xl.nys) <- attributes(Xl)
nys.check <- TRUE
} else {
nys.check <- FALSE
}
# What types of kernels? -----------------------------------------------------
if (length(kernel) < p && length(kernel) > 1) {
warning(paste0("Incomplete kernel specification (not of length ", p, ")"),
call. = FALSE)
}
if (length(kernel) > p && length(kernel) > 1) {
warning(paste0("Too many kernel options specification (not of length ", p,
")"),
call. = FALSE)
}
kernels <- rep(NA, p)
suppressWarnings(kernels[1:p] <- kernel)
# The next two lines ensure that the Pearson kernel is used for factors
which.pearson <- unlist(lapply(Xl, function(x) {is.factor(x) |
is.character(x)}))
kernels <- correct_pearson_kernel(kernels, which.pearson)
# Interactions ---------------------------------------------------------------
intr <- intr.3plus <- NULL
no.int.3plus <- no.int <- 0
if (!is.null(interactions)) {
if (isTRUE(formula.method)) {
intr <- interactions$intr
intr.3plus <- interactions$intr.3plus
} else {
ind.intr.3plus <- which_intr_3plus(interactions)
intr.3plus <- interactions[ind.intr.3plus]
intr <- interactions[!ind.intr.3plus]
intr <- sapply(strsplit(intr, ":"), as.numeric)
}
if (length(intr > 0)) no.int <- ncol(intr)
}
if (length(intr.3plus) == 0) intr.3plus <- NULL
if (!is.null(intr.3plus)) {
if (!isTRUE(formula.method)) intr.3plus <- tab_intr_3plus(intr.3plus)
no.int.3plus <- ncol(intr.3plus)
}
if (isTRUE(nys.check)) {
Hl <- get_Hl(Xl, Xl.nys, kernels, lambda)
} else {
Hl <- get_Hl(Xl, list(NULL), kernels, lambda)
}
kernels <- get_kernels_from_Hl(Hl)
if (!is.null(fixed.hyp)) {
if (isTRUE(fixed.hyp)) {
est.lambda <- est.hurst <- est.lengthscale <- est.offset <- est.psi <- FALSE
}
if (!isTRUE(fixed.hyp)) {
est.lambda <- est.hurst <- est.lengthscale <- est.offset <- est.psi <- TRUE
}
}
estl <- list(est.lambda = est.lambda, est.hurst = est.hurst,
est.lengthscale = est.lengthscale, est.offset = est.offset,
est.psi = est.psi)
names(lambda) <- names(psi) <- NULL # need to clean the names otherwise weird
# things happen
param <- kernel_to_param(kernels, lambda)
poly.deg <- param$deg
thetal <- param_to_theta(param, estl, log(psi))
thetal$n.theta <- length(thetal$theta)
nystroml <- NULL
if (isTRUE(as.logical(nystrom))) {
nystroml <- list(nys.samp = nys.samp, nys.seed = nys.seed,
nys.size = as.numeric(nystrom))
}
BlockBStuff <- NULL
# Only calculate BlockBStuff when using the closed form EM algorithm.
# | | EM.closed == TRUE |
# |-----------------|------------------:|
# | poly.deg | NULL |
# | est.lambda | TRUE |
# | est.hurst | FALSE |
# | est.lengthscale | FALSE |
# | est.offset | FALSE |
# | est.psi | TRUE |
# | intr.3plus | NULL |
BlockB.cond <- (
all(is.na(poly.deg)) & !isTRUE(est.hurst) & !isTRUE(est.lengthscale) &
!isTRUE(est.offset) & (isTRUE(est.lambda) | isTRUE(est.psi)) &
!isTRUE(nys.check) & is.null(intr.3plus)
)
if (isTRUE(BlockB.cond)) {
BlockBStuff <- BlockB_fn(Hl, intr, n, p)
}
res <- list(
# Data
y = y, Xl = Xl, Hl = Hl, intercept = intercept,
# Model
kernels = kernels, which.pearson = which.pearson,
poly.deg = poly.deg, thetal = thetal, estl = estl,
intr = intr, intr.3plus = intr.3plus, nystroml = nystroml,
BlockBStuff = BlockBStuff,
# Meta
n = n, p = p, no.int = no.int, no.int.3plus = no.int.3plus,
xname = xname, yname = yname, formula = NULL, terms = NULL,
y.levels = y.levels
)
if (isTRUE(train.check)) {
res$y.test <- y.test
res$Xl.test <- Xl.test
}
res$call <- fix_call_default(match.call(), "kernL") # fix function call
class(res) <- "ipriorKernel"
res
}
#' @rdname kernL
#' @export
kernL.formula <- function(formula, data, kernel = "linear", one.lam = FALSE,
est.lambda = TRUE, est.hurst = FALSE,
est.lengthscale = FALSE, est.offset = FALSE,
est.psi = TRUE, fixed.hyp = NULL, lambda = 1,
psi = 1, nystrom = FALSE, nys.seed = NULL,
model = list(), train.samp, test.samp, intercept,
...) {
list2env(formula_to_xy(formula = formula, data = data, one.lam = one.lam),
envir = environment())
res <- kernL.default(y = y, Xl.formula = Xl, interactions = interactions,
kernel = kernel, est.lambda = est.lambda,
est.hurst = est.hurst,
est.lengthscale = est.lengthscale,
est.offset = est.offset, est.psi = est.psi,
fixed.hyp = fixed.hyp, lambda = lambda, psi = psi,
nystrom = nystrom, nys.seed = nys.seed, model = model,
train.samp = train.samp, test.samp = test.samp,
intercept = intercept, ...)
res$yname <- yname
res$formula <- formula
res$terms <- tt
res$call <- fix_call_formula(match.call(), "kernL") # fix function call
res
}
#' @export
print.ipriorKernel <- function(x, units = "auto", standard = "SI", ...) {
tmp <- expand_Hl_and_lambda(x$Hl, seq_along(x$Hl), x$intr, x$intr.3plus)
# if (isTRUE(x$probit)) {
# cat("Categorical response variables\n")
# } else if (is.ipriorKernel_nys(x)) {
# cat("Nystrom kernel approximation ()\n")
# }
n <- x$n
if (is.ipriorKernel_cv(x)) {
n.test <- length(x$y.test)
n <- paste0(n + n.test, " (", n, " train + ", n.test, " test)")
}
cat("Sample size:", n, "\n")
cat("No. of covariates:", length(x$Xl), "\n")
# cat("No. of interactions:", x$no.int + x$no.int.3plus, "\n")
cat("Object size: ")
print(object.size(x), units = units, standard = standard)
cat("\n")
cat("Kernel matrices:\n")
for (i in seq_along(tmp$Hl)) {
cat("", i, print_kern(tmp$Hl[[i]], ...), "\n")
}
cat("\n")
cat("Hyperparameters to estimate:\n")
if (x$thetal$n.theta > 0)
cat(paste(names(x$thetal$theta), collapse = ", "))
else
cat("none")
cat("\n")
cat("\n")
methods <- c("direct", "em", "canonical", "mixed", "fixed")
poss.method <- NULL
for (i in seq_along(methods)) {
suppressWarnings(tmp <- iprior_method_checker(x, methods[i]))
poss.method <- c(poss.method, names(which(tmp)))
}
poss.method <- gsub("em.closed", "em", poss.method)
poss.method <- gsub("em.reg", "em", poss.method)
poss.method <- gsub("nystrom", "direct", poss.method)
if (is.nystrom(x)) {
poss.method <- paste(poss.method, "(Nystrom)")
}
if (is.categorical(x)) {
poss.method <- c(poss.method, "iprobit (recommended)")
}
poss.method <- paste0(unique(poss.method), collapse = ", ")
cat("Estimation methods available:\n")
cat(poss.method)
cat("\n")
}
print_kern <- function(x, ...) {
# Helper function to prettify the print output of kernel matrices.
#
# Args: x is a kernel matrix obtained from one of the kern_x() functions.
# Additional ... are passed to str().
#
# Returns: Prettified kernel matrix str() print ouput.
kern.type <- attr(x, "kernel")
res <- capture.output(str(x, ...))[1]
res <- gsub(" num", kern.type, res)
res
}
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.