運行lm具有多個響應和權重

Question

我必須將具有相同模型矩陣的線性模型擬合爲多個響應。通過將響應指定爲矩陣而不是矢量，可以很容易地在R中完成此操作。用這種方法計算速度非常快。

現在我還想給模型添加權重，以對應答案的準確性。因此，對於每個響應向量，我還需要不同的權重向量。但是，lm只允許將權重作爲向量而不是矩陣輸入。所以我無法批量輸入權重，需要分別運行lm。這樣計算會變得更慢。

有沒有辦法在批處理模式下運行這些類型的模型，而不會多次調用lm？

Answer 1

現在我還想給模型添加權重，以對應答案的準確性。因此，對於每個響應向量，我還需要不同的權重向量。但是，lm允許僅將權重作爲向量而不是矩陣輸入。所以我無法批量輸入權重，需要分別運行lm。這樣計算會變得更慢。

如在Fitting a linear model with multiple LHS中所解釋的，「mlm」的效率要求所有LHS響應的共享模型矩陣。然而，加權迴歸不提供模型矩陣的重用，因爲對於不同的權重組，響應y和模型矩陣X需要重新縮放。請閱讀R: lm() result differs when using weights argument and when using manually reweighted data瞭解加權迴歸的工作原理。

有沒有辦法在批處理模式下運行這些類型的模型，而不會重複調用lm？

這取決於你想要什麼。如果你需要一個完整的lmObject，那麼你必須每次打電話lm。如果你只想要係數，你可以使用.lm.fit。上面的第二個鏈接已經演示了使用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等。但線性模型的推斷是容易的，所以它是簡單的計算你自己（只要你知道背後的理論）：

1. 迴歸係數的t統計量表（從$qr）;
2. F統計量和方差分析表（來自$effects）;
3. 殘餘標準誤差，R平方和調整的R平方（從$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倍提高被認爲是（基於平均時間）。