leiv {leiv}

Bivariate Linear Errors-In-Variables Estimation

Description

Generates a linear errors-in-variables object.

Usage

leiv(formula, data, subset, prior = NULL,
      n = NULL, cor = NULL, sdRatio = NULL, xMean = 0, yMean = 0,
      probIntCalc = FALSE, level = 0.95, subdivisions = 100,
      rel.tol = .Machine$double.eps^0.25, abs.tol = 0.1*rel.tol, ...)

## S4 method for signature 'leiv'
print(x, digits = max(3, getOption("digits") - 3), ...)
## S4 method for signature 'leiv,missing'
plot(x, plotType = "density", xlim = NULL, ylim = NULL,
     xlab = NULL, ylab = NULL, col = NULL, lwd = NULL, ...)

Arguments

formula

an optional object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. The details of model specification are given in the ‘Details’ section of the documentation for lm. An intercept is always included and integrated out as a nuisance parameter: y ~ x, y ~ 0 + x, and y ~ x - 1 are equivalent. If not provided, the sufficient statistics n, cor, and sdRatio must be provided.

data

an optional data frame (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which leiv is called.

subset

an optional vector specifying a subset of observations to be used in the fitting process.

prior

an optional object of class leiv to use as the prior density of the scale invariant slope; otherwise the rotationally invariant Cauchy density is used.

n

an optional sample size (if formula is missing).

cor, sdRatio

optional sample correlation cor(x,y) and ratio sd(y)/sd(x) (if formula is missing).

xMean, yMean

optional sample means mean(x) and mean(y) (if formula is missing).

probIntCalc

logical; if TRUE returns the shortest (100*level)% probability intervals; if FALSE (the default) no probability intervals are returned.

level

the probability level requested (if probIntCalc = TRUE).

subdivisions

the maximum number of subintervals (see integrate).

rel.tol

the relative accuracy requested (see integrate).

abs.tol

the absolute accuracy requested (see integrate).

x

a leiv object.

digits

controls formating of numeric objects.

plotType

specifies the type of plot; if plotType = "density" (the default) then the posterior density of the slope is plotted; if plotType = "scatter" then a scatter plot with the fitted line.

xlim, ylim

x limits c(x1,x2) and y limits c(y1,y2) of the plot.

xlab, ylab

labels for the x and y axes of the plot.

col, lwd

color and width of plotted lines.

...

additional argument(s) for generic methods.

Details

Use leiv to estimate the slope and intercept of a bivariate linear relationship when both variables are observed with error. The method is exact when the true values and the errors are normally distributed. The posterior density depends on the data only through the correlation coefficient and ratio of standard deviations; it is invariant to interchange and scaling of the coordinates.

Value

leiv returns an object of class "leiv" with the following components:

slope

the (posterior median) slope estimate.

intercept

the (maximum likelihood) intercept estimate.

slopeInt

the shortest (100*level)% probability interval of the slope.

interceptInt

the shortest (100*level)% probability interval of the intercept.

density

the posterior probability density function.

n

the number of (x,y) pairs.

cor

the sample correlation cor(x,y).

sdRatio

the ratio sd(y)/sd(x).

xMean

the sample mean mean(x).

yMean

the sample mean mean(y).

call

the matched call.

probIntCalc

the logical probability interval request.

level

the probability level of the probability interval.

x

the x data.

y

the y data.

Note

Numerical integration is used to normalize the posterior density. When the data is nearly linear, normalization using the default tolerance parameters may fail. Specifying abs.tol = 1e-6 (or smaller) may help, but expect a longer run time. In general, rel.tol cannot be less than max(50*.Machine$double.eps, 0.5e-28) if abs.tol <= 0. In addition, when using a sharply peaked leiv object as a prior density, normalization may fail. In this case, an alternative is to first fit using the default Cauchy prior, then multiply by the appropriate ratio of prior densities and tackle the normalization outside of the leiv environment.

Author(s)

David Leonard

References

Leonard, David. (2011). “Estimating a Bivariate Linear Relationship.” Bayesian Analysis, 6:727-754. DOI:10.1214/11-BA627.

Zellner, Arnold. (1971). An Introduction to Bayesian Inference in Econometrics, Chapter 5. John Wiley & Sons.

See Also

lm for formula syntax; integrate for control parameters.

Examples

## generate artificial data 
set.seed(1123)
n <- 20
X <- rnorm(n, mean=5, sd=4) # true x
x <- X + rnorm(n, mean=0, sd=5) # observed x
Y <- 2 + X # true y
y <- Y + rnorm(n, mean=0, sd=3) # observed y

## fit with default options
fit <- leiv(y ~ x)
print(fit)
plot(fit) # density plot
dev.new()
plot(fit,plotType="scatter")

## calculate a density to use as an informative prior density of
## the scale invariant slope in a subsequent fit
fit0 <- leiv(n=10, cor=0.5, sdRatio=1.0)
print(fit0)

## refit the data using the informative prior density
fit1 <- leiv(y ~ x, prior=fit0, abs.tol=1e-6)
print(fit1)