次に、応答の精度に対応する重みをモデルに追加したいと思います。したがって、応答ベクトルごとに、異なる重みベクトルも必要になります。ただし、lm
では重みを行列ではなくベクトルとしてのみ入力できます。そのため、重みをバッチで入力できず、lm
応答ごとに個別に実行する必要がありました。この方法では、計算がはるかに遅くなります。
複数の LHS による線形モデルの適合で説明されているように、"mlm" の効率は、すべての LHS 応答に対して共有モデル行列を必要とします。ただし、重み付き回帰では、モデル マトリックスを再利用することはできません。重みのセットが異なる場合は、応答y
とモデル マトリックスの両方X
を再スケーリングする必要があります。R: lm()weights
の結果は、引数を使用する場合と、手動で再重み付けされたデータを使用する場合とで異なり、重み付けされた回帰がどのように機能するかを確認してください。
lm
繰り返し呼び出すことなく、これらのタイプのモデルをバッチ モードで実行する方法はありますか?
それはあなたが望むものに依存します。完全な が必要な場合は、毎回lmObject
呼び出す必要があります。lm
係数のみが必要な場合は、 を使用できます.lm.fit
。上記の 2 番目のリンクは の使用を示していますがlm.fit
、 の使用.lm.fit
はほとんど同じです。単純なテンプレートは次のようになります。
## weighted mlm, by specifying matrix directly
## `xmat`: non-weighted model matrix, manually created from `model.matrix`
## `ymat`: non-weighted response matrix
## `wmat`: matrix of weights
## all matrices must have the same number of rows (not checked)
## `ymat` and `wmat` must have the same number of columns (not checked)
## no `NA` values in any where is allowed (not checked)
## all elements of `wmat` must be strictly positive (not checked)
wmlm <- function (xmat, ymat, wmat) {
N <- ncol(ymat)
wmlmList <- vector("list", length = N)
for (j in 1:N) {
rw <- sqrt(wmat[, j])
wmlmList[[j]] <- .lm.fit(rw * xmat, rw * ymat[, j])
}
return(wmlmList)
}
その使用の簡単な例を考えてみましょう:
## a toy dataset of 200 data with 3 numerical variables and 1 factor variable
dat <- data.frame(x1 = rnorm(200), x2 = rnorm(200), x3 = rnorm(200), f = gl(5, 40, labels = letters[1:5]))
## consider a model `~ x1 + poly(x3, 3) + x2 * f`
## we construct the non-weighted model matrix
xmat <- model.matrix (~ x1 + poly(x3, 3) + x2 * f, dat)
## now let's assume we have 100 model responses as well as 100 sets of weights
ymat <- matrix(rnorm(200 * 100), 200)
wmat <- matrix(runif(200 * 100), 200)
## Let's call `wmlm`:
fit <- wmlm (xmat, ymat, wmat)
.lm.fit
さらにモデルを推論するための重要な情報を返し、完全lmObject
はこれらのエントリのほとんどを継承します。
## take the first fitted model as an example
str(fit[[1]])
#$ qr : num [1:200, 1:14] -10.4116 0.061 0.0828 0.0757 0.0698 ...
# ..- attr(*, "assign")= int [1:14] 0 1 2 2 2 3 4 4 4 4 ...
# ..- attr(*, "contrasts")=List of 1
# .. ..$ f: chr "contr.treatment"
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:200] "1" "2" "3" "4" ...
# .. ..$ : chr [1:14] "(Intercept)" "x1" "poly(x3, 3)1" "poly(x3, 3)2" ...
#$ coefficients: num [1:14] 0.1184 -0.0506 0.3032 0.1643 0.4269 ...
#$ residuals : num [1:200] -0.7311 -0.0795 -0.2495 0.4097 0.0495 ...
#$ effects : num [1:200] -0.351 -0.36 0.145 0.182 0.291 ...
#$ rank : int 14
#$ pivot : int [1:14] 1 2 3 4 5 6 7 8 9 10 ...
#$ qraux : num [1:14] 1.06 1.13 1.07 1.05 1.01 ...
#$ tol : num 1e-07
#$ pivoted : logi FALSE
の結果には、、、 、 など.lm.fit
の一般的な関数がサポートされていません。ただし、線形モデルの推論は簡単なので、自分で計算するのは簡単です (背後にある理論を知っている場合)。summary
anova
predict
plot
- 回帰係数の t 統計表 (から
$qr
);
- F 統計量と ANOVA テーブル (から
$effects
);
- 残差標準誤差、R-2 乗および調整済み R-2 乗 (
$residulas
およびから$rank
)。
最後に、ベンチマークを提供します。
library(microbenchmark)
microbenchmark(wmlm = {wmlm (xmat, ymat, wmat);},
lm = {for (i in 1:ncol(ymat))
lm(ymat[, i] ~ x1 + poly(x3, 3) + x2 * f, dat, weights = wmat[, i]);} )
#Unit: milliseconds
# expr min lq mean median uq max neval cld
# wmlm 20.84512 23.02756 27.29539 24.49314 25.9027 79.8894 100 a
# lm 400.53000 405.10622 430.09787 414.42152 442.2640 535.9144 100 b
したがって、17.25 倍のブーストが見られます (中央値に基づく)。