r - 执行无截距的戴明回归

Question

我想执行戴明回归(或任何 X 和 Y 变量都具有不确定性的回归方法的等效方法，例如约克回归)。

在我的应用程序中，我有很好的科学依据来故意将截距设置为零。但是，我找不到将其设置为零的方法，无论是在 R 包 deming 中，当我在公式中使用 -1 时，它都会出错:

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

在其他包中(例如 mcr::mcreg 或 IsoplotR::york 或 MethComp::Deming)，输入是两个向量x和y，所以我无法输入模型矩阵或修改公式。

您知道如何实现这一目标吗？谢谢。

Answer 1

当您删除拦截时，函数中存在错误。我要举报。

修复起来很容易，只需更改原始函数中的 2 行即可。 打印无法正常工作，但可以解释输出。

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
...

为了调整我所执行的功能，请执行以下步骤:

1 .- 加载包的内部函数。

deming.fit1 <- getAnywhere(deming.fit1)[[2]][[1]]
deming.fit2 <- getAnywhere(deming.fit2)[[2]][[1]]

2 .- 定位问题并解决，通过示例逐步执行该功能。

Y <- as.matrix(model.response(mf, type = "numeric"))
...
xx <- X[, ncol(X), drop = FALSE]

3 .- 修复因更改而可能产生的其他错误。

在这种情况下，请删除输出的类以避免打印错误。

错误报告:

Terry Therneau(deming的作者)将新版本上传到CRAN，解决了这个问题。