使用多个响应和权重运行 lm答案

【问题标题】：Run lm with multiple responses and weights使用多个响应和权重运行 lm
【发布时间】：2013-02-24 03:57:36
【问题描述】：

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

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

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

【问题讨论】：

欢迎来到本站。因为这个问题只是关于 R，所以在 StackOverflow 上会更好。我已将其标记为在那里迁移。
你看过lm.fit吗？
我认为不可能。在lm.wfit 的代码中，您会看到z <- .Call(C_Cdqrls, x * wts, y * wts, tol)。这意味着x 和y 都使用相同的权重进行转换。这是不可能的，如果你想为你的多个响应设置不同的权重。

标签： r regression linear-regression lm mlm

【解决方案1】：

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

正如Fitting a linear model with multiple LHS 中所述，“传销”的效率需要所有 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 等。但是线性模型的推断很容易，因此可以直接计算自己（提供你知道背后的理论）：

回归系数的t统计表（来自$qr）；
F-statistic 和 ANOVA 表（来自$effects）；
残差标准误差、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 倍的提升（基于中位时间）。

【讨论】：