perform Deming regression without intercept

Question

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.

Answer 1

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.