R/gmrf_smooth_spec.R

In Faskally/gmrf: Lightweight helper for fitting GMRF models in mgcv

#' Construct a smoother using a GMRF smoothing prior
#'
#' This function is used internally in mgcv when fitting a smoother 
#' of type gmrf.  More details to be added.
#'
#' @param object a smooth specification object, usually generated by a term
#'              `s(...,bs="gmrf", penalty = list(Q = ...))'. `x' is a factor
#'          variable giving labels for the nodes in the GMRF graph, and the `xt'
#'          argument is obligatory: see details.
#' @param data a list containing just the data (including any `by' variable)
#'          required by this term, with names corresponding to
#'          `object$term' (and `object$by'). The `by' variable is the
#'          last element.
#' @param knots not used.
#' @return An object of class `"gmrf.smooth"'.
#' @examples
#' require(gmrf)
#' n <- 100
#' idx <- 1:n
#' idy <- c(5:40, 60:n)
#' set.seed(64)
#' Q <- getQrw(n, order = 2)
#' x <- simQ(exp(1) * Q)
#' # simulate the first 3/4, say
#' y <- x + rnorm(n) * 3.5
#' y <- y[idy]
#' ## set up variables for smoothing
#' rownames(Q) <- colnames(Q) <- idx
#' ## fit an RW2 smoother with restricted df
#' g1 <- gam(y ~ s(idy, bs = "gmrf", xt = list(penalty = Q), k = length(y)-1), method="REML")
#' summary(g1)
#' plot(idy, y, xlim = range(idx), ylim = range(x,y))
#' lines(idx, x, col = "blue", lwd = 2)
#' pred <- predict(g1, newdata = list(idy = idx), se = TRUE)
#' lines(idx, pred$fit, col = "red", lwd = 2)
#' lines(idx, pred$fit + 2*pred$se.fit, col = "red", lty = 2)
#' lines(idx, pred$fit - 2*pred$se.fit, col = "red", lty = 2)
#' @export
smooth.construct.gmrf.smooth.spec <- function(object, data, knots) {
    k <- factor(rownames(object$xt$penalty), levels = rownames(object$xt$penalty))
    x <- data[[object$term]]
    x <- factor(x, levels = levels(k))
    #k <- factor(levels(x), levels = levels(x))
    if (object$bs.dim < 0) object$bs.dim <- length(levels(k))
    if (object$bs.dim > length(levels(k))) stop("GMRF basis dimension set too high")
    object$X <- model.matrix(~ x - 1, )

    # penalty
    object$S[[1]] <- object$xt$penalty
    if (ncol(object$S[[1]]) != nrow(object$S[[1]])) 
        stop("supplied penalty not square!")
    if (ncol(object$S[[1]]) != ncol(object$X)) 
        stop("supplied penalty wrong dimension!")
    if (!is.null(colnames(object$S[[1]]))) {
        a.name <- colnames(object$S[[1]])
        if (all.equal(levels(k), a.name) != TRUE) {
            stop("penalty column names don't match supplied area names!")
        }
        else {
            if (all.equal(sort(a.name), a.name) != TRUE) {
              object$S[[1]] <- object$S[[1]][levels(k), ]
              object$S[[1]] <- object$S[[1]][, levels(k)]
            }
        }
    }

    # get rank of full smoother
    if (!is.null(object $ xt $ rank)) {
      Qrank <- object $ xt $ rank
    }
    else {
        message("estimating the rank of the gmrf smooth via eigen decomp.  Supply via 'xt' to avoid this.")
        ev <- eigen(object$S[[1]], symmetric = TRUE, only.values = TRUE)$values
        Qrank <- sum(ev > .Machine$double.eps^0.8 * max(ev))
    }
    null.space.dim <- length(levels(k)) - Qrank

    # this section applies reduced rank reduction
    # note that the null space is always retained
    reduceRank <- object$bs.dim < length(levels(k))
    if (reduceRank) {
        if (object$bs.dim <= (null.space.dim + 1)) {
          object$bs.dim <- null.space.dim + 2
          message("k has been increased to 2 + null space dimension = ", object$bs.dim)
        }
        mi <- which(colSums(object$X) == 0)
        np <- ncol(object$X)
        # if there are missing observations?
        if (length(mi) > 0) {
            object$X <- rbind(matrix(0, length(mi), np), object$X)
            for (i in 1:length(mi)) object$X[i, mi[i]] <- 1
        }
        rp <- nat.param(object$X, object$S[[1]])
        ind <- (np - object$bs.dim + 1):np
        object$X <- if (length(mi)) 
            rp$X[-(1:length(mi)), ind]
        else rp$X[, ind]
        object$P <- rp$P[, ind]
        object$S[[1]] <- diag(c(rp$D[ind[ind <= rp$rank]], rep(0, sum(ind > rp$rank))))
        object$rank <- rp$rank
    } else {
        object$rank <- Qrank
    } 
    object$null.space.dim <- ncol(object$X) - object$rank
    object$knots <- k
    object$df <- ncol(object$X)
    object$plot.me <- FALSE
    #object$fixed <- FALSE # force penalty - no sense in allowing fixed to be false
    # specify the constraints
    object $ C <- object $ xt $ constraint
    # if reduced rank gmrfs are used - we need to also reduce the constraints!
    if (reduceRank & !is.null(object $ C)) {
      object $ C <- object $ C %*% t(MASS::ginv(object $ X))    
    }

    class(object) <- "gmrf.smooth"
    object
}
 

#' Predict from a gmrf smoother
#'
#' This function gives predictions from a gmrf smoother 
#' 
#'
#' @param object an object of class `"gmrf.smooth"' produced by the
#'          `smooth.construct' method.
#' @param data a list containing just the data (including any `by' variable)
#'          required by this term, with names corresponding to
#'          `object$term' (and `object$by'). The `by' variable is the
#'          last element.
#' @return A design matrix
#' @export
Predict.matrix.gmrf.smooth <- function (object, data) {
    x <- factor(data[[object$term]], levels = levels(object$knots))
    X <- model.matrix(~x - 1)
    if (!is.null(object$P)) 
        X <- X %*% object$P
    X
}



# not exported.  Copied from mgcv for now.
nat.param <- function (X, S, rank = NULL)  {
  qrx <- qr(X, tol = .Machine$double.eps^0.8)
  R <- qr.R(qrx)
  RSR <- forwardsolve(t(R), t(forwardsolve(t(R), t(S))))
  er <- eigen(RSR, symmetric = TRUE)
  if (is.null(rank) || rank < 1 || rank > ncol(S)) {
    rank <- sum(er$value > max(er$value) * .Machine$double.eps^0.8)
  }
  null.exists <- rank < ncol(X)
  D <- er$values[1:rank]
  X <- qr.Q(qrx, complete = FALSE) %*% er$vectors
  P <- backsolve(R, er$vectors)
  # scale
  ind <- 1:rank
  scale <- 1/sqrt(mean(X[, ind]^2))
  X[, ind] <- X[, ind] * scale
  P[, ind] <- P[, ind] * scale
  D <- D * scale^2
  if (null.exists) {
    ind <- (rank + 1):ncol(X)
    scalef <- 1/sqrt(mean(X[, ind]^2))
    X[, ind] <- X[, ind] * scalef
    P[, ind] <- P[, ind] * scalef
  }
  list(X = X, D = D, P = P, rank = rank, type = 0)
}