R/vows-mgcv.R

In vows: Voxelwise Semiparametrics

#' Utility functions related to the mgcv package
#' 
#' These internal functions are used by \code{\link{semipar.mix.mp}} (but can
#' also be used more generally) to customize the implementation of B-spline
#' smoothing by \code{\link[mgcv]{gam}}. Specifically, a B-spline smooth with
#' equispaced knots can be incorporated in a call to \code{\link[mgcv]{gam}}
#' using a term of the form \code{s(x,bs="be")}, whereas knots at equally
#' spaced quantiles of the data can be specified by \code{s(x,bs="bq")}.
#' 
#' These functions are not normally called directly. For further details,
#' please see \code{\link[mgcv]{smooth.construct.ps.smooth.spec}} and
#' \code{\link[mgcv]{Predict.matrix.cr.smooth}}.
#' 
#' @aliases smooth.construct.be.smooth.spec smooth.construct.be.smooth.spec
#'          Predict.matrix.bspline.smooth
#' @param object a \code{\link[mgcv]{gam}} smooth specification object
#' generated by a term such as \code{s(x,bs="be")} or \code{s(x,bs="bq")}.
#' @param data For \code{smooth.construct.be.smooth.spec} and
#' \code{smooth.construct.bq.smooth.spec}, a list containing just the data
#' (including any \code{by} variable) required by the given term, with names
#' corresponding to \code{object$term} (and \code{object$by}). The \code{by}
#' variable is the last element. For \code{Predict.matrix.bspline.smooth}, a
#' data frame containing the values of the (named) covariates at which the
#' smooth term is to be evaluated. Exact requirements are as for
#' \code{\link[mgcv]{smooth.construct}} and
#' \code{\link[mgcv]{smooth.construct2}}.
#' @param knots a list containing any knots supplied for basis setup, in the
#' same order and with the same names as \code{data}. If \code{NULL}, a default
#' set of knots is used.
#' @return Either \code{smooth.construct.be.smooth.spec} or
#' \code{smooth.construct.bq.smooth.spec} produces an object of class
#' \code{"bspline.smooth"}; see \code{\link[mgcv]{smooth.construct}} for the
#' elements that this object will contain. \code{Predict.matrix.bspline.smooth}
#' produces a matrix mapping the coefficients for the smooth term to its values
#' at the supplied data values.
#' @author Yin-Hsiu Chen \email{enjoychen0701@@gmail.com} and Philip Reiss
#' \email{phil.reiss@@nyumc.org}
#' @examples
#' 
#' x. = rnorm(20)
#' smoo.be <- smooth.construct.be.smooth.spec(s(x), data.frame(x = x.), NULL)
#' smoo.bq <- smooth.construct.bq.smooth.spec(s(x), data.frame(x = x.), NULL)
#' Predict.matrix.bspline.smooth(smoo.bq, data.frame(x = seq(min(x.),max(x.),,100)))
#' 
#' @name vows-mgcv

#' @export smooth.construct.bq.smooth.spec
#' @export smooth.construct.be.smooth.spec
#' @export Predict.matrix.bspline.smooth

NULL

#' @rdname vows-mgcv
# Adapted from smooth.construct.ps.smooth.spec
# B-splines with knots at quantiles of data
smooth.construct.bq.smooth.spec <- function(object, data, knots) {
    if (length(object$p.order) == 1) 
        m <- rep(object$p.order, 2)
    else m <- object$p.order
    m[is.na(m)] <- 2
    object$p.order <- m
    if (object$bs.dim < 0) 
        object$bs.dim <- max(10, m[1] + 1)
    if (object$bs.dim <= m[1]) 
        stop("basis dimension too small for B-spline order")
    if (length(object$term) != 1) 
        stop("Basis only handles 1D smooths")
    x <- data[[object$term]]
    k <- knots[[object$term]]
    if (is.null(k)) {
	     bsb <- create.bspline.basis(range(x), norder=m[1]+2, nbasis = object$bs.dim, breaks = quantile(x, seq(0,1,length.out=object$bs.dim-m[1])))
	     k <- c(bsb$rangeval[1], bsb$params, bsb$rangeval[2])
    }            
	else bsb <- create.bspline.basis(range(x), norder=m[1]+2, breaks=k)

    object$X <- eval.basis(x, bsb)
    object$S <- list(getbasispenalty(bsb, m[2]))
    object$rank <- object$bs.dim - m[2]
    object$null.space.dim <- m[2]
    object$knots <- k
    object$m <- m
    class(object) <- "bspline.smooth"
    object
}

#' @rdname vows-mgcv 
# B-splines with equally-spaced knots
smooth.construct.be.smooth.spec <- function(object, data, knots) {
    if (length(object$p.order) == 1) 
        m <- rep(object$p.order, 2)
    else m <- object$p.order
    m[is.na(m)] <- 2
    object$p.order <- m
    if (object$bs.dim < 0) 
        object$bs.dim <- max(10, m[1] + 1)
    if (object$bs.dim <= m[1]) 
        stop("basis dimension too small for B-spline order")
    if (length(object$term) != 1) 
        stop("Basis only handles 1D smooths")
    x <- data[[object$term]]
    k <- knots[[object$term]]
    if (is.null(k)) {
	     bsb <- create.bspline.basis(range(x), norder=m[1]+2, nbasis = object$bs.dim)
	     k <- c(bsb$rangeval[1], bsb$params, bsb$rangeval[2])
    }            
	else bsb <- create.bspline.basis(range(x), norder=m[1]+2, breaks=k)

    object$X <- eval.basis(x, bsb)
    object$S <- list(getbasispenalty(bsb, m[2]))
    object$rank <- object$bs.dim - m[2]
    object$null.space.dim <- m[2]
    object$knots <- k
    object$m <- m
    class(object) <- "bspline.smooth"
    object
}


#' @rdname vows-mgcv
Predict.matrix.bspline.smooth<-function (object, data) {
    x = data[[object$term]]
    k = object$knots
    m = object$m
    bsb <- create.bspline.basis(range(x), norder=m[1]+2, breaks=k)
    X <- eval.basis(x, bsb)
    X
}