我试图从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类中包含的qr (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矩阵的结果。
发布于 2012-11-06 13:57:25
两个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] 1https://stackoverflow.com/questions/13248960
复制相似问题