#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.