| Title: | Helper for Fitting GMRF Models on River Networks in MGCV |
|---|---|
| Description: | A set of functions to help with the fitting of gmrf models for river networks in mgcv. The package is intended for the fitting 1st and 2nd order smoothers on river networks. 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-10-30 04:45:36 UTC |
| Source: | https://github.com/faskally/rc |
This increases the size of the graph by the number of new nodes
add.node(g, from, to)add.node(g, from, to)
g |
an igraph graph |
from |
a vector of node IDs |
to |
a vector of node IDs |
Description - This function does stuff.
an igraph graph
Details
buildTopo(lines, from = "FNODE_", to = "TNODE_")buildTopo(lines, from = "FNODE_", to = "TNODE_")
lines |
a spatial lines data.frame containing a river network |
from |
the column name containing the 'from' node IDs |
to |
the column name containing the 'to' node IDs |
Description - This function does stuff.
what does it return
Returns the incidence matrix of a graph
Dgraph(g)Dgraph(g)
g |
an igraph graph specifying the dependence structure |
Description - This function does stuff.
a Matrix with a column for each node and a row for each connection / edge
Returns the laplacian matrix of a graph. If weights are supplied then the weighted lacplacian is returned. Sensible weights are based on the flows of the incoming upstream graph edges.
getQgraph(g, weights = NULL)getQgraph(g, weights = NULL)
g |
an igraph graph specifying the dependence structure |
weights |
weights to be applied to the edges of the graph (see details) |
Weights argument is passed onto the graph.laplacian function from the igraph package. If weights is NULL (default) and the graph has an edge attribute called weight, then it will be used automatically. Set this to NA if you want the unweighted Laplacian on a graph that has a weight edge attribute.
matrix
# make a simple river netork M <- rbind(c( 1, -1, 0, 0, 0, 0, 0, 0), c( 0, 1, -1, 0, 0, 0, 0, 0), c( 0, 1, 0, -1, 0, 0, 0, 0), c( 0, 0, 1, 0, -1, 0, 0, 0), c( 0, 0, 1, 0, 0, -1, 0, 0), c( 0, 0, 0, 1, 0, 0, -1, 0), c( 0, 0, 0, 1, 0, 0, 0, -1)) A <- -1 * t(M) %*% M diag(A) <- 0 g <- graph.adjacency(A, mode = "undirected") # plot graph plot(g) # add a node in between all other nodes g <- add.node(g, c(1,2,2,3,3,4,4), c(2,3,4,5,6,7,8)) # get precision matrix for river network Q <- getQgraph(g) # simulate river network effect x <- simQ(Q) # simulate observations, 1 per region in this case y <- x + rnorm(length(x))*0.5 # get node ids for observations # note node ID should not be character # it can either be numeric, or a factor. nid <- factor(rownames(Q)) # collect data in a list dat <- data.frame(y = y, nid = nid) # provide rank to avoid calculating this at every fit xtr <- list(penalty = Q, rank = nrow(Q) - 1) # plot simulation breaks <- seq(min(y)-0.001, max(y)+0.001, length = 11) par(mfrow = c(1,2)) plot(g, vertex.color = heat.colors(length(breaks)-1)[as.numeric(cut(x, breaks))]) # fit a model g1 <- gam(y ~ s(nid, bs = "gmrf", xt = xtr), method="REML", data = dat) summary(g1) # plot fitted values plot(g, vertex.color = heat.colors(length(breaks)-1)[as.numeric(cut(fitted(g1), breaks))])# make a simple river netork M <- rbind(c( 1, -1, 0, 0, 0, 0, 0, 0), c( 0, 1, -1, 0, 0, 0, 0, 0), c( 0, 1, 0, -1, 0, 0, 0, 0), c( 0, 0, 1, 0, -1, 0, 0, 0), c( 0, 0, 1, 0, 0, -1, 0, 0), c( 0, 0, 0, 1, 0, 0, -1, 0), c( 0, 0, 0, 1, 0, 0, 0, -1)) A <- -1 * t(M) %*% M diag(A) <- 0 g <- graph.adjacency(A, mode = "undirected") # plot graph plot(g) # add a node in between all other nodes g <- add.node(g, c(1,2,2,3,3,4,4), c(2,3,4,5,6,7,8)) # get precision matrix for river network Q <- getQgraph(g) # simulate river network effect x <- simQ(Q) # simulate observations, 1 per region in this case y <- x + rnorm(length(x))*0.5 # get node ids for observations # note node ID should not be character # it can either be numeric, or a factor. nid <- factor(rownames(Q)) # collect data in a list dat <- data.frame(y = y, nid = nid) # provide rank to avoid calculating this at every fit xtr <- list(penalty = Q, rank = nrow(Q) - 1) # plot simulation breaks <- seq(min(y)-0.001, max(y)+0.001, length = 11) par(mfrow = c(1,2)) plot(g, vertex.color = heat.colors(length(breaks)-1)[as.numeric(cut(x, breaks))]) # fit a model g1 <- gam(y ~ s(nid, bs = "gmrf", xt = xtr), method="REML", data = dat) summary(g1) # plot fitted values plot(g, vertex.color = heat.colors(length(breaks)-1)[as.numeric(cut(fitted(g1), breaks))])
Returns the laplacian matrix of a graph. If weights are supplied then the weighted lacplacian is returned. Sensible weights are based on the flows of the incoming upstream graph edges.
getWRW1Mat(g, weights = NULL)getWRW1Mat(g, weights = NULL)
g |
an igraph graph specifying the dependence structure |
weights |
weights to be applied to the edges of the graph (see details) |
Weights argument is passed onto the graph.laplacian function from the igraph package. If weights is NULL (default) and the graph has an edge attribute called weight, then it will be used automatically. Set this to NA if you want the unweighted Laplacian on a graph that has a weight edge attribute.
what does it return
Lightweight add-on for mgcv to allow slighly easier use of GMRF smoothing models