gplm0: Generalized power-law model with a constant variance

In bdrc: Bayesian Discharge Rating Curves

Generalized power-law model with a constant variance

Description

gplm0 is used to fit a discharge rating curve for paired measurements of stage and discharge using a generalized power-law model with a constant variance as described in Hrafnkelsson et al. (2022). See "Details" for a more elaborate description of the model.

Usage

gplm0(
  formula,
  data,
  c_param = NULL,
  h_max = NULL,
  parallel = TRUE,
  num_cores = NULL,
  forcepoint = rep(FALSE, nrow(data))
)

Arguments

formula	an object of class "formula", with discharge column name as response and stage column name as a covariate, i.e. of the form y~x where y is discharge in m^3/s and x is stage in m (it is very important that the data is in the correct units).
data	data.frame containing the variables specified in formula.
c_param	stage for which there is zero discharge. If NULL, it is treated as unknown in the model and inferred from the data.
h_max	maximum stage to which the rating curve should extrapolate to. If NULL, the maximum stage value in the data is selected as an upper bound.
parallel	logical value indicating whether to run the MCMC in parallel or not. Defaults to TRUE.
num_cores	integer between 1 and 4 (number of MCMC chains) indicating how many cores to use. Only used if parallel=TRUE. If NULL, the number of cores available on the device is detected automatically.
forcepoint	logical vector of the same length as the number of rows in data. If an element at index i is TRUE it indicates that the rating curve should be forced through the i-th measurement. Use with care, as this will strongly influence the resulting rating curve.

Details

The generalized power-law model is of the form

Q=a(h-c)^{f(h)}

where Q is discharge, h is stage, a and c are unknown constants and f is a function of h referred to as the generalized power-law exponent.

The generalized power-law model is here inferred by using a Bayesian hierarchical model. The function f is modeled at the latent level as a fixed constant $b$ plus a continuous stochastic process,\beta(h), which is assumed to be twice differentiable. The model is on a logarithmic scale

\log(Q_i) = \log(a) + (b + \beta(h_i)) \log(h_i - c) + \varepsilon, i = 1,...,n

where \varepsilon follows a normal distribution with mean zero and variance \sigma_\varepsilon^2, independent of stage. The stochastic process \beta(h) is assumed a priori to be a Gaussian process governed by a Matern covariance function with smoothness parameter \nu = 2.5. An efficient posterior simulation is achieved by sampling from the joint posterior density of the hyperparameters of the model, and then sampling from the density of the latent parameters conditional on the hyperparameters.

Bayesian inference is based on the posterior density and summary statistics such as the posterior mean and 95% posterior intervals are based on the posterior density. Analytical formulas for these summary statistics are intractable in most cases and thus they are computed by generating samples from the posterior density using a Markov chain Monte Carlo simulation.

Value

gplm0 returns an object of class "gplm0". An object of class "gplm0" is a list containing the following components:

rating_curve	a data frame with 2.5%, 50% and 97.5% percentiles of the posterior predictive distribution of the rating curve.
rating_curve_mean	a data frame with 2.5%, 50% and 97.5% percentiles of the posterior distribution of the mean of the rating curve.
param_summary	a data frame with 2.5%, 50% and 97.5% percentiles of the posterior distribution of latent- and hyperparameters. Additionally contains columns with r_hat and the effective number of samples for each parameter as defined in Gelman et al. (2013).
beta_summary	a data frame with 2.5%, 50% and 97.5% percentiles of the posterior distribution of \beta.
Deviance_summary	a data frame with 2.5%, 50% and 97.5% percentiles of the posterior distribution of the deviance.
rating_curve_posterior	a matrix containing the full thinned posterior samples of the posterior predictive distribution of the rating curve (excluding burn-in).
rating_curve_mean_posterior	a matrix containing the full thinned posterior samples of the posterior distribution of the mean of the rating curve (excluding burn-in).
a_posterior	a numeric vector containing the full thinned posterior samples of the posterior distribution of a.
b_posterior	a numeric vector containing the full thinned posterior samples of the posterior distribution of b.
c_posterior	a numeric vector containing the full thinned posterior samples of the posterior distribution of c.
sigma_eps_posterior	a numeric vector containing the full thinned posterior samples of the posterior distribution of \sigma_{\varepsilon}.
sigma_beta_posterior	a numeric vector containing the full thinned posterior samples of the posterior distribution of \sigma_{\beta}.
phi_beta_posterior	a numeric vector containing the full thinned posterior samples of the posterior distribution of \phi_{\beta}.
sigma_eta_posterior	a numeric vector containing the full thinned posterior samples of the posterior distribution of \sigma_{\eta}.
beta_posterior	a numeric vector containing the full thinned posterior samples of the posterior distribution of \beta.
Deviance_posterior	a numeric vector containing the full thinned posterior samples of the posterior distribution of the deviance excluding burn-in samples.
D_hat	deviance at the median value of the parameters.
effective_num_param_DIC	effective number of parameters, which is calculated as median(Deviance_posterior) minus D_hat.
DIC	Deviance Information Criterion for the model, calculated as D_hat plus 2*effective_num_parameters_DIC.
lppd	log pointwise predictive probability of the observed data under the model
effective_num_param_WAIC	effective number of parameters, which is calculated by summing up the posterior variance of the log predictive density for each data point.
WAIC	Watanabe-Akaike information criterion for the model, defined as -2*( lppd - effective_num_param_WAIC ).
autocorrelation	a data frame with the autocorrelation of each parameter for different lags.
acceptance_rate	proportion of accepted samples in the thinned MCMC chain (excluding burn-in).
formula	object of type "formula" provided by the user.
data	data provided by the user, ordered by stage.
run_info	information about the input arguments and the specific parameters used in the MCMC chain.

References

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013). Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science. Taylor & Francis.

Hrafnkelsson, B., Sigurdarson, H., and Gardarsson, S. M. (2022). Generalization of the power-law rating curve using hydrodynamic theory and Bayesian hierarchical modeling, Environmetrics, 33(2):e2711.

Spiegelhalter, D., Best, N., Carlin, B., Van Der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 64(4), 583–639.

Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. J. Mach. Learn. Res. 11, 3571–3594.

See Also

summary.gplm0 for summaries, predict.gplm0 for prediction. It is also useful to look at spread_draws and plot.gplm0 to help visualize the full posterior distributions.

Examples

data(krokfors)
set.seed(1)
gplm0.fit <- gplm0(formula=Q~W,data=krokfors,num_cores=2)
summary(gplm0.fit)

bdrc documentation built on March 31, 2023, 11:41 p.m.