1

我必须将具有相同模型矩阵的线性模型拟合到多个响应中。这可以通过将响应指定为矩阵而不是向量在 R 中轻松完成。以这种方式计算非常快。

现在我还想为模型添加与响应准确性相对应的权重。因此，对于每个响应向量，我还需要不同的权重向量。但是，lm只允许将权重作为向量而不是矩阵输入。所以我不能批量输入权重，必须分别lm为每个响应运行。这样计算会变得慢得多。

有没有办法以批处理模式运行这些类型的模型，而无需lm重复调用？

4

1 回答 1

2

现在我还想为模型添加与响应准确性相对应的权重。因此，对于每个响应向量，我还需要不同的权重向量。但是，lm允许仅将权重作为向量而不是矩阵输入。所以我不能批量输入权重，必须分别lm为每个响应运行。这样计算会变得慢得多。

如拟合具有多个 LHS 的线性模型中所述，“mlm”的效率需要所有 LHS 响应的共享模型矩阵。然而，加权回归没有给出模型矩阵的重用，因为对于一组不同的权重，响应y和模型矩阵都X需要重新调整。阅读R: lm() 结果在使用weights参数和使用手动重新加权数据以查看加权回归如何工作时有所不同。

有没有办法以批处理模式运行这些类型的模型，而无需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通用函数。但是线性模型的推断很容易，因此可以直接计算自己（前提是您知道背后的理论）：summaryanovapredictplot

  1. 回归系数的 t 统计表（来自$qr）；
  2. F 统计量和 ANOVA 表（来自$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 倍的提升（基于中位时间）。

于 2016-11-20T12:31:05.280 回答