bdpreg {rstanbdp}

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

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.

Value

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

References

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

Examples

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)