| Title: | Lightweight helper for fitting GMRF models in mgcv |
|---|---|
| Description: | A set of functions to help with the fitting of gmrf models in mgcv. The package is intended for the fitting GMRF models such as random walks and seasonal effects. The code is written as simply as possible to make is easy to see how GMRFs are constructed. |
| Authors: | Colin Millar [aut, cre] |
| Maintainer: | Colin Millar <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.0.1 |
| Built: | 2024-11-13 04:43:45 UTC |
| Source: | https://github.com/faskally/gmrf |
Details
Dar(n, phi)Dar(n, phi)
n |
the size of the GMRF |
phi |
the autoregressive parameters |
Description - This function does stuff.
a Matrix
Details
Dharm(n, wavelength = n, cyclic = FALSE)Dharm(n, wavelength = n, cyclic = FALSE)
n |
the size of the GMRF |
wavelength |
the wavelength of the harmonic. default is for this this to be n, i.e. one cycle. |
cyclic |
logical: should the smoother be cyclic, i.e. corrected for toroidal edge effects. |
Description - This function does stuff.
a Matrix
Details
Dnb(nb)Dnb(nb)
nb |
the neighbourhood structure of a regional layout |
Description - This function does stuff.
a Matrix with a column for each region and a row for each connection / edge
Details
Drw(n, order = 2, cyclic = FALSE)Drw(n, order = 2, cyclic = FALSE)
n |
the size of the GMRF. |
order |
the order of the random walk. |
cyclic |
logical: should the differences be cyclic, i.e. corrected for toroidal edge effects. |
Description - This function does stuff.
a Matrix
Details
Drw1(n, cyclic = FALSE)Drw1(n, cyclic = FALSE)
n |
the size of the GMRF. |
cyclic |
logical: should the differences be cyclic, i.e. corrected for toroidal edge effects . |
Description - This function does stuff.
a Matrix
This model treats the sum over m time points to be stationary
Dseasonal(n, m)Dseasonal(n, m)
n |
the size of the GMRF |
m |
the length of the seasonal period |
Description - This function does stuff.
a Matrix
Details. Each group of regions sums to zero.
getCnb(Q)getCnb(Q)
Q |
a precision matrix for a regional GMRF |
Description - This function does stuff.
a Matrix with a column for each region and a row for each distinct group
Details. Each group of regions sums to zero.
getFactorsnb(Q)getFactorsnb(Q)
Q |
a precision matrix for a regional GMRF |
Description - This function does stuff.
a Matrix with a column for each region and a row for each distinct group
Details
getQar(n, phi, weights = NULL)getQar(n, phi, weights = NULL)
n |
the size of the GMRF |
phi |
the autoregressive parameters |
weights |
weights to be applied to the node differences in effect allowing different variances at each time step. |
Description - This function does stuff.
a Matrix
Details
getQnb(nb, weights = NULL)getQnb(nb, weights = NULL)
nb |
the neighbourhood structure of a regional layout |
weights |
weights to be applied to the node differences in effect allowing different connections between regions to vary differently. An example of a weight could be the length of the shared border between regions |
Description - This function does stuff.
a Matrix with a column for each region and a row for each connection / edge
Details
getQpoly(poly, weights = NULL)getQpoly(poly, weights = NULL)
poly |
a spatial polygon |
weights |
weights to be applied to the node differences in effect allowing different connections between regions to vary differently. An example of a weight could be the length of the shared border between regions |
Description - This function does stuff.
a Matrix with a column for each region and a row for each connection / edge
Details
getQrw(n, order = 2, weights = NULL, cyclic = FALSE)getQrw(n, order = 2, weights = NULL, cyclic = FALSE)
n |
the size of the GMRF |
order |
the order of the random walk |
weights |
weights to be applied to the node differences (see details) |
cyclic |
logical: should the smoother be cyclic, i.e. corrected for toroidal edge effects. |
Description - This function does stuff.
what does it return
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)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)
If weights are supplied
getQrw1(n, weights = NULL, cyclic = FALSE)getQrw1(n, weights = NULL, cyclic = FALSE)
n |
the size of the GMRF |
weights |
weights to be applied to the node differences (see details) |
cyclic |
logical: should the smoother be cyclic, i.e. corrected for toroidal edge effects. |
Description - This function does stuff.
what does it return
This one is for back compatability
getRegionalGMRF(nbmat, weights = NULL)getRegionalGMRF(nbmat, weights = NULL)
nbmat |
the neighbourhood structure of a regional layout as a matrix |
weights |
weights to be applied to the node differences in effect allowing different connections between regions to vary differently. An example of a weight could be the length of the shared border between regions |
Description - This function does stuff.
a Matrix with a column for each region and a row for each connection / edge
Lightweight add-on for mgcv to allow slighly easier use of GMRF smoothing models
This function gives predictions from a gmrf smoother
## S3 method for class 'gmrf.smooth' Predict.matrix(object, data)## S3 method for class 'gmrf.smooth' Predict.matrix(object, data)
object |
an object of class ‘"gmrf.smooth"’ produced by the ‘smooth.construct’ method. |
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. |
A design matrix
Required the eigen decomposition of a precision matrix
simeQ(eQ, tol = 1e-09, rank = NULL, k = 1)simeQ(eQ, tol = 1e-09, rank = NULL, k = 1)
eQ |
an eigen decomposition of a symmetric positive semi-definate matrix corresponding to the precision matrix of an inhomogenous GMRF |
tol |
tolerance (relative to largest eigen value) for numerical lack of positive-definiteness in ‘Q’. |
rank |
the rank of the precision matrix. If this is not supplied it it is estimated by comparing the eigen values to tol * largest eigenvalue. |
k |
optional scaling of the precision matrix. This can be used to save recomputing the eigen decomposition for different smoothing parameters. |
Note if rank is suplied and is less than that true rank, the simulation will be from a reduced rank GMRF.
a single draw from with the appropriate conditional covariance structure
## create a rw2 GMRF precision matrix and simulate Q <- getQrw(100, order = 2) eQ <- eigen(Q) x <- simeQ(eQ, k = exp(5)) plot(x, type = "l")## create a rw2 GMRF precision matrix and simulate Q <- getQrw(100, order = 2) eQ <- eigen(Q) x <- simeQ(eQ, k = exp(5)) plot(x, type = "l")
Details
simQ(Q, tol = 1e-09, rank = NULL)simQ(Q, tol = 1e-09, rank = NULL)
Q |
a symmetric positive semi-definate matrix corresponding to the precision matrix of an inhomogenous GMRF |
tol |
tolerance (relative to largest eigen value) for numerical lack of positive-definiteness in ‘Q’. |
rank |
the rank of the precision matrix. If this is not supplied it it is estimated by comparing the eigen values to tol * largest eigenvalue |
This function does stuff.
a single draw from with the appropriate covariance structure
## create a rw2 GMRF precision matrix and simulate Q <- getQrw(100, order = 2) x <- simQ(Q) plot(x) ## ## simulate an AR1 GMRF with toroidal edge correction Q <- getQar(100, phi = 0.8) x <- simQ(exp(5)*Q) plot(x, type = "l")## create a rw2 GMRF precision matrix and simulate Q <- getQrw(100, order = 2) x <- simQ(Q) plot(x) ## ## simulate an AR1 GMRF with toroidal edge correction Q <- getQar(100, phi = 0.8) x <- simQ(exp(5)*Q) plot(x, type = "l")
This function is used internally in mgcv when fitting a smoother of type gmrf. More details to be added.
## S3 method for class 'gmrf.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'gmrf.smooth.spec' smooth.construct(object, data, knots)
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. |
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. |
knots |
not used. |
An object of class ‘"gmrf.smooth"’.
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)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)