| Title: | Hidden Markov Modelling of Hydrological State Change |
|---|---|
| Description: | Identifies regime changes in streamflow runoff not explained by variations in precipitation. The package builds a flexible set of Hidden Markov Models of annual, seasonal or monthly streamflow runoff with precipitation as a predictor. Suites of models can be built for a single site, ranging from one to three states and each with differing combinations of error models and auto-correlation terms. The most parsimonious model is easily identified by AIC, and useful for understanding catchment drought non-recovery: Peterson TJ, Saft M, Peel MC & John A (2021) <doi:10.1126/science.abd5085>. |
| Authors: | Tim Peterson [aut, cre, cph] (ORCID: <https://orcid.org/0000-0002-1885-0826>), Thomas Westfall [aut] (ORCID: <https://orcid.org/0009-0000-0529-881X>) |
| Maintainer: | Tim Peterson <[email protected]> |
| License: | GPL-3 |
| Version: | 0.2.0.0 |
| Built: | 2026-05-11 08:29:29 UTC |
| Source: | https://github.com/peterson-tim-j/hydrostate |
hydroState is an R-package that identifies regime changes in streamflow time-series that is not explained by variations in precipitation.
For details of the mathematics, and its use in understanding catchment drought non-recovery, see:
Peterson TJ, Saft M, Peel MC & John A (2021), Watersheds may not recover from drought, Science, DOI: 10.1126/science.abd5085
The package allows a flexible set of Hidden Markov Models (HMMs) of annual, seasonal or monthly time-step to be built and which includes precipitation as a predictor of streamflow. Suites of models can be build for a single catchment, ranging from from one to three states and each with differing combinations of error models and auto-correlation terms, allowing the most parsimonious model to easily be identified (by AIC). The entire package is written in R S4 object oriented code with user facing functions and vignettes.
build — Builds hydroState model
build.all — Builds all hydroState models
fit.hydroState — Fit hydroState model(s)
setInitialYear — Sets state names given initial year
plot.hydroState — Plot states or pseudo residuals over time
get.residuals — Get pseudo residuals
get.states — Get states
check — Check reliability of state predictions
get.AIC — Get AIC
get.seasons — Get data as seasons
The package contains a default hydroState model object that explains streamflow as a function of precipitation using a linear model. Once the model object is built (build), the model is fitted (fit.hydroState) to determine the most likely rainfall-runoff state at each time-step. To assess the adequacy of the fit, the residuals are plotted (plot, and an adequate fit requires the residuals to be normally distributed, uniform, with minimal correlation and minimal trends. The resulting runoff states from the fitted model can then be evaluated over time (plot) and even exported (get.states) with the state values, confidence intervals, and conditional probabilities at each time-step. Input data requires a dataframe with catchment average runoff and precipitation at annual, seasonal, or monthly timesteps, and gaps with missing data are permitted. An example of this workflow with the default model is demonstrated within the Getting Started article.
To better explain the rainfall-runoff relationship, the default model can be adjusted by selecting various items within the build function. These include:
data.transform: transform streamflow observations in order to reduce skew: 'boxcox', 'log', 'burbidge', or 'none'
parameters: account for auto-correlation through including the degree of auto-correlation: 'AR1', 'AR2', or 'AR3'
state.shift.parameters: assign either the intercept, slope, or auto-correlation parameter as state dependent parameter
error.distribution: adjust the error distribution with 'normal', 'gamma', or 'truc.normal'
seasonal.parameters: account for intra-annual varation within the rainfall-runoff relationship
transition.graph: set the number of possible states in the model (1, 2, or 3).
An example of how to adjust the default model is demonstrated within the Adjust the default state model article and Seasonal and monthly models article.
There is an additional option to construct all possible types of models using the build.all, and compare them using the same fit.hydroState function. The most likely model can be selected based on the AIC where the best model will have the lowest AIC. An example of this is demonstrated at the end of the Adjust the default state model article
and Seasonal and monthly models article. To get stated, it is recommended to evaluate the default model at first with one state and again with two states.
The package development was funded by the Victorian Government The Department of Energy, Environment, and Climate Action (https://www.water.vic.gov.au/).
Maintainer: Tim Peterson [email protected] (ORCID) [copyright holder]
Authors:
Thomas Westfall [email protected] (ORCID)
Useful links:
Report bugs at https://github.com/peterson-tim-j/HydroState/issues
build builds a hydrostate model object with either a default model or the model can be specified with options from below. Every model depends on a linear base model where streamflow, , is a function of precipitation, : . The default model is a variant of this base linear model with state shifts expected in the intercept, , and standard deviation, , of the rainfall-runoff relationship. There are additional options to adjust this model with auto-correlation terms and seasonal parameters that can be both independent of state or change with state. The number of states and assumed error distribution can also be selected. After the model is built, the hydroState model is ready to be fitted with fit.hydroState
build( input.data = data.frame(year = c(), flow = c(), precip = c()), data.transform = "boxcox", parameters = c("a0", "a1", "std"), seasonal.parameters = NULL, state.shift.parameters = c("a0", "std"), error.distribution = "truc.normal", flickering = FALSE, transition.graph = matrix(TRUE, 2, 2) )build( input.data = data.frame(year = c(), flow = c(), precip = c()), data.transform = "boxcox", parameters = c("a0", "a1", "std"), seasonal.parameters = NULL, state.shift.parameters = c("a0", "std"), error.distribution = "truc.normal", flickering = FALSE, transition.graph = matrix(TRUE, 2, 2) )
input.data |
dataframe of annual, seasonal, or monthly runoff and precipitation observations. Gaps with missing data in either streamflow or precipitation are permitted. Monthly data is required when using |
data.transform |
character sting with the method of transformation. The default is 'boxcox'. Other options: 'log', 'burbidge', 'none' |
parameters |
character vector of parameters to construct model. Required and default: |
seasonal.parameters |
character vector of one or all parameters ( |
state.shift.parameters |
character vector of one or all parameters ( |
error.distribution |
character string of the distribution in the HMM error. Default is 'truc.normal'. Others include: 'normal' or 'gamma' |
flickering |
logical |
transition.graph |
matrix given the number of states. Default is a 2-state matrix (2 by 2): |
build
There are a selection of items to consider when defining the rainfall-runoff relationship and investigating state shifts in this relationship. hydroState provides various options for modelling the rainfall-runoff relationship.
Data gaps with input.data: When there is missing input.data in either the dependent variable, streamflow, or independent variable, precipitation, the emissions probability of the missing time-step is set equal to one. This essentially ignores the missing periods. The time step after the missing period has a state probability dependent on the length of the gap. The larger the gap, the closer the state probability gets to approaching a finite probability near zero (as the transition probabilities are recursively multiplied). When the model has auto-correlation terms and there are gaps in the dependent variable, the auto-correlation function restarts at the beginning of each continuous period after the gap. This ignores auto-correlation at the first time steps after the gap. For instance, an 'AR1' model would ignore the contribution of the prior time step for the first (1) observation after the gap.
Transform Observations with data.transform: Transforms streamflow observations to remove heteroscedasticity. Often there is skew within hydrologic data. When defining relationships between rainfall-runoff, this skew results in an unequal variance in the residuals, heteroscedasticity. Transforming streamflow observations is often required. There are several options to transform observations. Since the degree of transformation is not typically known, boxcox is the default. Other options include: log, burbidge, and of course, none when no transformation is performed.
Model Structure with parameters and seasonal.parameters: The structure of the model depends on the parameters. hydroState simulates runoff, , as being in one of a finite states, , at every time-step, , depending on the distribution of states at prior time steps. This results in a runoff distribution for each state that can vary over time (). The model defines the relationship that is susceptible to state shifts with precipitation, , as a predictor. This takes the form as a simple linear model :
where and are constant parameters. These parameters and the model error, , establish the rainfall-runoff relationship and are required parameters for every built model object. These are the default parameters: c('a0', 'a1', 'std').
Auto-correlation parameters: The relationship may contain serial correlation and would be better defined with an auto-regressive term:
where AR1 is the lag-1 auto-correlation term. Either, lag-1: AR1, lag-2: AR2, and lag-3: AR3 auto-correlation coefficients are an option as additional parameters to better define the rainfall-runoff relationship.
Sub-annual analysis with seasonal.parameters: Additional options include explaining the seasonal rainfall-runoff relationship with a sinusoidal function that better defines either of the constant parameters or error () throughout the year, i.e:
where is an integer month at . Monthly streamflow and precipitation are required as input.data for the sub-annual analysis.
State Dependent Parameters with state.shift.parameters: These are state dependent parameters where they are subject to shift in order to better explain the state of streamflow over time. Any or all of the previously chosen parameters can be selected (a_0, a_1, std, AR1, AR2, AR3). The default model evaluates shifts in the rainfall-runoff relationship with a_0 std as state dependent parameters.
Distribution of the Residuals with error.distribution: The distribution of the residuals (error) within a state of the model can be chosen to reduce skew and assist with making models statistically adequate (see plot(pse.residuals = TRUE)). Either normal: normal, truncated normal: truc.normal, or gamma: gamma distributions are acceptable. These error distribution ensures streamflow is greater than zero Q > 0, and specifically for truc.normal greater than or equal to zero Q >= 0. The default is truc.normal. Sub-annual models are restricted to only a gamma distribution.
Markov flickering with flickering: When flickering is FALSE, the markov avoids state shifts for very short duration, and hence for a state shift to occur it should last for an extended period. The default is FALSE. If TRUE, flickering between more states is more sensitive. For further explanation on this method, see: Lambert et al., 2003. The current form of the Markov model is homogeneous where the transition probabilities are time-invariant.
Number of States with transition.graph: The number of possible states in the rainfall-runoff relationship and transition between the states is selected with the transition.graph. The default is a 2-state model in a 2 by 2 unstructured matrix with a TRUE transition to and from each state (i.e. matrix(TRUE,2,2)). hydroState accepts 1-state up to 3-states (i.e. for 3-state unstructured transition graph: matrix(TRUE,3,3)). The unstructured transition graph allows either state to remain in the current state or transition between any state. For the 3-state transition graph, one may want to assume the transitions can only occur in a particular order, as in (Very low -> Low -> Normal->) rather than (Very low <-> Low <-> Normal <-> Very low). Thus, a structured graph is also acceptable (3-state structured: matrix(c(TRUE,TRUE,FALSE,FALSE,TRUE,TRUE,TRUE,FALSE,TRUE),3,3)). More details in Supplementary Materials of (Peterson TJ, Saft M, Peel MC & John A (2021), Watersheds may not recover from drought, Science, DOI: doi:10.1126/science.abd5085).
A built hydroState model object ready to be fitted with fit.hydroState()
# Load data data(streamflow_annual_221201) ## Build default annual hydroState model model = build(input.data = streamflow_annual_221201) # OR ## Build annual hydroState model with specified objects # Build hydroState model with: 2-state, normal error distribution, # 1-lag of auto-correlation, and state dependent parameters ('a1', 'std') model = build(input.data = streamflow_annual_221201, data.transform = 'boxcox', parameters = c('a0','a1','std','AR1'), state.shift.parameters = c('a1','std'), error.distribution = 'normal', flickering = FALSE, transition.graph = matrix(TRUE,2,2))# Load data data(streamflow_annual_221201) ## Build default annual hydroState model model = build(input.data = streamflow_annual_221201) # OR ## Build annual hydroState model with specified objects # Build hydroState model with: 2-state, normal error distribution, # 1-lag of auto-correlation, and state dependent parameters ('a1', 'std') model = build(input.data = streamflow_annual_221201, data.transform = 'boxcox', parameters = c('a0','a1','std','AR1'), state.shift.parameters = c('a1','std'), error.distribution = 'normal', flickering = FALSE, transition.graph = matrix(TRUE,2,2))
build.all builds all possible combinations of hydroState models. The same fields are available as in build in order to specify the type of models to be built. After all models are built, they are fitted using the same fit.hydroState() function.
build.all( input.data = data.frame(year = c(), flow = c(), precip = c()), data.transform = NULL, parameters = NULL, seasonal.parameters = NULL, state.shift.parameters = NULL, error.distribution = NULL, flickering = FALSE, transition.graph = NULL, summary.table = NULL, siteID = NULL )build.all( input.data = data.frame(year = c(), flow = c(), precip = c()), data.transform = NULL, parameters = NULL, seasonal.parameters = NULL, state.shift.parameters = NULL, error.distribution = NULL, flickering = FALSE, transition.graph = NULL, summary.table = NULL, siteID = NULL )
input.data |
dataframe of annual, seasonal, or monthly runoff and precipitation observations. Gaps with missing data in either streamflow or precipitation are permitted, and the handling of them is further discussed in |
data.transform |
character string of method of transformation. If empty, the default builds all possible combinations of models with |
parameters |
character vector of parameters to determine model form. If empty, the default builds all possible combinations of model forms. |
seasonal.parameters |
character vector of parameters with sinusoidal function to represent seasonal variation. Requires monthly or seasonal data. If empty and monthly or seasonal data is given, the default builds all possible combinations of models with a seasonal parameter for each and all parameters. |
state.shift.parameters |
character vector of one or all parameters to identify state dependent parameters. Only one set of parameters permitted. If empty, the default builds all possible model combinations with |
error.distribution |
character string of the distribution in the HMM error. If empty, the default builds models with all possible combinations of error distribution: |
flickering |
logical |
transition.graph |
matrix given the number of states. If empty, the default builds models with all possible combinations of states:
1-state matrix (1 by 1): |
summary.table |
data frame with a table summarizing all built models and corresponding reference model. From function |
siteID |
character string of site identifier. |
build.all
All possible combinations of hydroState models are built for each auto-correlation lag and residual distribution from 1 to 3 states for a specified data transformation. This allows for investigation of state changes in the state.shift.parameters: the intercept c('a0', 'std') or slope c('a1', 'std'). To reduce the number of models in the search, specify which field(s) to remain constant. For example, to investigate the best model with the number of auto-correlation terms and number of states with a boxcox data transform and gamma distribution of the residuals, set data.transform to boxcox and error.distribution to gamma. If no fields are specified, all possible model combinations are built. If investigating state shifts in the intercept a0 and slope a1, it is recommended to build and fit the model combinations separately.
A list of built hydroState models with every combination of objects ready to be fitted
# Load data data(streamflow_annual_221201) # Build all annual models with state shift in intercept 'a0' all.annual.models = build.all(input.data = streamflow_annual_221201, state.shift.parameters = c('a0','std'), siteID = '221201') # OR # Build all annual models with state shift in slope 'a1' all.annual.models = build.all(input.data = streamflow_annual_221201, state.shift.parameters = c('a1','std'), siteID = '221201')# Load data data(streamflow_annual_221201) # Build all annual models with state shift in intercept 'a0' all.annual.models = build.all(input.data = streamflow_annual_221201, state.shift.parameters = c('a0','std'), siteID = '221201') # OR # Build all annual models with state shift in slope 'a1' all.annual.models = build.all(input.data = streamflow_annual_221201, state.shift.parameters = c('a1','std'), siteID = '221201')
check After fitting the model, the reliability of the estimated states can be assessed by generating synthetic state sequences and then assessing how well the model identified them.
check(model, n.samples = 1e+05)check(model, n.samples = 1e+05)
model |
fitted |
n.samples |
integer of samples to re-sample. Default is 100000. |
check
This validates the model's states at each time-step through re-sampling the input data and re-running the Viterbi algorithm. The input data is duplicated 100 times, and a synthetic series is generated from the model with sample states. This provides a time series of the transformed streamflow observations that can be compared with observations of the ‘known’ state. The Viterbi states of the re-sampled transformed observations are inferred, and the probability of the inferred state equaling the 'known' state is calculated.
A data frame is returned with a matrix depending on the number of states. For a 2 state model, a 2x2 matrix is returned. The diagonal cell estimates the probability of correctly identifying that state. The off diagonals estimate the probability of incorrectly identifying a state that in state 2.
## Check reliability of state predictions (>5s to run) check(model = model.annual.fitted.221201)## Check reliability of state predictions (>5s to run) check(model = model.annual.fitted.221201)
fit.hydroState fits a single hydroState model (build) or multiple models (build.all) using global optimization by differential evolution DEoptim library. If fitting all models be sure to install and load the parallelly library. The fitting of all models may take hours or days, but the calibration can occur in parallel if the parallelly library is installed and loaded.
fit.hydroState( model, pop.size.perParameter = 10, max.generations = 500, doParallel = FALSE, ... )fit.hydroState( model, pop.size.perParameter = 10, max.generations = 500, doParallel = FALSE, ... )
model |
built hydroState model object, hydroState.allModels object, or hydroState.subAnnual.allModels object |
pop.size.perParameter |
integer that should be greater than or equal to the number of parameters in the model. The default is '10' and is sufficient for all models. |
max.generations |
integer that will stop the optimizer when set number of generations are reached. The default is '500'. |
doParallel |
TRUE/FALSE to perform fitting in parallel on all computer cores. Default is FALSE |
... |
additional options to change the optimization settings: reltol, print.iterations, etc. from the DEoptim library |
fit.hydroState
After a hydroState model object is built, the model is ready to fit to the observed streamflow through minimizing the negative log-likelihood function too calibrate model parameters. The only required input is the given built hydroState model object or hydroState.allModels object (all models from build.all = TRUE). When fitting all models, the models are fitted from least to most complex (least to max amount of parameters). Each model has a minimum of 5 and maximum of 20 calibration attempts to outperform the prior reference model else the model is rejected. For instance 'model.1State.normal.log.AR0' contains 3-parameters and is the reference model for 'model.1State.normal.log.AR1' which contains 4-parameters. The objective function of 'model.1State.normal.log.AR1' must calibrate the model with a lower negative log-likelihood than 'model.1State.normal.log.AR0'. These reference models are pre-defined, but this function allows the user to edit the reference models in the data.frame if needed using summary. Details on the likelihood function is as follows:
The likelihood function is estimated as:
where:
is the initial state distribution, the initial probability of being in each state:
is the x diagonal emissions matrix of the probability density for each state using a lower tail truncated Gaussian distribution or a two-parameter Gamma distribution:
where is the probability density function for the standard normal distribution, is the cumulative distribution function for the standard normal distribution, is the shape parameter, is the scale parameter, and is the gamma function. For more details, refer to pg. 8–17 in Supplementary Materials of (Peterson TJ, Saft M, Peel MC & John A (2021), Watersheds may not recover from drought, Science, DOI: doi:10.1126/science.abd5085).
is the transition matrix
is the number of time-steps.
A fitted hydroState model
# Load data data(streamflow_annual_221201) ## Build default annual hydroState model model = build(input.data = streamflow_annual_221201) ## Fit built model (runtime ~ 14 sec) model = fit.hydroState(model) ## Fit all built models (runtime > several hours) # Load data data(streamflow_annual_221201) ## Build all annual models all.annual.models = build.all(input.data = streamflow_annual_221201, siteID = '221201') # Fit all (runtime > several hours) all.annual.models = fit.hydroState(all.annual.models)# Load data data(streamflow_annual_221201) ## Build default annual hydroState model model = build(input.data = streamflow_annual_221201) ## Fit built model (runtime ~ 14 sec) model = fit.hydroState(model) ## Fit all built models (runtime > several hours) # Load data data(streamflow_annual_221201) ## Build all annual models all.annual.models = build.all(input.data = streamflow_annual_221201, siteID = '221201') # Fit all (runtime > several hours) all.annual.models = fit.hydroState(all.annual.models)
get.AIC retrieves Akaike information criteria from a fitted hydroState model object or all models.
get.AIC(model)get.AIC(model)
model |
fitted |
get.AIC
The AIC is the negative log-likelihood of the model plus a penalty for model parameters. This function can be performed on a single model or a selection of models to find the lowest AIC of the set.
AIC value of a single model or a list variable of AIC values for al models
# Load fitted model data(model.annual.fitted.221201) ## AIC of a single model get.AIC(model.annual.fitted.221201) ## Lowest AIC of a model set get.AIC(all.models.annual.fitted.407211)# Load fitted model data(model.annual.fitted.221201) ## AIC of a single model get.AIC(model.annual.fitted.221201) ## Lowest AIC of a model set get.AIC(all.models.annual.fitted.407211)
The pseudo residuals were derived from the conditional probabilities of the observations. At each time-step, the pseudo residual is the probability of an observation occurring given the prior observations and latter observations.
get.residuals(model)get.residuals(model)
model |
fitted hydroState model object. |
get.residuals
get.residuals retrieves residuals from the fitted model and exports them as a data frame.
Data frame of residuals for each time-step
# Load fitted model data(model.annual.fitted.221201) ## Get residuals in a dataframe get.residuals(model = model.annual.fitted.221201)# Load fitted model data(model.annual.fitted.221201) ## Get residuals in a dataframe get.residuals(model = model.annual.fitted.221201)
Aggregates monthly data to 4 seasons in a year.
get.seasons( input.data = data.frame(year = c(), month = c(), flow = c(), precip = c()) )get.seasons( input.data = data.frame(year = c(), month = c(), flow = c(), precip = c()) )
input.data |
dataframe of monthly runoff and precipitation observations. Gaps with missing data in either streamflow or precipitation are permitted, and the handling of them is further discussed in |
get.seasons
This function takes sums monthly runoff and precipitation observations into 4 seasons of a year.
A dataframe of seasonal observations with an additional column counting the number of months in each season.
# Load data data(streamflow_monthly_221201) # aggregate monthly data to seasonal streamflow_seasonal_221201 = get.seasons(streamflow_monthly_221201)# Load data data(streamflow_monthly_221201) # aggregate monthly data to seasonal streamflow_seasonal_221201 = get.seasons(streamflow_monthly_221201)
get.states uses the Viterbi algorithm to globally decode the model and estimate the most probable sequence of states.
get.states(model)get.states(model)
model |
fitted |
get.states
These dataframe of results include:
time-step: year and possibly either season or month for subannual analysis
Viterbi State Number: state number (i.e. 1, 2, or 3) to differentiate states
Obs. flow: streamflow observations
Viterbi Flow: flow values of the Viterbi state including the 5\
Normal State Flow: flow values of the normal state including the 5\
Conditional Prob: conditional probabilities for each state show the probability of remaining in the given state. When the conditional probability is closer to 1, there is a higher probability that hydroState model remains in that state for the next time-step.
Emission Density: emission density for each state is the result of multiplying the conditional probabilities by the transition probabilities at each timestep.
data frame of results to evaluate the rainfall-runoff states over time
# Load fitted model data(model.annual.fitted.221201) ## Set initial year to set state names model.annual.fitted.221201 = setInitialYear(model = model.annual.fitted.221201, initial.year = 1990) ## Get states model.annual.fitted.221201.states = get.states(model = model.annual.fitted.221201)# Load fitted model data(model.annual.fitted.221201) ## Set initial year to set state names model.annual.fitted.221201 = setInitialYear(model = model.annual.fitted.221201, initial.year = 1990) ## Get states model.annual.fitted.221201.states = get.states(model = model.annual.fitted.221201)
plot produces several figures to visualize pseudo residuals or results of the markov states over time. setInitialYear is required before plot. It is recommend to evaluate the pseudo residuals before the markov states. The pseudo residuals are the probability of an observation occurring at each time-step given the prior observations and latter observations, and these are derived from the conditional probabilities of the observations. The markov states are from the Viterbi algorithm globally decoding the model to estimate the most probable sequence of states.
## S3 method for class 'hydroState' plot( x, ..., pse.residuals = FALSE, ind.variable = FALSE, dep.variable = FALSE, dep.variable.transformed = FALSE, cond.state.prob = FALSE, siteID = NULL, file = NULL )## S3 method for class 'hydroState' plot( x, ..., pse.residuals = FALSE, ind.variable = FALSE, dep.variable = FALSE, dep.variable.transformed = FALSE, cond.state.prob = FALSE, siteID = NULL, file = NULL )
x |
is the fitted hydroState model object. |
... |
additional arguments passed for plotting, none available at this time. |
pse.residuals |
option to plot pseudo residuals. Default is FALSE. |
ind.variable |
option to plot independent variable over time. |
dep.variable |
option to plot dependent variable and states over time. |
dep.variable.transformed |
option to plot transformed dependent variable and states over time. |
cond.state.prob |
option to plot the conditional state probabilities over time for each state. |
siteID |
character string of catchment identifier (i.e. gauge ID). Default is NULL. Only recommended when exporting (i.e, file = "filename.pdf"). |
file |
character string of file directory/name to export plots as a pdf: "flow.state.plots.407211.pdf". Default is NULL, no pdf file is exported. |
plot
plot produces five figures of psuedo residuals OR up to four figures of the results from the fitted hydroState model. When the pse.residuals is FALSE, the default plot produces all four markov state plots (ind.variable, dep.variable, dep.variable.transformed, cond.state.prob). Figures are more easily viewed as an exported pdf when directory/file name is given (i.e. file = "flow.state.plots.siteID.pdf").
psuedo residual figures
A) Time-series of normal-pseudo residuals to ensure the residuals each year are within the confidence intervals.
B) Auto-correlation function (ACF) of normal-pseudo residuals to ensure there is minimal serial correlation in residuals. Lag spikes should be below confidence interval at each lag (except 0).
C) Histogram of uniform-pseudo residuals should show uniform distribution (equal frequency for each residual value)
D) Histogram of normal-pseudo residuals should show normal distribution centered on zero and with no skew
E) Quantile-Quantile (Q-Q) plot where normal-pseudo residuals vs. theoretical quantities should align on the diagonal line. The last plot contains the Akaike information criterion (AIC) and Shapiro-Wilk p-value. The AIC is an estimator to determine the most parsimonious, best performing model given the number of parameters. When comparing models, the lowest AIC is the best performing model. Shapiro-Wilks test for normality in the residuals and a p-value greater than 0.05 (chosen alpha level) indicates the residuals are normally distributed; the null hypothesis that the residuals are normally distributed is not rejected.
markov state figures
A) independent variable: precipitation
B) dependent variable and states: streamflow observations, most likely state, and relative normal state estimate
C) transformed dependent variable and states: transformed streamflow observations and most likely state
D) conditional state probabilities for each state: probability of hydroState model remaining in given state
These figures are often large, and below are a few common errors when the plotting window is too small. Exporting the plots as a file is recommend for the pseudo residual figure (file = "file name").
"Error in plot.new() : figure margins too large": reset plot window with "dev.off()", enlarge plot area and re-run plot.residuals.
"Error in par(op) : invalid value specified for graphical parameter "pin" if the R plot window is not reset with "dev.off", an additional plot.residuals attempt will result in this error.
plots to evaluate rainfall-runoff states over time along with observations and the conditional probabilities of each state.
# Load fitted model data(model.annual.fitted.221201) ## Set initial year to set state names model.annual.fitted.221201 = setInitialYear(model = model.annual.fitted.221201, initial.year = 1990) ## Plot only residuals plot(model.annual.fitted.221201, pse.residuals = TRUE) ## Plot all markov state figures (default) plot(model.annual.fitted.221201) ## Plot only dependent variable transformed with markov states plot(model.annual.fitted.221201, dep.variable.transformed = TRUE)# Load fitted model data(model.annual.fitted.221201) ## Set initial year to set state names model.annual.fitted.221201 = setInitialYear(model = model.annual.fitted.221201, initial.year = 1990) ## Plot only residuals plot(model.annual.fitted.221201, pse.residuals = TRUE) ## Plot all markov state figures (default) plot(model.annual.fitted.221201) ## Plot only dependent variable transformed with markov states plot(model.annual.fitted.221201, dep.variable.transformed = TRUE)
sets the state names for each time-step relative to the initial year given
setInitialYear(model, initial.year)setInitialYear(model, initial.year)
model |
fitted hydroState model object. |
initial.year |
integer with year (YYYY). Default is first year in input.data. |
setInitialYear
hydroState assigns names to the computed states. This requires choosing an initial year where the state value from that year will be named 'Normal'. Other state values will be given names relative to the state value in the initial year. The choice of the initial year does not affect results. It is a means to more easily interpret the difference in state values relative to each other. It is best to choose a year based on the question being asked. For example, in testing the impact of drought, a year before the beginning of the drought, 1990, was selected as an initial year when conditions were considered 'Normal' (Peterson TJ, Saft M, Peel MC & John A (2021), Watersheds may not recover from drought, Science, DOI: doi:10.1126/science.abd5085)
A fitted hydroState model object with state names for each time-step ready for plot
# Load fitted model data(model.annual.fitted.221201) ## Set initial year to set state names model.annual.fitted.221201 = setInitialYear(model = model.annual.fitted.221201, initial.year = 1990)# Load fitted model data(model.annual.fitted.221201) ## Set initial year to set state names model.annual.fitted.221201 = setInitialYear(model = model.annual.fitted.221201, initial.year = 1990)
summary outputs a summary table of all built hydroState models and allows users to edit the reference models for calibration.
## S3 method for class 'hydroState.allModels' summary(object, ...)## S3 method for class 'hydroState.allModels' summary(object, ...)
object |
hydroState.allModels object with a list of models from |
... |
No additional input required |
summary
For every model object in build.all, there is a reference model for calibration. The reference model is a slightly simpler model with one less parameter. During calibration with fit, the model performance must exceed the the performance of the reference model else the model is rejected. This function is used to output the summary.table and adjust the reference models. Afterwards, all models can be re-build with including this summary.table in the build.all function.
A data.frame with a summary table of all models and reference models
# Show summary table of all fitted model details and reference models all.models.ref.table = summary(all.models.annual.fitted.407211)# Show summary table of all fitted model details and reference models all.models.ref.table = summary(all.models.annual.fitted.407211)