[R] finding a faster way to run lm on rows of predictor matrix

["Dénes TÓTH"]

Hi,

you can solve the task by simple matrix algebra.
Note that I find it really inconvenient to have a matrix with variables in
rows and cases in columns, so I transposed your predictors matrix.

# your data
regress.y = rnorm(150)
predictors = matrix(rnorm(6000*150), ncol=150, nrow=6000)
tpreds <- t(predictors)

# compute coefficients
coefs <- apply(tpreds,2,function(x)
solve(crossprod(x),crossprod(x,regress.y)))

# compute residuals
resids <- regress.y - sweep(tpreds,2,coefs,"*")

(Note that resids shall be transposed back if you really insist on your
original matrix format.)

HTH,
 Denes


> Hi, everyone.
> I need to run lm with the same response vector but with varying predictor
> vectors. (i.e. 1 response vector on each individual 6,000 predictor
> vectors)
> After looking through the R archive, I found roughly 3 methods that has
> been suggested.
> Unfortunately, I need to run this task multiple times(~ 5,000 times) and
> would like to find a faster way than the existing methods.
> All three methods I have bellow run on my machine timed with system.time
> 13~20 seconds.
>
> The opposite direction of 6,000 response vectors and 1 predictor vectors,
> that is supported with lm runs really fast ~0.5 seconds.
> They are pretty much performing the same number of lm fits, so I was
> wondering if there was a faster way, before I try to code this up in c++.
>
> thanks!!
>
> ## sample data ###
> regress.y = rnorm(150)
> predictors = matrix(rnorm(6000*150), ncol=150, nrow=6000)
>
> ## method 1 ##
> data.residuals = t(apply(predictors, 1, function(x)( lm(regress.y ~ -1 +
> as.vector(x))$residuals)))
>
> user  system elapsed
>  15.076   0.048  15.128
>
> ## method 2 ##
> data.residuals = matrix(rep(0, nrow(predictors) * ncol(predictors)),
> nrow=nrow(predictors), ncol=ncol(predictors) )
>
> for( i in 1:nrow(predictors)){
>     pred = as.vector(predictors[i,])
>     data.residuals[i, ] = lm(regress.y ~ -1 + pred )$residuals
> }
>
>  user  system elapsed
>  13.841   0.012  13.854
>
> ## method 3 ##
> library(nlme)
>
> all.data <- data.frame( y=rep(regress.y, nrow(predictors)),
> x=c(t(predictors)), g=gl(nrow(predictors), ncol(predictors)) )
> all.fits <- lmList( y ~ x | g, data=all.data)
> data.residuals = matrix( residuals(all.fits), nrow=nrow(predictors),
> ncol=ncol(predictors))
>
> user  system elapsed
>  36.407   0.484  36.892
>
>
> ## the opposite direction, supported by lm ###
> lm(t(predictors) ~ -1 + regress.y)$residuals
>
>  user  system elapsed
>  0.500   0.120   0.613
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>