4

我想在不使用 "lm"的情况下计算 R 中的普通最小二乘 ( OLS ) 估计值，这有几个原因。首先，考虑到数据大小在我的情况下是一个问题，“lm”还计算了很多我不需要的东西（例如拟合值）。其次，我希望能够在 R 中自己实现 OLS，然后再用另一种语言（例如，在 C 中使用 GSL）。

你可能知道，模型是：Y=Xb+E；与 E ~ N(0, sigma^2)。如下详述，b 是具有 2 个参数的向量，均值 (b0) 和另一个系数 (b1)。最后，对于我将要做的每个线性回归，我想要 b1（效应大小）的估计值、它的标准误差、sigma^2（残差）和 R^2（确定系数）的估计值。

这是我的数据。我有 N 个样本（例如个人，N~=100）。对于每个样本，我有 Y 个输出（响应变量，Y~=10^3）和 X 个点（解释变量，X~=10^6）。我想分别处理 Y 输出，即。我想启动 Y 线性回归：一个用于输出 1，一个用于输出 2，等等。此外，我想使用解释变量一个 y 一个：对于输出 1，我想在第 1 点回归，然后在第 2 点回归，然后......最后在X点上。（我希望它很清楚......！）

这是我的R 代码，用于检查“lm”的速度与通过矩阵代数计算 OLS 估计我自己的速度。

首先，我模拟虚拟数据：

nb.samples <-  10  # N
nb.points <- 1000  # X
x <- matrix(data=replicate(nb.samples,sample(x=0:2,size=nb.points, replace=T)),
            nrow=nb.points, ncol=nb.samples, byrow=F,
            dimnames=list(points=paste("p",seq(1,nb.points),sep=""),
              samples=paste("s",seq(1,nb.samples),sep="")))
nb.outputs <- 10  # Y
y <- matrix(data=replicate(nb.outputs,rnorm(nb.samples)),
            nrow=nb.samples, ncol=nb.outputs, byrow=T,
            dimnames=list(samples=paste("s",seq(1,nb.samples),sep=""),
              outputs=paste("out",seq(1,nb.outputs),sep="")))

这是我自己在下面使用的函数：

GetResFromCustomLinReg <- function(Y, xi){ # both Y and xi are N-dim vectors
  n <- length(Y)
  X <- cbind(rep(1,n), xi)  #
  p <- 1      # nb of explanatory variables, besides the mean
  r <- p + 1  # rank of X: nb of indepdt explanatory variables
  inv.XtX <- solve(t(X) %*% X)
  beta.hat <- inv.XtX %*% t(X) %*% Y
  Y.hat <- X %*% beta.hat
  E.hat <- Y - Y.hat
  E2.hat <- (t(E.hat) %*% E.hat)
  sigma2.hat <- (E2.hat / (n - r))[1,1]
  var.covar.beta.hat <- sigma2.hat * inv.XtX
  se.beta.hat <- t(t(sqrt(diag(var.covar.beta.hat))))
  Y.bar <- mean(Y)
  R2 <- 1 - (E2.hat) / (t(Y-Y.bar) %*% (Y-Y.bar))
  return(c(beta.hat[2], se.beta.hat[2], sigma2.hat, R2))
}

这是我使用内置“lm”的代码：

res.bi.all <- apply(x, 1, function(xi){lm(y ~ xi)})

这是我的自定义 OLS 代码：

res.cm.all <- apply(x, 1, function(xi){apply(y, 2, GetResFromCustomLinReg, xi)})

当我使用上面给出的值运行此示例时，我得到：

> system.time( res.bi.all <- apply(x, 1, function(xi){lm(y ~ xi)}) )
   user  system elapsed
  2.526   0.000   2.528
> system.time( res.cm.all <- apply(x, 1, function(xi){apply(y, 2, GetResFromCustomLinReg, xi)}) )
   user  system elapsed
  4.561   0.000   4.561

（当然，当增加 N、X 和 Y 时，情况会变得更糟。）

当然，“lm”具有“自动”分别拟合响应矩阵的每一列（y~xi）的好特性，而我必须使用“应用”Y 次（对于每个 yi~xi）。但这是我的代码变慢的唯一原因吗？你们当中有人知道如何改进吗？

（对不起，这么长的问题，但我真的试图提供一个最小但全面的例子。）

> sessionInfo()
R version 2.12.2 (2011-02-25)
Platform: x86_64-redhat-linux-gnu (64-bit)
4

2 回答 2

7

看看CRAN 上的RcppArmadillo包中的fastLm()函数。在此之前的RcppGSL中也有类似的内容——但您可能需要基于Armadillo的解决方案。我在较旧的演示文稿（使用 R 的 HPC 上）中有一些幻灯片，显示了速度的提升。fastLm()

还要注意帮助页面中关于比 X'X 的直接逆更好的“枢轴”方法的提示，这可能与退化模型矩阵有关。

于 2011-06-04T04:24:06.233 回答
4

根据 Marek 的评论，下面是比较内置函数“lm”和“lm.fit”，我自己的函数，“fastLm”和“fastLmPure”包 RcppArmadillo 的结果：

> system.time( res1 <- apply(x, 1, function(xi){lm(y ~ xi)}) )
   user  system elapsed
  2.859   0.005   2.865
> system.time( res2 <- apply(x, 1, function(xi){apply(y, 2, GetResFromCustomLinReg, xi)}) )
   user  system elapsed
  4.620   0.004   4.626
> system.time( res3 <- apply(x, 1, function(xi){lm.fit(x=cbind(1,xi), y=y)}) )
   user  system elapsed
  0.454   0.004   0.458
> system.time( res4 <- apply(x, 1, function(xi){apply(y, 2, fastLm, x=cbind(1,xi))}) )
   user  system elapsed
  2.279   0.005   2.283
> system.time( res5 <- apply(x, 1, function(xi){apply(y, 2, fastLmPure, cbind(1,xi))}) )
   user  system elapsed
  1.053   0.003   1.056

但是，在比较这些数字时要小心。差异不仅在于不同的实现，还在于有效计算的结果：

> names(res1$p1)
 [1] "coefficients"  "residuals"     "effects"       "rank"        
 [5] "fitted.values" "assign"        "qr"            "df.residual" 
 [9] "xlevels"       "call"          "terms"         "model"       
> # res2 (from my own custom function) returns the estimate of beta, its standard error, the estimate of sigma and the R^2
> names(res3$p1)
[1] "coefficients"  "residuals"     "effects"       "rank"        
[5] "fitted.values" "assign"        "qr"            "df.residual" 
> names(res4$p1$out1)
[1] "coefficients"  "stderr"        "df"            "fitted.values"
[5] "residuals"     "call"        
> names(res5$p1$out1)
[1] "coefficients" "stderr"       "df"

例如，我们可能更喜欢使用“lm.fit”而不是“lm”，但如果我们需要 R^2，我们将不得不自己计算它。同上，我们可能想要使用“fastLm”而不是“lm”，但是如果我们想要估计 sigma，我们将不得不自己计算它。并且使用自定义 R 函数计算此类事情可能不是很有效（与“lm”所做的相比）。

鉴于这一切，我暂时将继续使用“lm”，但确实 Dirk 关于“fastLm”的评论是一个很好的建议（这就是我选择他的答案的原因，因为它应该会引起其他人的兴趣）。

于 2011-06-27T16:24:37.387 回答