I would like to perform Deming regression (or any equivalent of a regression method with uncertainties in both X and Y variables, such as York regression).
In my application, I have a very good scientific justification to deliberately set the intercept as zero. However, I can't find way of setting it to zero, either in R package deming, it gives an error when I use -1 in the formula:
df=data.frame(x=rnorm(10), y=rnorm(10), sx=runif(10), sy=runif(10))
library(deming)
deming(y~x-1, df, xstd=sy, ystd=sy)
Error in lm.wfit(x, y, wt/ystd^2) : 'x' must be a matrix
In other packages (like mcr::mcreg or IsoplotR::york or MethComp::Deming), the input are two vectors x and y, so there is no way I can input a model matrix or modify the formula.
Do you have any idea on how to achieve that? Thanks.
There is a bug in the function when you remove the intercept. I'm going to report it.
It is easy to fix, you just have to change 2 lines in the original function. The print does not work correctly, but it is possible to interpret the output.
deming.aux <-
function (formula, data, subset, weights, na.action, cv = FALSE,
xstd, ystd, stdpat, conf = 0.95, jackknife = TRUE, dfbeta = FALSE,
x = FALSE, y = FALSE, model = TRUE)
{
deming.fit1 <- getAnywhere(deming.fit1)[[2]][[1]]
deming.fit2 <- getAnywhere(deming.fit2)[[2]][[1]]
Call <- match.call()
indx <- match(c("formula", "data", "weights", "subset", "na.action", "xstd", "ystd"), names(Call), nomatch = 0)
if (indx[1] == 0)
stop("A formula argument is required")
temp <- Call[c(1, indx)]
temp[[1]] <- as.name("model.frame")
mf <- eval(temp, parent.frame())
Terms <- terms(mf)
n <- nrow(mf)
if (n < 3)
stop("less than 3 non-missing observations in the data set")
xstd <- model.extract(mf, "xstd")
ystd <- model.extract(mf, "ystd")
Y <- as.matrix(model.response(mf, type = "numeric"))
if (is.null(Y))
stop("a response variable is required")
wt <- model.weights(mf)
if (length(wt) == 0)
wt <- rep(1, n)
usepattern <- FALSE
if (is.null(xstd)) {
if (!is.null(ystd))
stop("both of xstd and ystd must be given, or neither")
if (missing(stdpat)) {
if (cv)
stdpat <- c(0, 1, 0, 1)
else stdpat <- c(1, 0, 1, 0)
}
else {
if (any(stdpat < 0) || all(stdpat[1:2] == 0) || all(stdpat[3:4] ==
0))
stop("invalid stdpat argument")
}
if (stdpat[2] > 0 || stdpat[4] > 0)
usepattern <- TRUE
else {
xstd <- rep(stdpat[1], n)
ystd <- rep(stdpat[3], n)
}
}
else {
if (is.null(ystd))
stop("both of xstd and ystd must be given, or neither")
if (!is.numeric(xstd))
stop("xstd must be numeric")
if (!is.numeric(ystd))
stop("ystd must be numeric")
if (any(xstd <= 0))
stop("xstd must be positive")
if (any(ystd <= 0))
stop("ystd must be positive")
}
if (conf < 0 || conf >= 1)
stop("invalid confidence level")
if (!is.logical(dfbeta))
stop("dfbeta must be TRUE or FALSE")
if (dfbeta & !jackknife)
stop("the dfbeta option only applies if jackknife=TRUE")
X <- model.matrix(Terms, mf)
if (ncol(X) != (1 + attr(Terms, "intercept")))
stop("Deming regression requires a single predictor variable")
xx <- X[, ncol(X), drop = FALSE]
if (!usepattern)
fit <- deming.fit1(xx, Y, wt, xstd, ystd, intercept = attr(Terms, "intercept"))
else fit <- deming.fit2(xx, Y, wt, stdpat, intercept = attr(Terms, "intercept"))
yhat <- fit$coefficients[1] + fit$coefficients[2] * xx
fit$residuals <- Y - yhat
if (x)
fit$x <- X
if (y)
fit$y <- Y
if (model)
fit$model <- mf
na.action <- attr(mf, "na.action")
if (length(na.action))
fit$na.action <- na.action
fit$n <- length(Y)
fit$terms <- Terms
fit$call <- Call
fit
}
deming.aux(y ~ x + 0, df, xstd=sy, ystd=sy)
$`coefficients`
[1] 0.000000 4.324481
$se
[1] 0.2872988 0.7163073
$sigma
[1] 2.516912
$residuals
[,1]
1 9.19499513
2 2.13037957
3 3.00064886
4 2.16751905
5 0.00168729
6 4.75834265
7 3.44108236
8 6.40028085
9 6.63531039
10 -1.48624851
$model
y x (xstd) (ystd)
1 2.1459817 -1.6300251 0.48826221 0.48826221
2 1.3163362 -0.1882407 0.46002166 0.46002166
3 1.5263967 -0.3409084 0.55771660 0.55771660
4 -0.9078000 -0.7111417 0.81145673 0.81145673
5 -1.6768719 -0.3881527 0.01563191 0.01563191
6 -0.6114545 -1.2417205 0.41675425 0.41675425
7 -0.7783790 -0.9757150 0.82498713 0.82498713
8 1.1240046 -1.2200946 0.84072712 0.84072712
9 -0.3091330 -1.6058442 0.35926078 0.35926078
10 0.7215432 0.5105333 0.23674788 0.23674788
$n
[1] 10
$terms
y ~ x + 0
...
To adapt the function I have performed, these steps:
1 .- Load the internal functions of the package.
deming.fit1 <- getAnywhere(deming.fit1)[[2]][[1]]
deming.fit2 <- getAnywhere(deming.fit2)[[2]][[1]]
2 .- Locate the problem and solve it, executing the function step by step with an example.
Y <- as.matrix(model.response(mf, type = "numeric"))
...
xx <- X[, ncol(X), drop = FALSE]
3 .- Fix other possible error generated by the changes.
In this case, delete the class of the output to avoid an error in the print.
Bug report:
Terry Therneau (the author of deming) uploaded a new version to CRAN, with this problem solved.