Nothing
#
# Check Deming regression on a simple data set.
#
# When var(y_i) = k var(x_i) = constant there is a closed form solution,
# see
#
closed <- function(x, y, k, w) {
# if w is missing assume a vector of ones.
n <- length(x)
if (missing(w)) wt <- rep(1/n, length(x))
else wt <- (1/w) / sum(1/w)
xbar <- sum(wt * x)
ybar <- sum(wt * y)
xmat <- cbind(y-ybar, x-xbar)
ss <- t(xmat) %*% diag(wt) %*% xmat
temp <- (ss[1,1] - k*ss[2,2])
beta <- (temp + sqrt(temp^2 + 4*k*ss[1,2]^2))/ (2 * ss[1,2])
beta
}
require(deming)
aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
x <- 1:10
y <- 1:10 *2.3 + c(0, -1, 2, -3, 4, -5, 6, -7, 8, -5)/2
# The deming routine uses optimize, which defaults to a tolerance of
# .Machine$double.eps^.25, less than the default for all.equal.
tol <- .Machine$double.eps^.25
true <- closed(x,y, 1)
dfit1 <- deming(y ~ x)
aeq(coef(dfit1)[2], true, tol=tol)
dfit2 <- deming(x ~ y)
aeq(1/true, coef(dfit2)[2], tol=tol)
# verify an alternate form: the Deming angle is that rotation which gives
# a least-squares regression coef of 0
rotate <- function(theta, x, y) {
data.frame(x = x*cos(theta) + y* sin(theta),
y = y*cos(theta) - x* sin(theta))
}
lfit <- lm(y ~ x, data= rotate(atan(true), x, y))
aeq(coef(lfit)[2], 0)
# Regression through the origin
dfit3 <- deming(y ~ x-1)
ofun <- function(slope) {
theta <- atan(slope)
newdata <- rotate(theta, x, y)
mean(newdata$y^2)
}
ofit <- optimize(ofun, lower=.1, upper=3)
aeq(dfit3$coef[2], ofit$minimum, tol=1e-5)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.