bdpreg: Bayesian Deming Pioda Regression for two method comparison with Rstan

Description

bdpreg is used to compare two measurement methods by means of a Bayesian regression analysis.

Usage

bdpreg(
  X,
  Y,
  ErrorRatio = 1,
  df = NULL,
  trunc = TRUE,
  heteroscedastic = c("homo", "linear"),
  slopeMu = 1,
  slopeSigma = 0.3,
  slopeTruncMin = 0.3333,
  slopeTruncMax = 10,
  interceptMu = 0,
  interceptSigma = 30,
  sigmaLambda = 0.3,
  AlphaMu = 1,
  AlphaSigma = 10,
  BetaMu = 0.1,
  BetaSigma = 0.5,
  BetaTruncMin = -1,
  BetaTruncMax = 1,
  ...
)

Value

An object of class bdpreg which contains out a stanfit object returned by rstan::sampling and standata as list of input parameters.

Arguments

X: Numeric vector of input values.
Y: Numeric vector of output values.
ErrorRatio: Deming variance ratio between reference and test method. Default = 1.
df: Degree of freedom. Must be df >= 1 (robust Cauchy regression). Default is \(N-2\), For robust regression set it to \(df < N-2\)
trunc: Boolean. Default TRUE. Use truncated slope prior for stability with extreme ErrorRatios. See slopeTruncMin.
heteroscedastic: Bayesian Deming model choice. Alternatives are: "homo" - Homoscedastic model. Default.
"linear" - Heteroscedastic with linear growth of the variance. Highly experimental model.
"exponential" - Heteroscedastic with exponential growth of the variance. Highly experimental model.
slopeMu: Slope normal Mu prior value. Default 1.
slopeSigma: Slope normal Sigma prior value. Default 0.3.
slopeTruncMin: slope normal lower truncation limit. Default 0.3333.
slopeTruncMax: slope normal higher truncation limit. Default 10.
interceptMu: Intercept normal Mu prior value. Default 0.
interceptSigma: Intercept normal Sigma prior value. Default 30.
sigmaLambda: sigma exponential prior lambda. Default 0.3.
AlphaMu: Lin. heterosc. intercept normal mu prior. Must be > 0. Default 1.
AlphaSigma: Lin. heterosc. intercept normal sigma prior. Default 10.
BetaMu: Lin. heterosc. slope normal prior. Default 0.1.
BetaSigma: Lin. heterosc. slope normal prior. Default 0.5.
BetaTruncMin: Lin. heterosc. slope normal prior truncation min. Default -1.
BetaTruncMax: Lin. heterosc. slope normal prior truncation min. Default 1.
...: Arguments passed to rstan::sampling (e.g. iter, chains)

Details

The Bayesian Deming regression can be run in a traditional fashion. In this case the error term is sampled from a \(T\) distribution with \(N-2\) degree of freedom (\(N\) sample size).

The Bayesian Deming regression can be run as a robust regression specifying a decreased \(df\) parameter. It is possible to set \(df = 1\) and perform the sampling from an extremely robust Cauchy distribution to suppress leveraged outliers. For moderate robustness a reasonably low value of \(df\) in the interval \([6;10]\) can be an appropriated choice.

ErrorRatio can be set as usual for classical Deming regression. Default is 1. Strong ErrorRatio can lead to instability in the chains that may not converge after the burn in. For this purpose the trunc parameter can be used. In this way the normal distribution for the slope gets truncated at a minimum of 0.3333 (default). The parameter slopeTruncMin can override this value.

With the parameter heteroschedastic it is possible to use an alternative regression which models the heteroscedasticity with a linear growing variance. Alpha and Beta are the intercept and the slope for the variance variation. Alpha must be > 0. Beta is usually zero if no real heteroscedasticity is detected. Alternatively Beta shows low positive values, typically below 0.5 if heteroscedasticity is successfully modeled. The CI of Beta could indeed act as a test for heteroscedasticity. According to these empiric observations, Beta is also truncated to avoid erratic behavior of the Hamiltonian sampler.

The Bayesian Deming regression is recommended in many cases where traditional and non parametric method fail. It is particularly convenient with very small data set and/or with data set with low digit precision. In fact Bayesian Deming regression has no problem with ties.

The method with linear heteroscedastic fitting can be a meaningful answer to heteroscedastic data set. The CI are much narrower and the trade off between robustness and power can find a natural solution. It must be considered as highly experimental but also highly promising method. Users are advised to carefully check the sampled output for undesirable correlation between Alpha and/or Beta vs the slope and/or intercept. A plot with pairs() highly recommended.

Stan is usually good enough that init values for the chains must not be specified. In extreme cases it is anyway possible to set init values as a list of list.

References

G. Pioda (2014) https://piodag.github.io/bd1/

Examples

Run this code

library(rstanbdp)
data(glycHem)

# Bayesian Deming Regression, for example with  df=10
fit.1 <-bdpreg(glycHem$Method1,glycHem$Method2,heteroscedastic="homo",
              df=10,chain=1,iter=1000)

# Print results
bdpPrint(fit.1,digits_summary = 4)

# Plot 2D intercepts /slopes pairs with CI and MD distance
bdpPlotBE(fit.1,cov.method="MCD",ci=0.95)

# Plot regression with CI
bdpPlot(fit.1,ci=0.95)

# Calculate response, plot histogram and CI
bdpCalcResponse(fit.1,Xval = 6)

# Extract Xhat, Yhat and Residuals
bdpExtract(fit.1)

# Plot a traceplot of the sampled chains
bdpTraceplot(fit.1)

# Plot standardized residuals
bdpPlotResiduals(fit.1)

# Plot posterior samples pairwise
bdpPairs(fit.1)

Run the code above in your browser using DataLab