Я пытаюсь восстановить матрицу R из разложения QR, используемого в biglm. Для этого я использую часть кода в 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
Я не могу понять, почему это так. Я хотел бы иметь возможность получить тот же результат для матрицы R, что и lm из biglm.
lmнаbiglmв каком-то существующем коде, так как у меня есть большой набор данных, для которого я хочу выполнить регрессию, и я не смогу использоватьlm. Моя программа использует эту матрицу R дальше по линии, и было бы неплохо, если бы я мог просто извлечь ту же матрицу, а не придумывать новую формулу. Если до этого дойдет, так тому и быть :-) просто решил сначала поспрашивать. 06.11.2012