Nothing
#' Restricted cubic splines with cure
#'
#' Function for computing the basis matrix for restricted cubic splines which are constant beyond the last knot
#'
#' @param x Values to evaluate the basis functions in.
#' @param df Degrees of freedom. One can supply \code{df} rather than knots; \code{cbc} then chooses \code{df + 1} knots
#' at suitably chosen quantiles of \code{x} (which will ignore missing values) and adds an additional knot at the 95th
#' quantile of \code{x}.
#' @param knots Chosen knots for the spline.
#' @param ortho Logical. If \code{TRUE} orthogonalization of the basis matrix is carried out.
#' @param R.inv Matrix or vector containing the values of the R matrix from the QR decomposition of the basis matrix.
#' This is used for making new predictions based on the initial orthogonalization.
#' Therefore the default is \code{NULL}.
#' @param intercept Logical. If \code{FALSE}, the intercept of the restricted cubic spline is removed.
#' @return A matrix with containing the basis functions evaluated in \code{x}.
#' @references Andersson T.M.-L., et al. (2011) Estimating and modelling cure in population-based cancer
#' studies within the framework of flexible parametric survival models.
#' \emph{BMC Medical Research Methodology}, 11:96.
#' @export
cbc <- function(x, df = NULL, knots = NULL, ortho = FALSE, R.inv = NULL, intercept = FALSE){
if(is.null(knots)){
nIknots <- df
qs <- seq.int(0, 1, length.out = nIknots)
knots <- quantile(x, sort(c(qs, 0.95)))
}
nk <- length(knots)
b <- matrix(nrow = length(x), ncol = nk - 1)
knots_rev <- rev(knots)
if (nk > 0) {
b[, 1] <- 1
}
if (nk > 2) {
for (j in 2:(nk - 1)) {
lam <- (knots_rev[nk - j + 1] - knots_rev[1])/(knots_rev[nk] - knots_rev[1])
b[, j] <- pmax(knots_rev[nk - j + 1] - x, 0)^3 - lam * pmax(knots_rev[nk] - x, 0)^3 -
(1 - lam) * pmax(knots_rev[1] - x, 0)^3
}
}
if(!intercept) b <- b[,-1, drop = F]
if(ortho){
if(is.null(R.inv)){
qr_decom <- qr(b)
b <- qr.Q(qr_decom)
R.inv <- solve(qr.R(qr_decom))
} else{
R.inv <- matrix(R.inv, nrow = ncol(b))
b <- b %*% R.inv
}
} else {
R.inv <- diag(ncol(b))
}
nx <- names(x)
dimnames(b) <- list(nx, 1L:ncol(b))
a <- list(knots = knots, ortho = ortho, R.inv = R.inv, intercept = intercept)
attributes(b) <- c(attributes(b), a)
class(b) <- c("cbc", "basis", "matrix")
b
}
#Predict function associated with bsx.
#' @export
#' @method predict cbc
predict.cbc <- function (object, newx, ...)
{
if (missing(newx))
return(object)
a <- c(list(x = newx), attributes(object)[c("knots", "ortho",
"R.inv", "intercept")])
do.call("cbc", a)
}
#Additional function needed to fix the knot location in cases where df is only specified
# #' @export
# makepredictcall.cbc <- function (var, call)
# {
# if (as.character(call)[1L] != "cbc")
# return(call)
# at <- attributes(var)[c("knots", "ortho", "R.inv", "intercept")]
# xxx <- call[1L:2]
# xxx[names(at)] <- at
# xxx
# }
# #' @export
# makepredictcall <- function(var, call) UseMethod("makepredictcall")
# dbasis.cure <- function(knots, x, ortho = TRUE, R.inv = NULL, intercept = TRUE){
# nk <- length(knots)
# b <- matrix(nrow = length(x), ncol = nk - 1)
# knots_rev <- rev(knots)
# if (nk > 0) {
# b[, 1] <- 0
# }
# if (nk > 2) {
# for (j in 2:(nk - 1)) {
# lam <- (knots_rev[nk - j + 1] - knots_rev[1])/(knots_rev[nk] - knots_rev[1])
# b[, j] <- - 3 * pmax(knots_rev[nk - j + 1] - x, 0)^2 + 3 * lam * pmax(knots_rev[nk] - x, 0)^2 +
# 3 * (1 - lam) * pmax(knots_rev[1] - x, 0)^2
# }
# }
#
# if(!intercept) b <- b[, -1]
#
# if(ortho){
# b <- b %*% R.inv
# }
# b
# }
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.