3

我正在尝试从 biglm 中使用的 QR 分解中恢复 R 矩阵。为此，我使用了 vcov.biglm 中的一部分代码并将其放入如下函数中：

qr.R.biglm <- function (object, ...) {
    # Return the qr.R matrix from a biglm object
    object$qr <- .Call("singcheckQR", object$qr)
    p <- length(object$qr$D)
    R <- diag(p)
    R[row(R) > col(R)] <- object$qr$rbar
    R <- t(R)
    R <- sqrt(object$qr$D) * R
    dimnames(R) <- list(object$names, object$names)
    return(R)
}

更具体地说，我试图从基本包中获得与使用 qr.R 相同的结果，该包用于类“qr”的 QR 分解，例如 lm 类 (lm$qr) 中包含的那些。基本函数的代码如下：

qr.R <- function (qr, complete = FALSE) {
    if (!is.qr(qr)) 
        stop("argument is not a QR decomposition")
    R <- qr$qr
    if (!complete) 
        R <- R[seq.int(min(dim(R))), , drop = FALSE]
    R[row(R) > col(R)] <- 0
    R
}

除了符号之外，我设法为样本回归获得了相同的结果。

x <- as.data.frame(matrix(rnorm(100 * 10), 100, 10))
y <- seq.int(1, 100)
fit.lm <- lm("y ~ .", data =  cbind(y, x))
R.lm <- qr.R(fit.lm$qr)

library(biglm)
fmla <- as.formula(paste("y ~ ", paste(colnames(x), collapse = "+")))
fit.biglm <- biglm(fmla, data = cbind(y, x))
R.biglm <- qr.R.biglm(fit.biglm)

比较两者，很明显绝对值匹配，但符号不匹配。

mean(abs(R.lm) - abs(R.biglm) < 1e-6)
[1] 1
mean(R.lm - R.biglm < 1e-6)
[1] 0.9338843

我不太明白为什么会这样。我希望能够从 biglm 获得与 lm 相同的 R 矩阵结果。

4

1 回答 1

2

两个 R 矩阵之间的区别在于biglm显然执行它的旋转使得 R 的对角线元素都是正数，而lm（或者，实际上，它调用的例程）没有施加这样的约束。（一种策略或另一种策略应该没有数字优势，所以区别只是约定之一，AFAIKT。）

您可以通过自己施加额外的约束来使lm的结果与biglm的结果相同。我会使用一个反射矩阵，将列乘以 1 或 -1，这样对角线元素最终都是正数：

## Apply the necessary reflections
R.lm2 <- diag(sign(diag(R.lm))) %*% R.lm

## Show that they did the job
mean(R.lm2 - R.biglm < 1e-6)
# [1] 1
于 2012-11-06T13:57:25.497 回答