我想在不使用"lm"的情况下计算R中的普通最小二乘(OLS)估计,这有几个原因。首先,考虑到数据大小在我的情况下是一个问题,"lm"还计算了很多我不需要的东西(比如拟合值)。其次,我希望能够在用另一种语言(例如,用C和GSL)实现OLS之前,先用R自己实现OLS。
正如您所知,模型是:Y=Xb+E;其中E~N(0,sigma^2)。如下所述,b是具有两个参数的向量,即平均值(b0)和另一个系数(b1)。最后,对于我要做的每一个线性回归,我想要b1(效应大小)、其标准误差、sigma^2(残差方差)和R^2(决定系数)的估计。
这是我的数据。我有N个样本(例如个人,N~=100)。对于每个样本,我有Y个输出(响应变量,Y~=10^3)和X个点(解释变量,X~=10^6)。我想单独处理Y输出,即我想启动Y线性回归:一个用于输出1,一个用于输入2,等等。此外,我想使用解释变量一个:对于输出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~neneneba xi)的每一列的良好特性,而我必须使用"应用"y次(对于每个yi~nenenebb xi)。但这是我的代码变慢的唯一原因吗?你们中有人知道如何改进吗?
(很抱歉问了这么长的问题,但我确实试图提供一个最小但全面的例子。)
> sessionInfo()
R version 2.12.2 (2011-02-25)
Platform: x86_64-redhat-linux-gnu (64-bit)
查看CRAN上RcppArmadillo包中的fastLm()
函数。在此之前,RcppGSL中也有类似的fastLm()
,但您可能想要基于Armadillo的解决方案。我在以前的演示中有一些幻灯片(关于带R的HPC),展示了速度的提高。
此外,请注意帮助页面中的提示,即与退化模型矩阵相关的X'X的直逆相比,更好的"枢轴"方法。
根据Marek的评论,以下是比较RppArmadillo:包中内置函数"lm"one_answers"lm.fit"、我自己的函数"fastLm"one_answers"fastLmPure"的结果
> 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,我们将不得不自己计算。Idem,我们可能想在"lm"上使用"fastLm",但如果我们想要西格玛的估计,我们必须自己计算。用自定义的R函数计算这些事情可能不是很有效(与"lm"所做的相比)。
考虑到这一切,我暂时会继续使用"lm",但实际上德克对"fastLm"的评论是一个很好的建议(这就是我选择他的答案的原因,因为其他人应该对此感兴趣)。