1

同じモデル マトリックスを使用して線形モデルを複数の応答に適合させる必要があります。これは、応答をベクトルではなく行列として指定することで、R で簡単に実行できます。このように計算は非常に高速です。

次に、応答の精度に対応する重みをモデルに追加したいと思います。したがって、応答ベクトルごとに、異なる重みベクトルも必要になります。ただし、lm行列ではなくベクトルとしてのみ重みを入力できます。そのため、重みをバッチで入力できず、lm応答ごとに個別に実行する必要がありました。この方法では、計算がはるかに遅くなります。

lm繰り返し呼び出すことなく、これらのタイプのモデルをバッチ モードで実行する方法はありますか?

4

1 に答える 1

2

次に、応答の精度に対応する重みをモデルに追加したいと思います。したがって、応答ベクトルごとに、異なる重みベクトルも必要になります。ただし、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の一般的な関数がサポートされていません。ただし、線形モデルの推論は簡単なので、自分で計算するのは簡単です (背後にある理論を知っている場合)。summaryanovapredictplot

  1. 回帰係数の t 統計表 (から$qr);
  2. F 統計量と ANOVA テーブル (から$effects);
  3. 残差標準誤差、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 倍のブーストが見られます (中央値に基づく)。

于 2016-11-20T12:31:05.280 に答える