在没有截距的情况下执行戴明回归
perform Deming regression without intercept
我想执行戴明回归(或 X 和 Y 变量都具有不确定性的任何等效回归方法,例如 York 回归)。
在我的申请中,我有很好的科学依据故意将截距设置为零。但是,我找不到将它设置为零的方法,无论是在 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,所以我无法输入模型矩阵或修改公式。
您知道如何实现吗?谢谢。
删除拦截时函数中存在错误。我要举报了
修复起来很简单,只需要在原来的函数中改两行即可。
打印不正确,但可以解释输出。
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 .- 修复更改可能产生的其他错误。
在这种情况下,删除输出的class以避免打印错误。
错误报告:
Terry Therneau(deming 的作者)上传了一个新版本到 CRAN,解决了这个问题。
我想执行戴明回归(或 X 和 Y 变量都具有不确定性的任何等效回归方法,例如 York 回归)。
在我的申请中,我有很好的科学依据故意将截距设置为零。但是,我找不到将它设置为零的方法,无论是在 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,所以我无法输入模型矩阵或修改公式。
您知道如何实现吗?谢谢。
删除拦截时函数中存在错误。我要举报了
修复起来很简单,只需要在原来的函数中改两行即可。 打印不正确,但可以解释输出。
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 .- 修复更改可能产生的其他错误。
在这种情况下,删除输出的class以避免打印错误。
错误报告:
Terry Therneau(deming 的作者)上传了一个新版本到 CRAN,解决了这个问题。