Title: Hidden Markov Modelling of Hydrological State Change
Version: 0.2.0.0
Maintainer: Tim Peterson <tim.peterson@monash.edu>
Depends: R (≥ 3.5.0)
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>.
Imports: methods, DEoptim, sn, truncnorm, diagram, padr, zoo, graphics, checkmate
BugReports: https://github.com/peterson-tim-j/HydroState/issues
URL: https://github.com/peterson-tim-j/HydroState, https://peterson-tim-j.github.io/HydroState/
License: GPL-3
Encoding: UTF-8
LazyData: true
RoxygenNote: 7.3.1
ByteCompile: true
Suggests: knitr, rmarkdown, testthat (≥ 3.0.0), V8
Config/testthat/edition: 3
Collate: 'abstracts.R' 'parameters.R' 'Qhat.boxcox.R' 'Qhat.burbidge.R' 'Qhat.log.R' 'Qhat.none.R' 'QhatModel.homo.normal.linear.R' 'QhatModel.homo.normal.linear.AR1.R' 'QhatModel.homo.gamma.linear.R' 'QhatModel.homo.gamma.linear.AR1.R' 'QhatModel.homo.normal.linear.AR2.R' 'QhatModel.homo.gamma.linear.AR2.R' 'QhatModel.homo.normal.linear.AR3.R' 'QhatModel.homo.gamma.linear.AR3.R' 'QhatModel.homo.skewedNormal.linear.R' 'QhatModel.homo.skewedNormal.linear.AR1.R' 'QhatModel.homo.skewedNormal.linear.AR2.R' 'QhatModel.homo.skewedNormal.linear.AR3.R' 'QhatModel.subAnnual.homo.gamma.linear.R' 'QhatModel.subAnnual.homo.gamma.linear.AR1.R' 'QhatModel.subAnnual.homo.gamma.linear.AR2.R' 'QhatModel.subAnnual.homo.gamma.linear.AR3.R' 'RhatModel.homo.normal.linear.R' 'hydroState-package.R' 'markov.annualHomogeneous.R' 'hydroState.R' 'hydroState.allModels.R' 'hydroState.subAnnual.allModels.R' 'markov.annualHomogeneous.flickering.R' 'wrapper.R'
NeedsCompilation: no
Packaged: 2025-09-16 01:32:24 UTC; twes0006
Author: Tim Peterson ORCID iD [aut, cre, cph], Thomas Westfall ORCID iD [aut]
Repository: CRAN
Date/Publication: 2025-09-22 07:40:02 UTC

hydroState: Hidden Markov modeling for hydrological state change

Description

hydroState is an R-package that identifies regime changes in streamflow time-series that is not explained by variations in precipitation.

Details

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.

Functions

Vignettes

Usage

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:

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.

Acknowledgments

The package development was funded by the Victorian Government The Department of Energy, Environment, and Climate Action (https://www.water.vic.gov.au/).

Author(s)

Maintainer: Tim Peterson tim.peterson@monash.edu (ORCID) [copyright holder]

Authors:

See Also

Useful links:


Builds hydroState model

Description

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, Q, is a function of precipitation, P: \hat{Q} = Pa_1+a_0. The default model is a variant of this base linear model with state shifts expected in the intercept, a_0, and standard deviation, std, 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

Usage

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)
)

Arguments

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 seasonal.parameters.

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: a0, a1, std. Auto-correlation terms optional: AR1, AR2, or AR3.

seasonal.parameters

character vector of one or all parameters (a0, a1, std) defined as a sinusoidal function to represent seasonal variation. Requires monthly or seasonal data. Default is empty, no seasonal parameters.

state.shift.parameters

character vector of one or all parameters (a0, a1, std, AR1, AR2, AR3) able to shift as dependent on state. Default is a0 and std.

error.distribution

character string of the distribution in the HMM error. Default is 'truc.normal'. Others include: 'normal' or 'gamma'

flickering

logical TRUE/FALSE. TRUE = allows more sensitive markov flickering between states over time. When FALSE (default), state needs to persist for at least three time steps before state shift can occur.

transition.graph

matrix given the number of states. Default is a 2-state matrix (2 by 2): matrix(TRUE,2,2).

Details

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.

Value

A built hydroState model object ready to be fitted with fit.hydroState()

Examples

# 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))


Builds all hydroState models

Description

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.

Usage

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
)

Arguments

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 build. Monthly data is required when using seasonal.parameters that assumes selected model parameters are better defined with a sinusoidal function.

data.transform

character string of method of transformation. If empty, the default builds all possible combinations of models with boxcox data transformation.

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 c('a0','std') as state shift parameters.

error.distribution

character string of the distribution in the HMM error. If empty, the default builds models with all possible combinations of error distribution: c('truc.normal', 'normal','gamma')

flickering

logical TRUE/FALSE. TRUE = allows more sensitive markov flickering between states over time. When FALSE (default), state needs to persist for at least three time steps before state shift can occur.

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): matrix(TRUE,1,1), 2-state matrix (2 by 2): matrix(TRUE,2,2), 3-state matrix (3 by 3): matrix(TRUE,3,3).

summary.table

data frame with a table summarizing all built models and corresponding reference model. From function summary(). If empty, summary table will be built automatically.

siteID

character string of site identifier.

Details

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.

Value

A list of built hydroState models with every combination of objects ready to be fitted

Examples

# 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 reliability of state predictions

Description

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.

Usage

check(model, n.samples = 1e+05)

Arguments

model

fitted hydroState model.

n.samples

integer of samples to re-sample. Default is 100000.

Details

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.

Value

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.

Examples

## Check reliability of state predictions (>5s to run)

check(model = model.annual.fitted.221201)


Fit hydroState model

Description

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.

Usage

fit.hydroState(
  model,
  pop.size.perParameter = 10,
  max.generations = 500,
  doParallel = FALSE,
  ...
)

Arguments

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

Details

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:

L_{T} = \delta P(x_{1}) + \Gamma \delta P(x_{2})...\Gamma \delta P(x_{T})1'

where:

Value

A fitted hydroState model

Examples


# 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

Description

get.AIC retrieves Akaike information criteria from a fitted hydroState model object or all models.

Usage

get.AIC(model)

Arguments

model

fitted hydroState model object.

Details

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.

Value

AIC value of a single model or a list variable of AIC values for al models

Examples

# 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)


Get pseudo residuals

Description

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.

Usage

get.residuals(model)

Arguments

model

fitted hydroState model object.

Details

get.residuals

get.residuals retrieves residuals from the fitted model and exports them as a data frame.

Value

Data frame of residuals for each time-step

Examples

# Load fitted model
data(model.annual.fitted.221201)

## Get residuals in a dataframe
get.residuals(model = model.annual.fitted.221201)



Get seasons

Description

Aggregates monthly data to 4 seasons in a year.

Usage

get.seasons(
  input.data = data.frame(year = c(), month = c(), flow = c(), precip = c())
)

Arguments

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 build. Monthly data is required when using seasonal.parameters that assumes selected model parameters are better defined with a sinusoidal function.

Details

get.seasons

This function takes sums monthly runoff and precipitation observations into 4 seasons of a year.

Value

A dataframe of seasonal observations with an additional column counting the number of months in each season.

Examples


# Load data
data(streamflow_monthly_221201)

# aggregate monthly data to seasonal
streamflow_seasonal_221201 = get.seasons(streamflow_monthly_221201)



Get states

Description

get.states uses the Viterbi algorithm to globally decode the model and estimate the most probable sequence of states.

Usage

get.states(model)

Arguments

model

fitted hydroState model object.

Details

get.states

These dataframe of results include:

Value

data frame of results to evaluate the rainfall-runoff states over time

Examples

# 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)


hydroState-internal

Description

There are several undocumented internal functions not intended for the user. These functions are not listed here but are shown within the file: hydroState-internal.Rd


Plot states or pseudo residuals over time

Description

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.

Usage

## 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
)

Arguments

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.

Details

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").

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").

Value

plots to evaluate rainfall-runoff states over time along with observations and the conditional probabilities of each state.

Examples

# 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 state names given initial year

Description

sets the state names for each time-step relative to the initial year given

Usage

setInitialYear(model, initial.year)

Arguments

model

fitted hydroState model object.

initial.year

integer with year (YYYY). Default is first year in input.data.

Details

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)

Value

A fitted hydroState model object with state names for each time-step ready for plot

Examples

# 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)



Summarize all models

Description

summary outputs a summary table of all built hydroState models and allows users to edit the reference models for calibration.

Usage

## S3 method for class 'hydroState.allModels'
summary(object, ...)

Arguments

object

hydroState.allModels object with a list of models from build.all

...

No additional input required

Details

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.

Value

A data.frame with a summary table of all models and reference models

Examples

# Show summary table of all fitted model details and reference models

all.models.ref.table = summary(all.models.annual.fitted.407211)