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)