plm: Power-law model with variance that varies with stage.
In sor16/RCmodels: Bayesian Discharge Rating Curves
plm R Documentation
Power-law model with variance that varies with stage.
Description
plm is used to fit a discharge rating curve for paired measurements of stage and discharge using a power-law model with variance that varies with stage as described in Hrafnkelsson et al. (2022). See "Details" for a more elaborate description of the model.
Usage
plm(
formula,
data,
c_param = NULL,
h_max = NULL,
parallel = TRUE,
num_cores = NULL,
forcepoint = rep(FALSE, nrow(data)),
verbose = TRUE
)
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
A data.frame containing the variables specified in formula.
c_param
The largest stage value for which there is zero discharge. If NULL, it is treated as unknown in the model and inferred from the data.
h_max
The 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
A logical value indicating whether to run the MCMC in parallel or not. Defaults to TRUE.
num_cores
An 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
A 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.
verbose
A logical value indicating whether to print progress and diagnostic information. If 'TRUE', the function will print messages as it runs. If 'FALSE', the function will run silently. Default is 'TRUE'.
Details
The power-law model, which is commonly used in hydraulic practice, is of the form
Q=a(h-c)^{b}
where Q is discharge, h is stage and a, b and c are unknown constants.
The power-law model is here inferred by using a Bayesian hierarchical model. The model is on a logarithmic scale
\log(Q_i) = \log(a) + b \log(h_i - c) + \varepsilon_i, i = 1,...,n
where \varepsilon_i follows a normal distribution with mean zero and variance \sigma_\varepsilon(h_i)^2 that varies with stage. The error variance, \sigma_\varepsilon^2(h), of the log-discharge data is modeled as an exponential of a B-spline curve, that is, a linear combination of six B-spline basis functions that are defined over the range of the stage observations. 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
plm returns an object of class "plm". An object of class "plm" is a list containing the following components:
rating_curveA data frame with 2.5%, 50% and 97.5% percentiles of the posterior predictive distribution of the rating curve.
rating_curve_meanA data frame with 2.5%, 50% and 97.5% percentiles of the posterior distribution of the mean of the rating curve. Additionally contains columns with r_hat and the effective number of samples for each parameter as defined in Gelman et al. (2013).
param_summaryA data frame with 2.5%, 50% and 97.5% percentiles of the posterior distribution of latent- and hyperparameters.
sigma_eps_summaryA data frame with 2.5%, 50% and 97.5% percentiles of the posterior of \sigma_{\varepsilon}.
posterior_log_likelihood_summaryA data frame with 2.5%, 50% and 97.5% percentiles of the posterior log-likelihood values.
rating_curve_posteriorA matrix containing the full thinned posterior samples of the posterior predictive distribution of the rating curve (excluding burn-in).
rating_curve_mean_posteriorA matrix containing the full thinned posterior samples of the posterior distribution of the mean of the rating curve (excluding burn-in).
a_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of a.
b_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of b.
c_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of c.
sigma_eps_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of \sigma_{\varepsilon}.
eta_1_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of \eta_1.
eta_2_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of \eta_2.
eta_3_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of \eta_3.
eta_4_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of \eta_4.
eta_5_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of \eta_5.
eta_6_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of \eta_6.
posterior_log_likelihoodA numeric vector containing the full thinned posterior log-likelihood values, excluding burn-in samples.
D_hatA statistic defined as -2 times the log-likelihood evaluated at the median value of the parameters.
effective_num_param_DICThe effective number of parameters, which is calculated as median(-2*posterior_log_likelihood) minus D_hat.
DICThe Deviance Information Criterion for the model, calculated as D_hat plus 2*effective_num_parameters_DIC.
lppdThe log pointwise predictive density of the observed data under the model.
WAICThe Widely Applicable Information Criterion for the model, defined as -2*( lppd - effective_num_param_WAIC ).
WAIC_iThe pointwise WAIC values, where WAIC := sum(WAIC_i).
effective_num_param_WAICThe effective number of parameters, which is calculated by summing up the posterior variance of the log predictive density for each data point.
autocorrelationA data frame with the autocorrelation of each parameter for different lags.
acceptance_rateThe proportion of accepted samples in the thinned MCMC chain (excluding burn-in).
formulaAn object of type "formula" provided by the user.
dataThe data provided by the user, ordered by stage.
run_infoThe 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. doi: https://doi.org/10.1201/b16018
Hrafnkelsson, B., Sigurdarson, H., Rögnvaldsson, S., Jansson, A. Ö., Vias, R. D., and Gardarsson, S. M. (2022). Generalization of the power-law rating curve using hydrodynamic theory and Bayesian hierarchical modeling, Environmetrics, 33(2):e2711. doi: https://doi.org/10.1002/env.2711
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. doi: https://doi.org/10.1111/1467-9868.00353
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.plm for summaries, predict.plm for prediction. It is also useful to look at spread_draws and plot.plm to help visualize the full posterior distributions.
Examples
data(spanga)
set.seed(1)
plm.fit <- plm(formula=Q~W,data=spanga,num_cores=2)
summary(plm.fit)
sor16/RCmodels documentation built on Sept. 13, 2024, 4:51 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| plm | R Documentation |
Power-law model with variance that varies with stage.
Description
plm is used to fit a discharge rating curve for paired measurements of stage and discharge using a power-law model with variance that varies with stage as described in Hrafnkelsson et al. (2022). See "Details" for a more elaborate description of the model.
Usage
plm(
formula,
data,
c_param = NULL,
h_max = NULL,
parallel = TRUE,
num_cores = NULL,
forcepoint = rep(FALSE, nrow(data)),
verbose = TRUE
)
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 |
data |
A data.frame containing the variables specified in formula. |
c_param |
The largest stage value for which there is zero discharge. If NULL, it is treated as unknown in the model and inferred from the data. |
h_max |
The 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 |
A logical value indicating whether to run the MCMC in parallel or not. Defaults to TRUE. |
num_cores |
An 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 |
A logical vector of the same length as the number of rows in data. If an element at index |
verbose |
A logical value indicating whether to print progress and diagnostic information. If 'TRUE', the function will print messages as it runs. If 'FALSE', the function will run silently. Default is 'TRUE'. |
Details
The power-law model, which is commonly used in hydraulic practice, is of the form
Q=a(h-c)^{b}
where Q is discharge, h is stage and a, b and c are unknown constants.
The power-law model is here inferred by using a Bayesian hierarchical model. The model is on a logarithmic scale
\log(Q_i) = \log(a) + b \log(h_i - c) + \varepsilon_i, i = 1,...,n
where \varepsilon_i follows a normal distribution with mean zero and variance \sigma_\varepsilon(h_i)^2 that varies with stage. The error variance, \sigma_\varepsilon^2(h), of the log-discharge data is modeled as an exponential of a B-spline curve, that is, a linear combination of six B-spline basis functions that are defined over the range of the stage observations. 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
plm returns an object of class "plm". An object of class "plm" is a list containing the following components:
rating_curveA data frame with 2.5%, 50% and 97.5% percentiles of the posterior predictive distribution of the rating curve.
rating_curve_meanA data frame with 2.5%, 50% and 97.5% percentiles of the posterior distribution of the mean of the rating curve. Additionally contains columns with r_hat and the effective number of samples for each parameter as defined in Gelman et al. (2013).
param_summaryA data frame with 2.5%, 50% and 97.5% percentiles of the posterior distribution of latent- and hyperparameters.
sigma_eps_summaryA data frame with 2.5%, 50% and 97.5% percentiles of the posterior of
\sigma_{\varepsilon}.posterior_log_likelihood_summaryA data frame with 2.5%, 50% and 97.5% percentiles of the posterior log-likelihood values.
rating_curve_posteriorA matrix containing the full thinned posterior samples of the posterior predictive distribution of the rating curve (excluding burn-in).
rating_curve_mean_posteriorA matrix containing the full thinned posterior samples of the posterior distribution of the mean of the rating curve (excluding burn-in).
a_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of
a.b_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of
b.c_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of
c.sigma_eps_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of
\sigma_{\varepsilon}.eta_1_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of
\eta_1.eta_2_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of
\eta_2.eta_3_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of
\eta_3.eta_4_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of
\eta_4.eta_5_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of
\eta_5.eta_6_posteriorA numeric vector containing the full thinned posterior samples of the posterior distribution of
\eta_6.posterior_log_likelihoodA numeric vector containing the full thinned posterior log-likelihood values, excluding burn-in samples.
D_hatA statistic defined as -2 times the log-likelihood evaluated at the median value of the parameters.
effective_num_param_DICThe effective number of parameters, which is calculated as median(-2*posterior_log_likelihood) minus D_hat.
DICThe Deviance Information Criterion for the model, calculated as D_hat plus 2*effective_num_parameters_DIC.
lppdThe log pointwise predictive density of the observed data under the model.
WAICThe Widely Applicable Information Criterion for the model, defined as -2*( lppd - effective_num_param_WAIC ).
WAIC_iThe pointwise WAIC values, where WAIC := sum(WAIC_i).
effective_num_param_WAICThe effective number of parameters, which is calculated by summing up the posterior variance of the log predictive density for each data point.
autocorrelationA data frame with the autocorrelation of each parameter for different lags.
acceptance_rateThe proportion of accepted samples in the thinned MCMC chain (excluding burn-in).
formulaAn object of type "formula" provided by the user.
dataThe data provided by the user, ordered by stage.
run_infoThe 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. doi: https://doi.org/10.1201/b16018
Hrafnkelsson, B., Sigurdarson, H., Rögnvaldsson, S., Jansson, A. Ö., Vias, R. D., and Gardarsson, S. M. (2022). Generalization of the power-law rating curve using hydrodynamic theory and Bayesian hierarchical modeling, Environmetrics, 33(2):e2711. doi: https://doi.org/10.1002/env.2711
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. doi: https://doi.org/10.1111/1467-9868.00353
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.plm for summaries, predict.plm for prediction. It is also useful to look at spread_draws and plot.plm to help visualize the full posterior distributions.
Examples
data(spanga)
set.seed(1)
plm.fit <- plm(formula=Q~W,data=spanga,num_cores=2)
summary(plm.fit)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.