| Title: | Bayesian Spatio-Temporal Analysis in Stream Networks |
|---|---|
| Description: | Fits Bayesian spatio-temporal models and makes predictions on stream networks using the approach by Santos-Fernandez, Edgar, et al. (2022)."Bayesian spatio-temporal models for stream networks" and Santos-Fernandez, Edgar, et al. (2023). "SSNbayes: An R Package for Bayesian Spatio-Temporal Modelling on Stream Networks". In these models, spatial dependence is captured using stream distance and flow connectivity, while temporal autocorrelation is modelled using vector autoregression methods. |
| Authors: | Edgar Santos-Fernandez [aut, cre, cph] (ORCID: <https://orcid.org/0000-0001-5962-5417>) |
| Maintainer: | Edgar Santos-Fernandez <[email protected]> |
| License: | GPL-2 |
| Version: | 0.1.2 |
| Built: | 2026-10-24 06:22:30 UTC |
| Source: | https://github.com/edgarsantos-fernandez/ssnbayes |
Fits a Bayesian spatio-temporal model with autoregressive temporal structure and stream-network-based spatial modeling using Stan. The method automatically selects spatial knots based on iterative RMSE performance and returns predictions with uncertainty and anomaly indicators.
barst( formula = yobs ~ X1 + X2 + X3, data, path, time_method = list("ar", "date"), space_method = list("use_ssn", c("Exponential.taildown")), iter = 3000, warmup = 1500, chains = 3, refresh = 100, loglik = FALSE, ppd = TRUE, seed = 123, total_iter = 10 )barst( formula = yobs ~ X1 + X2 + X3, data, path, time_method = list("ar", "date"), space_method = list("use_ssn", c("Exponential.taildown")), iter = 3000, warmup = 1500, chains = 3, refresh = 100, loglik = FALSE, ppd = TRUE, seed = 123, total_iter = 10 )
formula |
A formula specifying the response and covariates, e.g., |
data |
A data frame with the full space-time dataset, including covariates and response. |
path |
Path to the folder containing SSN or spatial network objects and distance matrices. |
time_method |
A list specifying temporal structure, e.g., |
space_method |
A list specifying spatial covariance function and structure, e.g., |
iter |
Number of iterations per Stan run (default 3000). |
warmup |
Number of warmup iterations (default 1500). |
chains |
Number of Stan chains (default 3). |
refresh |
Frequency of Stan progress printout (default 100). |
loglik |
Logical, whether to return pointwise log-likelihood (default |
ppd |
Logical, whether to return posterior predictive distributions (default |
seed |
Random seed (default: 123). |
total_iter |
Number of RMSE-based refinement iterations (default: 10). |
A list with elements:
fitThe final Stan model object.
predData frame with posterior predictions, intervals, and anomaly indicators.
data_listList of Stan inputs used for the final model.
replicMatrix of knot selections across iterations.
Example stream-network data from the Boise River Basin case study in
Santos-Fernandez et al. (2022). The example has two parts, both stored
under inst/extdata:
a spatial stream network directory, boise.ssn;
an observation data frame, boise_obs.RDS.
boise_obs.RDS is the observation data used in
Santos-Fernandez et al. (2022), Computational Statistics & Data
Analysis. boise.ssn is the spatial stream network supplied
with those observations. It contains network 3, the 42 sites used in
that paper, and other river networks as well.
boise_obs.RDS is a data frame with 3,654 rows and 49 columns.
Columns used in the published analysis, or identical to a defined site attribute:
Date. Observation date. Eighty-seven dates from
2010-12-01 to 2015-11-11. Each of the 42 locID values
occurs once on every date.
Numeric location identifier. 42 locations.
Numeric daily mean stream temperature in degrees Celsius.
Present on all 2,620 rows with dataset == "train" and
missing on all 1,034 rows with dataset == "test".
Character. "train" (2,620 rows) or
"test" (1,034 rows).
Numeric daily air temperature in degrees Celsius, the air-temperature covariate in the case study.
Numeric stream-slope covariate used in the case study.
It differs from SLOPE by at most 0.0005.
Integer elevation in meters. Identical to Elev_M.
Numeric cumulative drainage area in square kilometers.
Identical to CUMDRAINAG.
Numeric first Fourier terms. Within each date,
sin takes one value and cos takes one value.
Numeric latitude and longitude in decimal degrees.
Numeric. Identical to S1_93_11.
Numeric. Differs from MA_HistCMS by at most
0.00005.
Integer unique point identifier.
Factor network identifier. Every row is network 3. The factor also has unused levels 1, 2, and 5.
Character. Every row is "Boise_river".
Character stream name (30 streams).
Character stream name, as defined for the sites
shapefile.
Integer elevation in meters.
Numeric site identifier.
Numeric stream-reach slope, in percent.
Numeric cumulative drainage area in square kilometers.
Numeric mean August air temperature in degrees Celsius. One value is stored in this table (19.15). The workbook describes it as a basin-wide August mean.
Numeric mean August stream discharge. One value is stored in this table (4.15).
Numeric historical composite: 19-year average August mean stream temperature for 1993-2011.
Numeric mean annual flow from the historic NorWeST model, in cubic meters per second.
Numeric latitude and longitude.
Integer. FID of the line segment the site was snapped to.
Numeric distance the site was moved when it was snapped to the network.
Numeric coordinates of the snapped site.
Numeric angle from the site to the nearest point on the line segment.
Integer. Edge the site was snapped to.
Numeric position of the site along that edge.
Numeric upstream distance from the network outlet to the site.
Numeric additive function value.
The variable-definitions workbook in boise.ssn is the source
for the sites attribute descriptions above.
These columns are also present. Their names are listed without a
further definition: UTM_Xcoord, UTM_Ycoord,
locID_old, locID_con, LOCID, temp2,
locID_backup, id, temp_backup,
temp_pred_lm, and temp_pred_lm2.
locID_backup is not the same vector as locID.
temp_pred_lm and temp_pred_lm2 are missing on every
training row and present on every test row.
boise.ssn is the directory passed to SSN2::ssn_import().
It contains:
edges: 32,852 stream segments (LINESTRING);
sites: 226 observed sites (POINT);
preds: 50,694 prediction locations (POINT);
netID1.dat through netID8.dat;
binaryID.db;
Variable_Definitions.xlsx, a workbook of shapefile
attribute definitions for sites, preds, and
edges.
The prediction-site dataset is named preds. Pass
predpts = "preds" to import it. binaryID.db is already
present, so the default overwrite = FALSE is used and that file
is left unchanged.
The stream features use the projected coordinate reference system stored in the shapefiles (Albers, GNLCC, meters).
boise_obs.RDS is the data frame used in the Computational
Statistics \& Data Analysis (CSDA) paper: daily
records for the 42 locations on network 3 (netID 3), on each of
87 dates (3,654 rows). Every row in this file has netID 3.
Santos-Fernandez et al. (2022) describe these 87 dates as a systematic
21-day subsample of mean daily temperatures measured from 2010-12-01 to
2015-12-01. In this file the dates run from 2010-12-01 to 2015-11-11.
boise.ssn contains other river networks in addition to network 3.
In sites, network 3 has the 42 sites used in the paper.
Networks 1, 2, and 5 contain the other observed sites (68, 105, and 11).
The 6,422 prediction locations in the paper are the preds
features with netID 3. That shapefile has 50,694 features in
total because it also includes prediction locations on other networks.
The response in the case study is daily mean stream temperature
(temp, degrees Celsius). The published mean function uses stream
slope (slope), elevation (elev), cumulative drainage area
(drainage), daily air temperature (air_temp), and the
first Fourier pair (sin, cos).
Santos-Fernandez, E., Ver Hoef, J. M., Peterson, E. E., McGree, J., Isaak, D. J., and Mengersen, K. (2022). Bayesian spatio-temporal models for stream networks. Computational Statistics & Data Analysis, 170, 107446. doi:10.1016/j.csda.2022.107446
library(SSNbayes) # Path to the Boise stream-network directory boise_path <- system.file( "extdata", "boise.ssn", package = "SSNbayes" ) # Import the stream network and the preds prediction sites boise <- SSN2::ssn_import( boise_path, predpts = "preds" ) # Load the observation data boise_obs <- readRDS( system.file( "extdata", "boise_obs.RDS", package = "SSNbayes" ) ) head(boise_obs)library(SSNbayes) # Path to the Boise stream-network directory boise_path <- system.file( "extdata", "boise.ssn", package = "SSNbayes" ) # Import the stream network and the preds prediction sites boise <- SSN2::ssn_import( boise_path, predpts = "preds" ) # Load the observation data boise_obs <- readRDS( system.file( "extdata", "boise_obs.RDS", package = "SSNbayes" ) ) head(boise_obs)
Collapses a SpatialStreamNetwork object into a data frame
collapse(ssn, par = "afvArea")collapse(ssn, par = "afvArea")
ssn |
An S4 SpatialStreamNetwork object created with SSN2 package. |
par |
A spatial parameter such as the computed_afv (additive function value). |
The parameters (par) has to be present in the observed data frame via ssn_get_data(n, name = "obs"). More details of the argument par can be found in the additive.function() from SSN .
A data frame with the lat and long of the line segments in the network. The column line_id refers to the ID of the line.
#require("SSN2") #path <- system.file("extdata/clearwater.ssn", package = "SSNbayes") #ssn <- SSN2::ssn_import(path, predpts = "preds", overwrite = TRUE) #t.df <- collapse(ssn, par = 'afvArea')#require("SSN2") #path <- system.file("extdata/clearwater.ssn", package = "SSNbayes") #ssn <- SSN2::ssn_import(path, predpts = "preds", overwrite = TRUE) #t.df <- collapse(ssn, par = 'afvArea')
Combines two 'sf' objects, the sensor locations and river network, in to an 'sfnetwork'.
convert_network_to_graph(river_network, sensor_locations)convert_network_to_graph(river_network, sensor_locations)
river_network |
An 'sf' object of 'MULTILINESTRING'(s) containing the streams and rivers, cropped to the sensor locations. |
sensor_locations |
An 'sf' dataframe containing sensor locations and 'locID'. |
This function "snaps" sensor locations to the nearest point on the network. Additionally, the sensor locations are stored in the graph, 'g', and can be accessed via 'V(g)$type', as "sensor".
An 'sfnetwork' graph containing the river network and sensors.
Sean Francis
Creates a list containing the stream distances and weights
dist_weight_mat(path = path, net = 1, addfunccol = "addfunccol")dist_weight_mat(path = path, net = 1, addfunccol = "addfunccol")
path |
Path to the files |
net |
(optional) A network from the SSN2 object |
addfunccol |
(optional) A parameter to compute the spatial weights |
A list of matrices
path <- system.file("extdata/clearwater.ssn", package = "SSNbayes") mat_all <- dist_weight_mat(path, net = 2, addfunccol='afvArea')path <- system.file("extdata/clearwater.ssn", package = "SSNbayes") mat_all <- dist_weight_mat(path, net = 2, addfunccol='afvArea')
The output matrices are symmetric except the hydrologic distance matrix D.
dist_weight_mat_preds(path = path, net = 1, addfunccol = "addfunccol")dist_weight_mat_preds(path = path, net = 1, addfunccol = "addfunccol")
path |
Path with the name of the SpatialStreamNetwork object |
net |
(optional) A network from the SpatialStreamNetwork object |
addfunccol |
(optional) A parameter to compute the spatial weights |
A list of matrices
## Not run: path <- system.file("extdata/clearwater.ssn", package = "SSNbayes") mat_all_pred <- dist_weight_mat_preds(path, net = 2, addfunccol='afvArea') ## End(Not run)## Not run: path <- system.file("extdata/clearwater.ssn", package = "SSNbayes") mat_all_pred <- dist_weight_mat_preds(path, net = 2, addfunccol='afvArea') ## End(Not run)
Converts a dataframe containing sensor data in to an 'sf' object
format_sensor_data(sensor_data, lon_name, lat_name)format_sensor_data(sensor_data, lon_name, lat_name)
sensor_data |
A dataframe containing sensor locations and a column called 'locID' denoting the location IDs. |
lon_name |
The column name containing the longitude of the sensor. |
lat_name |
The column name containing the latitude of the sensor. |
An 'sf' dataframe containing sensor locations and 'locID'.
Sean Francis
Creates a list containing the stream distances and weights
gen_distance_matrices(network, sensor_locations, root_loc)gen_distance_matrices(network, sensor_locations, root_loc)
network |
An 'sfnetwork' containing the river network and sensor locations. |
sensor_locations |
An 'sf' dataframe containing sensor locations and 'locID'. |
root_loc |
A sensor location which all water flows towards. |
A list of matrices.
Sean Francis
It will take a dataframe containing sensor locations, with a coulmn 'locID' denoting location IDs (can be string or numeric). Requires a root sensor location, in which all sensor locations are connected to, and flow towards.
generate_osm_ssn( sensor_data, lon_name, lat_name, root_loc, plot_network = TRUE, gen_pred_sites = FALSE, num_pred_sites = NA )generate_osm_ssn( sensor_data, lon_name, lat_name, root_loc, plot_network = TRUE, gen_pred_sites = FALSE, num_pred_sites = NA )
sensor_data |
A dataframe containing sensor locations and a column called 'locID' denoting the location IDs. |
lon_name |
The column name containing the longitude of the sensor. |
lat_name |
The column name containing the latitude of the sensor. |
root_loc |
A sensor location which all water flows towards. |
plot_network |
A bool indicating whether to plot the network. Useful for determining potential errors. Will be saved in object. |
gen_pred_sites |
A bool indicating whether to generate new prediction sites. |
num_pred_sites |
Required if 'gen_pred_sites=TRUE'. An integer indicating the number of prediction sites to generate. |
Assumes longitude and latitude data are in the 4326 coordinate reference system. If they are not, data will need to be converted beforehand using 'st_transform()'. Predictions on sites using gen_pred_sites is under development specially for tail-down models.
An 'osm_ssn' object, containing: the graph network as an 'sfnetwork'; the distance adjacency matrices; if 'plot_network=TRUE', an interactive leaflet plot, if 'gen_pred_sites=TRUE', an 'sf' object containing the locations of the predictions sites.
Sean Francis
# require('SSNdata') # clear <- readRDS(system.file("extdata/clear_obs.RDS", package = "SSNdata")) # # # Generate osm ssn # clear_osm_ssn <- generate_osm_ssn(clear, # "long", "lat", # root_loc = 12, # plot_network = TRUE, # gen_pred_sites = TRUE, # num_pred_sites = 20) # # Show plot # clear_osm_ssn$plot# require('SSNdata') # clear <- readRDS(system.file("extdata/clear_obs.RDS", package = "SSNdata")) # # # Generate osm ssn # clear_osm_ssn <- generate_osm_ssn(clear, # "long", "lat", # root_loc = 12, # plot_network = TRUE, # gen_pred_sites = TRUE, # num_pred_sites = 20) # # Show plot # clear_osm_ssn$plot
Creates an 'sf' dataframe containing equally spaced points on a river network
generate_prediction_locations(river_network, num_pred_sites)generate_prediction_locations(river_network, num_pred_sites)
river_network |
An 'sf' object of 'MULTILINESTRING'(s) containing the streams and rivers, cropped to the sensor locations. |
num_pred_sites |
An integer indicating the number of prediction sites to generate. |
An 'sf' dataframe of 'POINT'(s) containing prediction sites.
Sean Francis
Takes an 'sf' object as input, and queries Open Street Maps for rivers and streams within the bounds of the sensor locations.
generate_river_network(sensor_locations, root_loc)generate_river_network(sensor_locations, root_loc)
sensor_locations |
An 'sf' dataframe containing sensor locations and 'locID'. |
root_loc |
A sensor location which all water flows towards. |
This function queries Open Street Maps using 'osmdata' package in R. Crops the data to the bounding box of the sensor locations. Removes all lines not connected to the outlet/root_loc.
An 'sf' object of 'MULTILINESTRING'(s) containing the streams and rivers, cropped to the sensor location, with all lines not containing the outlet/root_loc removed.
Sean Francis
Use predict.ssnbayes() instead. It will take an observed and a prediction data frame. It requires the same number of observation/locations per day. It requires location id (locID) and points id (pid). The locID are unique for each site. The pid is unique for each observation. Missing values are allowed in the response but not in the covariates.
krig( object = object, mat_all_preds = mat_all_preds, nsamples = 10, start = 1, chunk_size = 50, obs_data = obs_data, pred_data = pred_data, net = net, seed = seed )krig( object = object, mat_all_preds = mat_all_preds, nsamples = 10, start = 1, chunk_size = 50, obs_data = obs_data, pred_data = pred_data, net = net, seed = seed )
object |
A stanfit object returned from ssnbayes |
mat_all_preds |
A list with the distance/weights matrices |
nsamples |
The number of samples to draw from the posterior distributions. (nsamples <= iter) |
start |
(optional) The starting location id |
chunk_size |
(optional) the number of locID to make prediction from |
obs_data |
The observed data frame |
pred_data |
The predicted data frame |
net |
(optional) Network from the SSN object |
seed |
(optional) A seed for reproducibility |
A data frame
Edgar Santos-Fernandez
Use predict.ssnbayes() instead. It will take an observed and a prediction data frame. It requires the same number of observation/locations per day. It requires location id (locID) and points id (pid). The locID are unique for each site. The pid is unique for each observation. Missing values are allowed in the response but not in the covariates.
krig2( object = object, mat_all_preds = mat_all_preds, nsamples = 10, start = 1, chunk_size = 50, obs_data = obs_data, pred_data = pred_data, net = net, seed = seed )krig2( object = object, mat_all_preds = mat_all_preds, nsamples = 10, start = 1, chunk_size = 50, obs_data = obs_data, pred_data = pred_data, net = net, seed = seed )
object |
A stanfit object returned from ssnbayes |
mat_all_preds |
A list with the distance/weights matrices |
nsamples |
The number of samples to draw from the posterior distributions. (nsamples <= iter) |
start |
(optional) The starting location id |
chunk_size |
(optional) the number of locID to make prediction from |
obs_data |
The observed data frame |
pred_data |
The predicted data frame |
net |
(optional) Network from the SSN object |
seed |
(optional) A seed for reproducibility |
A data frame
Edgar Santos-Fernandez
A simple modeling function using a formula and data
mylm(formula, data)mylm(formula, data)
formula |
A formula as in lm() |
data |
A data.frame containing the elements specified in the formula |
A list of matrices
Jay ver Hoef
options(na.action='na.pass') data("iris") out_list = mylm(formula = Petal.Length ~ Sepal.Length + Sepal.Width, data = iris)options(na.action='na.pass') data("iris") out_list = mylm(formula = Petal.Length ~ Sepal.Length + Sepal.Width, data = iris)
Use predict.ssnbayes() instead. It will take an observed and a prediction data frame. It requires the same number of observation/locations per day. It requires location id (locID) and points id (pid). The locID are unique for each site. The pid is unique for each observation. Missing values are allowed in the response but not in the covariates.
pred_ssnbayes( object = object, use_osm_ssn = TRUE, osm_ssn = osm_ssn, path = path, obs_data = obs_data, pred_data = pred_data, net = 1, nsamples = 100, addfunccol = "afvArea", locID_pred = locID_pred, chunk_size = chunk_size, seed = seed )pred_ssnbayes( object = object, use_osm_ssn = TRUE, osm_ssn = osm_ssn, path = path, obs_data = obs_data, pred_data = pred_data, net = 1, nsamples = 100, addfunccol = "afvArea", locID_pred = locID_pred, chunk_size = chunk_size, seed = seed )
object |
A stanfit object returned from ssnbayes |
use_osm_ssn |
Use a supplied osm_ssn instead a SpatialStreamNetwork object. |
osm_ssn |
The osm_ssn to be used. |
path |
If not using an osm_ssn, path with the name of the SpatialStreamNetwork object. |
obs_data |
The observed data frame |
pred_data |
The predicted data frame |
net |
(optional) Network from the SSN object |
nsamples |
The number of samples to draw from the posterior distributions. (nsamples <= iter) |
addfunccol |
If not using an osm_ssn, the variable used for spatial weights. |
locID_pred |
(optional) the location id for the predictions. Used when the number of pred locations is large. |
chunk_size |
(optional) the number of locID to make prediction from |
seed |
(optional) A seed for reproducibility |
A data frame
Edgar Santos-Fernandez
#pred <- pred_ssnbayes(path = path, #obs_data = clear, #stanfit = fit_ar, #pred_data = preds, #net = 2, #nsamples = 100, # number of samples to use from the posterior in the stanfit object #addfunccol = 'afvArea') # variable used for spatial weights#pred <- pred_ssnbayes(path = path, #obs_data = clear, #stanfit = fit_ar, #pred_data = preds, #net = 2, #nsamples = 100, # number of samples to use from the posterior in the stanfit object #addfunccol = 'afvArea') # variable used for spatial weights
It will take an observed and a prediction data frame. It requires the same number of observation/locations per day. It requires location id (locID) and points id (pid). The locID are unique for each site. The pid is unique for each observation. Missing values are allowed in the response but not in the covariates.
## S3 method for class 'ssnbayes' predict( object = object, ..., path = path, obs_data = obs_data, pred_data = pred_data, net = net, nsamples = nsamples, addfunccol = addfunccol, locID_pred = locID_pred, chunk_size = chunk_size, seed = seed )## S3 method for class 'ssnbayes' predict( object = object, ..., path = path, obs_data = obs_data, pred_data = pred_data, net = net, nsamples = nsamples, addfunccol = addfunccol, locID_pred = locID_pred, chunk_size = chunk_size, seed = seed )
object |
A stanfit object returned from ssnbayes |
... |
Other parameters |
path |
Path with the name of the SpatialStreamNetwork object |
obs_data |
The observed data frame |
pred_data |
The predicted data frame |
net |
(optional) Network from the SSN object |
nsamples |
The number of samples to draw from the posterior distributions. (nsamples <= iter) |
addfunccol |
The variable used for spatial weights |
locID_pred |
(optional) the location id for the predictions. Used when the number of pred locations is large. |
chunk_size |
(optional) the number of locID to make prediction from |
seed |
(optional) A seed for reproducibility |
The returned data frame is melted to produce a long dataset. See examples. Currently, the predict() function produces predictions for normal random variables. However, this can be easily transformed in to counts (Poisson distributed) and presence/absence (binomial distributed).
A data frame with the location (locID), time point (date), plus the MCMC draws from the posterior from 1 to the number of iterations. The locID0 column is an internal consecutive location ID (locID) produced in the predictions, starting at max(locID(observed data)) + 1. It is used internally in the way predictions are made in chunks.
Edgar Santos-Fernandez
#require('SSNdata') #clear_preds <- readRDS(system.file("extdata/clear_preds.RDS", package = "SSNdata")) #clear_preds$y <- NA #pred <- predict(object = fit_ar, # path = path, # obs_data = clear, # pred_data = clear_preds, # net = 2, # nsamples = 100, # numb of samples from the posterior # addfunccol = 'afvArea', # var for spatial weights # locID_pred = locID_pred, # chunk_size = 60)#require('SSNdata') #clear_preds <- readRDS(system.file("extdata/clear_preds.RDS", package = "SSNdata")) #clear_preds$y <- NA #pred <- predict(object = fit_ar, # path = path, # obs_data = clear, # pred_data = clear_preds, # net = 2, # nsamples = 100, # numb of samples from the posterior # addfunccol = 'afvArea', # var for spatial weights # locID_pred = locID_pred, # chunk_size = 60)
It will take an observed and a prediction data frame. It requires the same number of observation/locations per day. It requires location id (locID) and points id (pid). The locID are unique for each site. The pid is unique for each observation. Missing values are allowed in the response but not in the covariates.
## S3 method for class 'ssnbayes2' predict( object = object, ..., use_osm_ssn = TRUE, osm_ssn = osm_ssn, path = path, obs_data = obs_data, pred_data = pred_data, net = net, nsamples = nsamples, addfunccol = addfunccol, locID_pred = locID_pred, chunk_size = chunk_size, seed = seed )## S3 method for class 'ssnbayes2' predict( object = object, ..., use_osm_ssn = TRUE, osm_ssn = osm_ssn, path = path, obs_data = obs_data, pred_data = pred_data, net = net, nsamples = nsamples, addfunccol = addfunccol, locID_pred = locID_pred, chunk_size = chunk_size, seed = seed )
object |
A stanfit object returned from ssnbayes |
... |
Other parameters |
use_osm_ssn |
Use a supplied osm_ssn instead a SpatialStreamNetwork object. |
osm_ssn |
The osm_ssn to be used. |
path |
If not using an osm_ssn, path with the name of the SpatialStreamNetwork object. |
obs_data |
The observed data frame |
pred_data |
The predicted data frame |
net |
(optional) Network from the SSN object |
nsamples |
The number of samples to draw from the posterior distributions. (nsamples <= iter) |
addfunccol |
If not using an osm_ssn, the variable used for spatial weights |
locID_pred |
(optional) the location id for the predictions. Used when the number of pred locations is large. |
chunk_size |
(optional) the number of locID to make prediction from |
seed |
(optional) A seed for reproducibility |
The returned data frame is melted to produce a long dataset. See examples.
A data frame with the location (locID), time point (date), plus the MCMC draws from the posterior from 1 to the number of iterations. The locID0 column is an internal consecutive location ID (locID) produced in the predictions, starting at max(locID(observed data)) + 1. It is used internally in the way predictions are made in chunks.
Edgar Santos-Fernandez
# require('SSNdata') # require('sf') # # clear <- readRDS(system.file("extdata/clear_obs.RDS", package = "SSNdata")) # # clear_osm_ssn <- generate_osm_ssn(clear, "long", "lat", root_loc = 12, plot_network = TRUE) # # formula = y ~ SLOPE + elev + air_temp + sin + cos # # family = "gaussian" # # # Note - missing data must be imputed before using the predict() function. # # This can be done using: # data_impute <- mtsdi::mnimput( # formula = temp ~ SLOPE + elev + h2o_area + air_temp + sin + cos, # dataset = clear, # eps = 1e-3, # ts = FALSE, # method = "glm" # )$filled.dataset # # clear <- left_join(clear, data_impute, by = c("elev", "air_temp", "sin", "cos", "SLOPE")) # # clear$y <- clear$temp.y # # fit_ar <- ssnbayes2(formula = formula, # data = clear, # osm_ssn = clear_osm_ssn, # family = family, # time_method = list("ar", "date"), # space_method = list('use_osm_ssn', c("Exponential.taildown")), # iter = 2000, # warmup = 1000, # chains = 3, # cores = 3) # # # Get the coordinate reference system of the near_X and near_Y variables # data_crs <- system.file("extdata/clearwater.ssn", package = "SSNbayes") %>% # SSN2::ssn_import(predpts = "preds", overwrite = TRUE) %>% # SSN2::ssn_get_data(name = "preds") %>% # st_crs # # # Load in the predictions, and convert to 4326 CRS # clear_preds <- readRDS(system.file("extdata/clear_preds.RDS", package = "SSNbayes")) %>% # st_as_sf(coords = c("NEAR_X","NEAR_Y"), # crs = data_crs) %>% # st_transform(crs = 4326) # # # Extract the long and lat columns # xy <- st_coordinates(clear_preds) # # colnames(xy) <- c("long", "lat") # # Merge back in with the prediction dataset # clear_preds <- cbind(clear_preds, xy) %>% # data.frame() %>% # select(-geometry) # # same_names <- intersect(names(clear), names(clear_preds)) # # # Combine all observed and prediction sites into 1 dataframe # all_sites <- rbind(clear %>% select(same_names), # data.frame(clear_preds) %>% select(same_names) # ) # # # clear_preds_osm_ssn <- generate_osm_ssn(sensor_data = all_sites, # lon_name = "long", lat_name = "lat", # root_loc = 12, # plot_network = TRUE, # gen_pred_sites = FALSE) # # # # locs <- clear_preds_osm_ssn$dist_mat_all$e %>% colnames # # # clear_krig <- clear %>% # filter(locID %in% locs) # # clear_krig_preds <- clear_preds %>% # filter(locID %in% locs) # # # preds <- predict(object = fit_ar, # use_osm_ssn = TRUE, # osm_ssn = clear_preds_osm_ssn, # obs_data = clear_krig, # pred_data = clear_krig_preds, # seed = seed, # nsamples = 25, # chunk_size = length(unique(clear_krig_preds$locID)) # ) # # # Condense data to posterior point estimates # ys <- reshape2::melt(preds, id.vars = c('locID0', 'locID', 'date'), value.name ='y') # ys$iter <- gsub("[^0-9.-]", "", ys$variable) # ys$variable <- NULL # # ys <- data.frame(ys) %>% dplyr::group_by(date, locID, locID0) %>% # dplyr::summarise("sd" = sd(y, na.rm=T), # "y_pred" = mean(y, na.rm=T)) # # ys <- dplyr::arrange(ys, locID) #'# require('SSNdata') # require('sf') # # clear <- readRDS(system.file("extdata/clear_obs.RDS", package = "SSNdata")) # # clear_osm_ssn <- generate_osm_ssn(clear, "long", "lat", root_loc = 12, plot_network = TRUE) # # formula = y ~ SLOPE + elev + air_temp + sin + cos # # family = "gaussian" # # # Note - missing data must be imputed before using the predict() function. # # This can be done using: # data_impute <- mtsdi::mnimput( # formula = temp ~ SLOPE + elev + h2o_area + air_temp + sin + cos, # dataset = clear, # eps = 1e-3, # ts = FALSE, # method = "glm" # )$filled.dataset # # clear <- left_join(clear, data_impute, by = c("elev", "air_temp", "sin", "cos", "SLOPE")) # # clear$y <- clear$temp.y # # fit_ar <- ssnbayes2(formula = formula, # data = clear, # osm_ssn = clear_osm_ssn, # family = family, # time_method = list("ar", "date"), # space_method = list('use_osm_ssn', c("Exponential.taildown")), # iter = 2000, # warmup = 1000, # chains = 3, # cores = 3) # # # Get the coordinate reference system of the near_X and near_Y variables # data_crs <- system.file("extdata/clearwater.ssn", package = "SSNbayes") %>% # SSN2::ssn_import(predpts = "preds", overwrite = TRUE) %>% # SSN2::ssn_get_data(name = "preds") %>% # st_crs # # # Load in the predictions, and convert to 4326 CRS # clear_preds <- readRDS(system.file("extdata/clear_preds.RDS", package = "SSNbayes")) %>% # st_as_sf(coords = c("NEAR_X","NEAR_Y"), # crs = data_crs) %>% # st_transform(crs = 4326) # # # Extract the long and lat columns # xy <- st_coordinates(clear_preds) # # colnames(xy) <- c("long", "lat") # # Merge back in with the prediction dataset # clear_preds <- cbind(clear_preds, xy) %>% # data.frame() %>% # select(-geometry) # # same_names <- intersect(names(clear), names(clear_preds)) # # # Combine all observed and prediction sites into 1 dataframe # all_sites <- rbind(clear %>% select(same_names), # data.frame(clear_preds) %>% select(same_names) # ) # # # clear_preds_osm_ssn <- generate_osm_ssn(sensor_data = all_sites, # lon_name = "long", lat_name = "lat", # root_loc = 12, # plot_network = TRUE, # gen_pred_sites = FALSE) # # # # locs <- clear_preds_osm_ssn$dist_mat_all$e %>% colnames # # # clear_krig <- clear %>% # filter(locID %in% locs) # # clear_krig_preds <- clear_preds %>% # filter(locID %in% locs) # # # preds <- predict(object = fit_ar, # use_osm_ssn = TRUE, # osm_ssn = clear_preds_osm_ssn, # obs_data = clear_krig, # pred_data = clear_krig_preds, # seed = seed, # nsamples = 25, # chunk_size = length(unique(clear_krig_preds$locID)) # ) # # # Condense data to posterior point estimates # ys <- reshape2::melt(preds, id.vars = c('locID0', 'locID', 'date'), value.name ='y') # ys$iter <- gsub("[^0-9.-]", "", ys$variable) # ys$variable <- NULL # # ys <- data.frame(ys) %>% dplyr::group_by(date, locID, locID0) %>% # dplyr::summarise("sd" = sd(y, na.rm=T), # "y_pred" = mean(y, na.rm=T)) # # ys <- dplyr::arrange(ys, locID) #'
It requires the same number of observation/locations per day. It requires location id (locID) and points id (pid). The locID are unique for each site. The pid is unique for each observation. Missing values are allowed in the response but not in the covariates.
ssnbayes( formula = formula, data = data, path = path, time_method = time_method, space_method = space_method, iter = 3000, warmup = 1500, chains = 3, refresh = max(iter/100, 1), net = 1, addfunccol = addfunccol, loglik = FALSE, ppd = FALSE, seed = seed )ssnbayes( formula = formula, data = data, path = path, time_method = time_method, space_method = space_method, iter = 3000, warmup = 1500, chains = 3, refresh = max(iter/100, 1), net = 1, addfunccol = addfunccol, loglik = FALSE, ppd = FALSE, seed = seed )
formula |
A formula as in lm() |
data |
A long data frame containing the locations, dates, covariates and the response variable. It has to have the locID and date. No missing values are allowed in the covariates. The order in this data.fame MUST be: spatial locations (1 to S) at time t=1, then locations (1 to S) at t=2 and so on. |
path |
Path with the name of the SpatialStreamNetwork object |
time_method |
A list specifying the temporal structure (ar = Autorregressive; var = Vector autorregression) and coumn in the data with the time variable. |
space_method |
A list defining if use or not of an SSN object and the spatial correlation structure. The second element is the spatial covariance structure. A 3rd element is a list with the lon and lat for Euclidean distance models. |
iter |
Number of iterations |
warmup |
Warm up samples |
chains |
Number of chains |
refresh |
Sampler refreshing rate |
net |
The network id (optional). Used when the SSN object contains multiple networks. |
addfunccol |
Variable to compute the additive function. Used to compute the spatial weights. |
loglik |
Logic parameter denoting if the loglik will be computed by the model. |
ppd |
Produce the posterior predictive distribution |
seed |
(optional) A seed for reproducibility |
Missing values are not allowed in the covariates and they must be imputed before using ssnbayes(). Many options can be found in https://cran.r-project.org/web/views/MissingData.html The pid in the data has to be consecutive from 1 to the number of observations. Users can use the SpatialStreamNetwork created with the SSN package. This will provide the spatial stream information used to compute covariance matrices. If that is the case, the data has to have point ids (pid) matching the ones in SSN distance matrices, so that a mapping can occur.
A list with the model fit
It returns a ssnbayes object (similar to stan returns). It includes the formula used to fit the model. The output can be transformed into the stanfit class using class(fits) <- c("stanfit").
Edgar Santos-Fernandez
## Not run: #options(mc.cores = parallel::detectCores()) # Import SpatialStreamNetwork object #path <- system.file("extdata/clearwater.ssn", package = "SSNbayes") #n <- SSN2::ssn_import(path, predpts = "preds", overwrite = TRUE) ## Imports a data.frame containing observations and covariates #clear <- readRDS(system.file("extdata/clear_obs.RDS", package = "SSNbayes")) #fit_ar <- ssnbayes(formula = y ~ SLOPE + elev + h2o_area + air_temp + sin + cos, # data = clear, # path = path, # time_method = list("ar", "date"), # space_method = list('use_ssn', c("Exponential.taildown")), # iter = 2000, # warmup = 1000, # chains = 3, # net = 2, # second network on the ssn object # addfunccol='afvArea') #space_method options examples #use list('no_ssn', 'Exponential.Euclid', c('lon', 'lat')) if no ssn object is available ## End(Not run)## Not run: #options(mc.cores = parallel::detectCores()) # Import SpatialStreamNetwork object #path <- system.file("extdata/clearwater.ssn", package = "SSNbayes") #n <- SSN2::ssn_import(path, predpts = "preds", overwrite = TRUE) ## Imports a data.frame containing observations and covariates #clear <- readRDS(system.file("extdata/clear_obs.RDS", package = "SSNbayes")) #fit_ar <- ssnbayes(formula = y ~ SLOPE + elev + h2o_area + air_temp + sin + cos, # data = clear, # path = path, # time_method = list("ar", "date"), # space_method = list('use_ssn', c("Exponential.taildown")), # iter = 2000, # warmup = 1000, # chains = 3, # net = 2, # second network on the ssn object # addfunccol='afvArea') #space_method options examples #use list('no_ssn', 'Exponential.Euclid', c('lon', 'lat')) if no ssn object is available ## End(Not run)
It requires the same number of observation/locations per day. It requires location id (locID) and points id (pid). The locID are unique for each site. The pid is unique for each observation. Missing values are allowed in the response and in the covariates. Missing values are imputed using the 'mtsdi' pacakge, with a 'glm' family.
ssnbayes2( formula = formula, family = family, data = data, osm_ssn = osm_ssn, path = NA, time_method = time_method, space_method = space_method, iter = 3000, warmup = 1500, chains = 3, cores = 3, refresh = max(iter/100, 1), net = NA, addfunccol = NA, loglik = FALSE, seed = seed )ssnbayes2( formula = formula, family = family, data = data, osm_ssn = osm_ssn, path = NA, time_method = time_method, space_method = space_method, iter = 3000, warmup = 1500, chains = 3, cores = 3, refresh = max(iter/100, 1), net = NA, addfunccol = NA, loglik = FALSE, seed = seed )
formula |
A formula as in lm() |
family |
A description of the response distribution and link function to be used in the model. Must be one of: 'gaussian', 'binomial', 'poisson". |
data |
A long data frame containing the locations, dates, covariates and the response variable. It has to have the locID and date. No missing values are allowed in the covariates. The order in this data.fame MUST be: spatial locations (1 to S) at time t=1, then locations (1 to S) at t=2 and so on. |
osm_ssn |
An 'osm_ssn' generated using 'generate_osm_ssn()'. |
path |
File path with the name of the SpatialStreamNetwork object |
time_method |
A list specifying the temporal structure (ar = Autorregressive; var = Vector autorregression) and column in the data with the time variable. |
space_method |
A list, first element must be one of 'use_ssn', 'use_osm_ssn', 'no_ssn. Whether to use SSN or osm ssn, and the spatial correlation structure. The second element is the spatial covariance structure. A 3rd element is a list with the lon and lat for Euclidean distance models. |
iter |
Number of iterations |
warmup |
Warm up samples |
chains |
Number of chains |
cores |
Number of cores |
refresh |
Sampler refreshing rate |
net |
The network id (optional). Used when the SSN object cotains multiple networks. |
addfunccol |
Variable to compute the additive function. Used to compute the spatial weights. |
loglik |
Logic parameter denoting if the loglik will be computed by the model. |
seed |
(optional) A seed for reproducibility |
Missing values on the covariates and response can be imputed by default using mtsdi::mnimput(). We strongly recommend the documentation for this function be read before use. The pid in the data has to be consecutive from 1 to the number of observations. Users can use the SpatialStreamNetwork created with the SSN package. This will provide the spatial stream information used to compute covariance matrices. If that is the case, the data has to have point ids (pid) matching the ones in SSN distance matrices, so that a mapping can occur.
A list with the model fit
It returns a ssnbayes object (similar to stan returns). It includes the formula used to fit the model. The output can be transformed into the stanfit class using class(fits) <- c("stanfit").
Edgar Santos-Fernandez
## Not run: #options(mc.cores = parallel::detectCores()) #require(SSNdata) #formula = temp ~ SLOPE + elev + h2o_area + air_temp + sin + cos ## Imports a data.frame containing observations and covariates #clear <- readRDS(system.file("extdata/clear_obs.RDS", package = "SSNdata")) #family <- "gaussian" # # If using osm_ssn: # # Generate osm_ssn # clear_osm_ssn <- generate_osm_ssn(clear, "long", "lat", root_loc = 12, plot_network = TRUE) # fit_ar <- ssnbayes2(formula, # data = clear, # osm_ssn = clear_osm_ssn, # family = family, # time_method = list("ar", "date"), # space_method = list('use_osm_ssn', c("Exponential.taildown")), # iter = 2000, # warmup = 1000, # chains = 3, # cores = 3) # # # # If not using osm_ssn, and instead using SSN: # Import SpatialStreamNetwork object # path <- system.file("extdata/clearwater.ssn", package = "SSNbayes") # fit_ar <- ssnbayes2(formula = formula, # data = clear, # path = path, # family = family, # time_method = list("ar", "date"), # space_method = list('use_ssn', c("Exponential.taildown")), # iter = 2000, # warmup = 1000, # chains = 3, # cores = 3, # net = 2, # second network on the ssn object # addfunccol='afvArea') # # # # #space_method options examples # #use list('no_ssn', 'Exponential.Euclid', c('lon', 'lat')) if no ssn object is available ## End(Not run)## Not run: #options(mc.cores = parallel::detectCores()) #require(SSNdata) #formula = temp ~ SLOPE + elev + h2o_area + air_temp + sin + cos ## Imports a data.frame containing observations and covariates #clear <- readRDS(system.file("extdata/clear_obs.RDS", package = "SSNdata")) #family <- "gaussian" # # If using osm_ssn: # # Generate osm_ssn # clear_osm_ssn <- generate_osm_ssn(clear, "long", "lat", root_loc = 12, plot_network = TRUE) # fit_ar <- ssnbayes2(formula, # data = clear, # osm_ssn = clear_osm_ssn, # family = family, # time_method = list("ar", "date"), # space_method = list('use_osm_ssn', c("Exponential.taildown")), # iter = 2000, # warmup = 1000, # chains = 3, # cores = 3) # # # # If not using osm_ssn, and instead using SSN: # Import SpatialStreamNetwork object # path <- system.file("extdata/clearwater.ssn", package = "SSNbayes") # fit_ar <- ssnbayes2(formula = formula, # data = clear, # path = path, # family = family, # time_method = list("ar", "date"), # space_method = list('use_ssn', c("Exponential.taildown")), # iter = 2000, # warmup = 1000, # chains = 3, # cores = 3, # net = 2, # second network on the ssn object # addfunccol='afvArea') # # # # #space_method options examples # #use list('no_ssn', 'Exponential.Euclid', c('lon', 'lat')) if no ssn object is available ## End(Not run)
Fits a deep learning spatio-temporal model using LSTM layers with attention mechanisms and spatial embeddings. This model is designed for detecting anomalies in high-frequency sensor data by capturing both spatial dependencies (via distance-based embeddings) and temporal dynamics (via recurrent layers with attention). Uncertainty is quantified via Monte Carlo dropout during inference.
starn( formula, data, spat_matrix, time_col, loc_col, timesteps = 11, spat_embedding_dim = 10, time_method = "lstm", space_method = "attention", batch_size = 64, epochs = 100, validation_split = 0.2, n_iter = 100 )starn( formula, data, spat_matrix, time_col, loc_col, timesteps = 11, spat_embedding_dim = 10, time_method = "lstm", space_method = "attention", batch_size = 64, epochs = 100, validation_split = 0.2, n_iter = 100 )
formula |
A formula specifying the response and predictors, e.g., |
data |
A data frame containing the spatio-temporal sensor data. |
spat_matrix |
A spatial distance matrix between sensor locations (square and symmetric). |
time_col |
Name of the column indicating temporal ordering (unquoted or quoted). |
loc_col |
Name of the column indicating spatial location IDs (unquoted or quoted). |
timesteps |
Length of the temporal window used to construct sequences for the LSTM model (default: 11). |
spat_embedding_dim |
Dimension of the spatial embedding extracted via classical MDS (default: 10). |
time_method |
Method used for temporal modeling (currently supports "lstm"). |
space_method |
Method used for spatial integration (currently supports "attention"). |
batch_size |
Training batch size (default: 64). |
epochs |
Number of training epochs (default: 100). |
validation_split |
Proportion of training data used for validation (default: 0.2). |
n_iter |
Number of forward passes with dropout for uncertainty estimation (default: 100). |
The function uses bidirectional LSTMs followed by an attention mechanism to extract temporal features, and incorporates spatial information via embeddings obtained through classical multidimensional scaling (MDS). Uncertainty is quantified using Monte Carlo dropout, and anomalies are flagged when observed values fall outside the estimated 90
A list with components:
modelThe trained Keras model object.
historyTraining history returned by Keras.
dataThe input data frame augmented with predictions, prediction intervals, and anomaly indicators.
Santos-Fernandez, E., et al. (2025). *New Bayesian and deep learning spatio-temporal models can reveal anomalies in sensor data more effectively*. (in review).
## Not run: # Load example data from package obs_path <- system.file("extdata", "obs_data.RDS", package = "SSNbayes") dm_path <- system.file("extdata", "Cov_mat.rds", package = "SSNbayes") obs_data <- readRDS(obs_path) dm <- readRDS(dm_path) spat_matrix <- dm$H # Fit the STARN model fit_starn <- starn( formula = y_impute ~ X1 + X2 + X3, data = obs_data, time_col = 'date', loc_col = 'locID', spat_matrix = spat_matrix, timesteps = 11, time_method = "lstm", space_method = "attention", batch_size = 64, epochs = 5, validation_split = 0.2 ) ## End(Not run)## Not run: # Load example data from package obs_path <- system.file("extdata", "obs_data.RDS", package = "SSNbayes") dm_path <- system.file("extdata", "Cov_mat.rds", package = "SSNbayes") obs_data <- readRDS(obs_path) dm <- readRDS(dm_path) spat_matrix <- dm$H # Fit the STARN model fit_starn <- starn( formula = y_impute ~ X1 + X2 + X3, data = obs_data, time_col = 'date', loc_col = 'locID', spat_matrix = spat_matrix, timesteps = 11, time_method = "lstm", space_method = "attention", batch_size = 64, epochs = 5, validation_split = 0.2 ) ## End(Not run)