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)
}

より具体的には、lm クラス (lm$qr) に含まれるクラス "qr" の QR 分解で使用される基本パッケージから qr.R を使用するのと同じ結果を得ようとしています。ベース関数のコードは次のとおりです。

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 行列に対して biglm の lm と同じ結果が得られるようにしたいと考えています。

4

1 に答える 1

2

2 つの R 行列の違いは、biglmは明らかに R の対角要素がすべて正になるように回転を実行するのに対し、lm (または実際にはそれが呼び出すルーチン) はそのような制約を課さないことです。（どちらの戦略にも数値的な利点はないはずなので、その違いは慣例の 1 つにすぎません。）

追加の制約を自分で課すことにより、 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 に答える