Hello, A collegue of mine has compared the runtime of a linear model + anova in SAS and S+. He got the same results, but SAS took a bit more than a minute whereas S+ took 17 minutes. I've tried it in R (1.9.0) and it took 15 min. Neither machine run out of memory, and I assume that all machines have similar hardware, but the S+ and SAS machines are on windows whereas the R machine is Redhat Linux 7.2. My question is if I'm doing something wrong (technically) calling the lm routine, or (if not), how I can optimize the call to lm or even using an alternative to lm. I'd like to run about 12,000 of these models in R (for a gene expression experiment - one model per gene, which would take far too long). I've run the follwong code in R (and S+):> options(contrasts=c('contr.helmert', 'contr.poly'))The 1st colum is the value to be modeled, and the others are factors.> names(df.gene1data) <- c("Va", "Ba", "Ti", "Do", "Ar", "Pr") > df[c(1:2,1343:1344),]Va Do Ti Ba Ar Pr 1 2.317804 000mM 24h NEW 1 1 2 2.495390 000mM 24h NEW 2 1 8315 2.979641 025mM 04h PRG 83 16 8415 4.505787 000mM 04h PRG 84 16 this is a dataframe with 1344 rows. x <- Sys.time(); wlm <- lm(Va ~ Ba+Ti+Do+Pr+Ba:Ti+Ba:Do+Ba:Pr+Ti:Do+Ti:Pr+Do:Pr+Ba:Ti:Do+Ba:Ti:Pr+Ba:Do:Pr+Ti:Do:Pr+Ba:Ti:Do:Pr+(Ba:Ti:Do)/Ar, data=df, singular=T); difftime(Sys.time(), x) Time difference of 15.33333 mins> anova(wlm)Analysis of Variance Table Response: Va Df Sum Sq Mean Sq F value Pr(>F) Ba 2 0.1 0.1 0.4262 0.653133 Ti 1 2.6 2.6 16.5055 5.306e-05 *** Do 4 6.8 1.7 10.5468 2.431e-08 *** Pr 15 5007.4 333.8 2081.8439 < 2.2e-16 *** Ba:Ti 2 3.2 1.6 9.8510 5.904e-05 *** Ba:Do 7 2.8 0.4 2.5054 0.014943 * Ba:Pr 30 80.6 2.7 16.7585 < 2.2e-16 *** Ti:Do 4 8.7 2.2 13.5982 9.537e-11 *** Ti:Pr 15 2.4 0.2 1.0017 0.450876 Do:Pr 60 10.2 0.2 1.0594 0.358551 Ba:Ti:Do 7 1.4 0.2 1.2064 0.296415 Ba:Ti:Pr 30 5.6 0.2 1.1563 0.259184 Ba:Do:Pr 105 14.2 0.1 0.8445 0.862262 Ti:Do:Pr 60 14.8 0.2 1.5367 0.006713 ** Ba:Ti:Do:Pr 105 15.8 0.2 0.9382 0.653134 Ba:Ti:Do:Ar 56 26.4 0.5 2.9434 2.904e-11 *** Residuals 840 134.7 0.2 The corresponding SAS program from my collegue is: proc glm data = "the name of the data set"; class B T D A P; model V = B T D P B*T B*D B*P T*D T*P D*P B*T*D B*T*P B*D*P T*D*P B*T*D*P A(B*T*D); run; Note, V = Va, B = Ba, T = Ti, D = Do, P = Pr, A = Ar of the R-example kind regards + thanks a lot for your help, Arne
The way to time things in R is system.time(). Without knowing much more about your problem we can only guess where R is spending the time. But you can find out by profiling -- see `Writing R Extensions'. If you want multiple fits with the same design matrix (do you?) you could look at the code of lm and call lm.fit repeatedly yourself. On Mon, 10 May 2004 Arne.Muller at aventis.com wrote:> Hello, > > A collegue of mine has compared the runtime of a linear model + anova in SAS and S+. He got the same results, but SAS took a bit more than a minute whereas S+ took 17 minutes. I've tried it in R (1.9.0) and it took 15 min. Neither machine run out of memory, and I assume that all machines have similar hardware, but the S+ and SAS machines are on windows whereas the R machine is Redhat Linux 7.2. > > My question is if I'm doing something wrong (technically) calling the lm routine, or (if not), how I can optimize the call to lm or even using an alternative to lm. I'd like to run about 12,000 of these models in R (for a gene expression experiment - one model per gene, which would take far too long). > > I've run the follwong code in R (and S+): > > > options(contrasts=c('contr.helmert', 'contr.poly')) > > The 1st colum is the value to be modeled, and the others are factors. > > > names(df.gene1data) <- c("Va", "Ba", "Ti", "Do", "Ar", "Pr") > > df[c(1:2,1343:1344),] > Va Do Ti Ba Ar Pr > 1 2.317804 000mM 24h NEW 1 1 > 2 2.495390 000mM 24h NEW 2 1 > 8315 2.979641 025mM 04h PRG 83 16 > 8415 4.505787 000mM 04h PRG 84 16 > > this is a dataframe with 1344 rows. > > x <- Sys.time(); > wlm <- lm(Va ~ > Ba+Ti+Do+Pr+Ba:Ti+Ba:Do+Ba:Pr+Ti:Do+Ti:Pr+Do:Pr+Ba:Ti:Do+Ba:Ti:Pr+Ba:Do:Pr+Ti:Do:Pr+Ba:Ti:Do:Pr+(Ba:Ti:Do)/Ar, data=df, singular=T); > difftime(Sys.time(), x) > > Time difference of 15.33333 mins > > > anova(wlm) > Analysis of Variance Table > > Response: Va > Df Sum Sq Mean Sq F value Pr(>F) > Ba 2 0.1 0.1 0.4262 0.653133 > Ti 1 2.6 2.6 16.5055 5.306e-05 *** > Do 4 6.8 1.7 10.5468 2.431e-08 *** > Pr 15 5007.4 333.8 2081.8439 < 2.2e-16 *** > Ba:Ti 2 3.2 1.6 9.8510 5.904e-05 *** > Ba:Do 7 2.8 0.4 2.5054 0.014943 * > Ba:Pr 30 80.6 2.7 16.7585 < 2.2e-16 *** > Ti:Do 4 8.7 2.2 13.5982 9.537e-11 *** > Ti:Pr 15 2.4 0.2 1.0017 0.450876 > Do:Pr 60 10.2 0.2 1.0594 0.358551 > Ba:Ti:Do 7 1.4 0.2 1.2064 0.296415 > Ba:Ti:Pr 30 5.6 0.2 1.1563 0.259184 > Ba:Do:Pr 105 14.2 0.1 0.8445 0.862262 > Ti:Do:Pr 60 14.8 0.2 1.5367 0.006713 ** > Ba:Ti:Do:Pr 105 15.8 0.2 0.9382 0.653134 > Ba:Ti:Do:Ar 56 26.4 0.5 2.9434 2.904e-11 *** > Residuals 840 134.7 0.2 > > The corresponding SAS program from my collegue is: > > proc glm data = "the name of the data set"; > > class B T D A P; > > model V = B T D P B*T B*D B*P T*D T*P D*P B*T*D B*T*P B*D*P T*D*P B*T*D*P A(B*T*D); > > run; > > Note, V = Va, B = Ba, T = Ti, D = Do, P = Pr, A = Ar of the R-example-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
Hello, thanks for your reply. I've now done the profiling, and I interpret that the most time is spend in the fortran routine(s): Each sample represents 0.02 seconds. Total run time: 920.219999999453 seconds. Total seconds: time spent in function and callees. Self seconds: time spent in function alone. % total % self total seconds self seconds name 100.00 920.22 0.02 0.16 "lm" 99.96 919.88 0.10 0.88 "lm.fit" 99.74 917.84 99.74 917.84 ".Fortran" 0.07 0.66 0.02 0.14 "storage.mode<-" 0.06 0.52 0.00 0.00 "eval" 0.06 0.52 0.04 0.34 "as.double" 0.02 0.22 0.02 0.22 "colnames<-" 0.02 0.20 0.02 0.20 "structure" 0.02 0.18 0.02 0.18 "model.matrix.default" 0.02 0.18 0.02 0.18 "as.double.default" 0.02 0.18 0.00 0.00 "model.matrix" 0.01 0.08 0.01 0.08 "list" % self % total self seconds total seconds name 99.74 917.84 99.74 917.84 ".Fortran" 0.10 0.88 99.96 919.88 "lm.fit" 0.04 0.34 0.06 0.52 "as.double" 0.02 0.22 0.02 0.22 "colnames<-" 0.02 0.20 0.02 0.20 "structure" 0.02 0.18 0.02 0.18 "as.double.default" 0.02 0.18 0.02 0.18 "model.matrix.default" 0.02 0.16 100.00 920.22 "lm" 0.02 0.14 0.07 0.66 "storage.mode<-" 0.01 0.08 0.01 0.08 "list" I guess this actually means I cannot do anything about it ... other than maybe splitting the problem into different (independaent parts - which I actually may be able to). Regarding the usage of lm.fit instead of lm, this might be a good idea, since I am using the same model.matrix for all fits! However, I'd need to recreate an lm object from the output, because I'd like to run the anova function on this. I'll first do some profiling on lm versus lm.fit for the 12,000 models ... kind regards + thanks again for your help, Arne -- Arne Muller, Ph.D. Toxicogenomics, Aventis Pharma arne dot muller domain=aventis com> -----Original Message----- > From: Prof Brian Ripley [mailto:ripley at stats.ox.ac.uk] > Sent: 11 May 2004 09:08 > To: Muller, Arne PH/FR > Cc: r-help at stat.math.ethz.ch > Subject: Re: [R] R versus SAS: lm performance > > > The way to time things in R is system.time(). > > Without knowing much more about your problem we can only > guess where R is > spending the time. But you can find out by profiling -- see > `Writing R > Extensions'. > > If you want multiple fits with the same design matrix (do you?) you > could look at the code of lm and call lm.fit repeatedly yourself. > > On Mon, 10 May 2004 Arne.Muller at aventis.com wrote: > > > Hello, > > > > A collegue of mine has compared the runtime of a linear > model + anova in SAS and S+. He got the same results, but SAS > took a bit more than a minute whereas S+ took 17 minutes. > I've tried it in R (1.9.0) and it took 15 min. Neither > machine run out of memory, and I assume that all machines > have similar hardware, but the S+ and SAS machines are on > windows whereas the R machine is Redhat Linux 7.2. > > > > My question is if I'm doing something wrong (technically) > calling the lm routine, or (if not), how I can optimize the > call to lm or even using an alternative to lm. I'd like to > run about 12,000 of these models in R (for a gene expression > experiment - one model per gene, which would take far too long). > > > > I've run the follwong code in R (and S+): > > > > > options(contrasts=c('contr.helmert', 'contr.poly')) > > > > The 1st colum is the value to be modeled, and the others > are factors. > > > > > names(df.gene1data) <- c("Va", "Ba", "Ti", "Do", "Ar", "Pr") > > > df[c(1:2,1343:1344),] > > Va Do Ti Ba Ar Pr > > 1 2.317804 000mM 24h NEW 1 1 > > 2 2.495390 000mM 24h NEW 2 1 > > 8315 2.979641 025mM 04h PRG 83 16 > > 8415 4.505787 000mM 04h PRG 84 16 > > > > this is a dataframe with 1344 rows. > > > > x <- Sys.time(); > > wlm <- lm(Va ~ > > > Ba+Ti+Do+Pr+Ba:Ti+Ba:Do+Ba:Pr+Ti:Do+Ti:Pr+Do:Pr+Ba:Ti:Do+Ba:Ti > :Pr+Ba:Do:Pr+Ti:Do:Pr+Ba:Ti:Do:Pr+(Ba:Ti:Do)/Ar, data=df, singular=T); > > difftime(Sys.time(), x) > > > > Time difference of 15.33333 mins > > > > > anova(wlm) > > Analysis of Variance Table > > > > Response: Va > > Df Sum Sq Mean Sq F value Pr(>F) > > Ba 2 0.1 0.1 0.4262 0.653133 > > Ti 1 2.6 2.6 16.5055 5.306e-05 *** > > Do 4 6.8 1.7 10.5468 2.431e-08 *** > > Pr 15 5007.4 333.8 2081.8439 < 2.2e-16 *** > > Ba:Ti 2 3.2 1.6 9.8510 5.904e-05 *** > > Ba:Do 7 2.8 0.4 2.5054 0.014943 * > > Ba:Pr 30 80.6 2.7 16.7585 < 2.2e-16 *** > > Ti:Do 4 8.7 2.2 13.5982 9.537e-11 *** > > Ti:Pr 15 2.4 0.2 1.0017 0.450876 > > Do:Pr 60 10.2 0.2 1.0594 0.358551 > > Ba:Ti:Do 7 1.4 0.2 1.2064 0.296415 > > Ba:Ti:Pr 30 5.6 0.2 1.1563 0.259184 > > Ba:Do:Pr 105 14.2 0.1 0.8445 0.862262 > > Ti:Do:Pr 60 14.8 0.2 1.5367 0.006713 ** > > Ba:Ti:Do:Pr 105 15.8 0.2 0.9382 0.653134 > > Ba:Ti:Do:Ar 56 26.4 0.5 2.9434 2.904e-11 *** > > Residuals 840 134.7 0.2 > > > > The corresponding SAS program from my collegue is: > > > > proc glm data = "the name of the data set"; > > > > class B T D A P; > > > > model V = B T D P B*T B*D B*P T*D T*P D*P B*T*D B*T*P B*D*P > T*D*P B*T*D*P A(B*T*D); > > > > run; > > > > Note, V = Va, B = Ba, T = Ti, D = Do, P = Pr, A = Ar of the > R-example > > -- > Brian D. Ripley, ripley at stats.ox.ac.uk > Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ > University of Oxford, Tel: +44 1865 272861 (self) > 1 South Parks Road, +44 1865 272866 (PA) > Oxford OX1 3TG, UK Fax: +44 1865 272595 > >
<Arne.Muller at aventis.com> writes:> Hello, > > A collegue of mine has compared the runtime of a linear model + anova in SAS and S+. He got the same results, but SAS took a bit more than a minute whereas S+ took 17 minutes. I've tried it in R (1.9.0) and it took 15 min. Neither machine run out of memory, and I assume that all machines have similar hardware, but the S+ and SAS machines are on windows whereas the R machine is Redhat Linux 7.2. > > My question is if I'm doing something wrong (technically) calling the lm routine, or (if not), how I can optimize the call to lm or even using an alternative to lm. I'd like to run about 12,000 of these models in R (for a gene expression experiment - one model per gene, which would take far too long). > > I've run the follwong code in R (and S+):... As Brian Ripley mentioned, you could save the model matrix and use it with each of your responses. Versions 0.8-1 and later of the Matrix package have a vignette that provides comparative timings of various ways of obtaining the least squares estimates. If you use the classes from the Matrix package and create and save the crossproduct of the model matrix mm = as(model.matrix(Va ~ Ba+Ti..., df), "geMatrix") cprod = crossprod(mm) then successive calls to coef = solve(cprod, crossprod(mm, df$Va)) will produce the coefficient estimates much faster than will calls to lm, which each do all the work of generating and decomposing the very large model matrix. Note that this method only produces the coefficient estimates, which may be enough for your purposes. Also, this method will not handle missing data or rank-deficient model matrices in the elegant way that lm does. If you are doing this 12,000 times it may be worthwhile checking if the sparse matrix formulation mmS = as(mm, "cscMatrix") cprodS = crossprod(mmS) is faster. The dense matrix formulation (but not the sparse) can benefit from installation of optimized BLAS routines such as Atlas or Goto's BLAS. -- Douglas Bates bates at stat.wisc.edu Statistics Department 608/262-2598 University of Wisconsin - Madison http://www.stat.wisc.edu/~bates/
I tried the following on an Opteron 248, R-1.9.0 w/Goto's BLAS:> y <- matrix(rnorm(14000*1344), 1344) > x <- matrix(runif(1344*503),1344) > system.time(fit <- lm(y~x))[1] 106.00 55.60 265.32 0.00 0.00 The resulting fit object is over 600MB. (The coefficient compoent is a 504 x 14000 matrix.) If I'm not mistaken, SAS sweeps on the extended cross product matrix to fit regression models. That, I believe, in usually faster than doing QR decomposition on the model matrix itself, but there are trade-offs. You could try what Prof. Bates suggested. Andy> From: Arne.Muller at aventis.com > > Hello, > > thanks for your reply. I've now done the profiling, and I > interpret that the most time is spend in the fortran routine(s): > > Each sample represents 0.02 seconds. > Total run time: 920.219999999453 seconds. > > Total seconds: time spent in function and callees. > Self seconds: time spent in function alone. > > % total % self > total seconds self seconds name > 100.00 920.22 0.02 0.16 "lm" > 99.96 919.88 0.10 0.88 "lm.fit" > 99.74 917.84 99.74 917.84 ".Fortran" > 0.07 0.66 0.02 0.14 "storage.mode<-" > 0.06 0.52 0.00 0.00 "eval" > 0.06 0.52 0.04 0.34 "as.double" > 0.02 0.22 0.02 0.22 "colnames<-" > 0.02 0.20 0.02 0.20 "structure" > 0.02 0.18 0.02 0.18 "model.matrix.default" > 0.02 0.18 0.02 0.18 "as.double.default" > 0.02 0.18 0.00 0.00 "model.matrix" > 0.01 0.08 0.01 0.08 "list" > > % self % total > self seconds total seconds name > 99.74 917.84 99.74 917.84 ".Fortran" > 0.10 0.88 99.96 919.88 "lm.fit" > 0.04 0.34 0.06 0.52 "as.double" > 0.02 0.22 0.02 0.22 "colnames<-" > 0.02 0.20 0.02 0.20 "structure" > 0.02 0.18 0.02 0.18 "as.double.default" > 0.02 0.18 0.02 0.18 "model.matrix.default" > 0.02 0.16 100.00 920.22 "lm" > 0.02 0.14 0.07 0.66 "storage.mode<-" > 0.01 0.08 0.01 0.08 "list" > > I guess this actually means I cannot do anything about it ... > other than maybe splitting the problem into different > (independaent parts - which I actually may be able to). > > Regarding the usage of lm.fit instead of lm, this might be a > good idea, since I am using the same model.matrix for all > fits! However, I'd need to recreate an lm object from the > output, because I'd like to run the anova function on this. > I'll first do some profiling on lm versus lm.fit for the > 12,000 models ... > > kind regards + thanks again for your help, > > Arne > > -- > Arne Muller, Ph.D. > Toxicogenomics, Aventis Pharma > arne dot muller domain=aventis com > > > -----Original Message----- > > From: Prof Brian Ripley [mailto:ripley at stats.ox.ac.uk] > > Sent: 11 May 2004 09:08 > > To: Muller, Arne PH/FR > > Cc: r-help at stat.math.ethz.ch > > Subject: Re: [R] R versus SAS: lm performance > > > > > > The way to time things in R is system.time(). > > > > Without knowing much more about your problem we can only > > guess where R is > > spending the time. But you can find out by profiling -- see > > `Writing R > > Extensions'. > > > > If you want multiple fits with the same design matrix (do you?) you > > could look at the code of lm and call lm.fit repeatedly yourself. > > > > On Mon, 10 May 2004 Arne.Muller at aventis.com wrote: > > > > > Hello, > > > > > > A collegue of mine has compared the runtime of a linear > > model + anova in SAS and S+. He got the same results, but SAS > > took a bit more than a minute whereas S+ took 17 minutes. > > I've tried it in R (1.9.0) and it took 15 min. Neither > > machine run out of memory, and I assume that all machines > > have similar hardware, but the S+ and SAS machines are on > > windows whereas the R machine is Redhat Linux 7.2. > > > > > > My question is if I'm doing something wrong (technically) > > calling the lm routine, or (if not), how I can optimize the > > call to lm or even using an alternative to lm. I'd like to > > run about 12,000 of these models in R (for a gene expression > > experiment - one model per gene, which would take far too long). > > > > > > I've run the follwong code in R (and S+): > > > > > > > options(contrasts=c('contr.helmert', 'contr.poly')) > > > > > > The 1st colum is the value to be modeled, and the others > > are factors. > > > > > > > names(df.gene1data) <- c("Va", "Ba", "Ti", "Do", "Ar", "Pr") > > > > df[c(1:2,1343:1344),] > > > Va Do Ti Ba Ar Pr > > > 1 2.317804 000mM 24h NEW 1 1 > > > 2 2.495390 000mM 24h NEW 2 1 > > > 8315 2.979641 025mM 04h PRG 83 16 > > > 8415 4.505787 000mM 04h PRG 84 16 > > > > > > this is a dataframe with 1344 rows. > > > > > > x <- Sys.time(); > > > wlm <- lm(Va ~ > > > > > Ba+Ti+Do+Pr+Ba:Ti+Ba:Do+Ba:Pr+Ti:Do+Ti:Pr+Do:Pr+Ba:Ti:Do+Ba:Ti > > :Pr+Ba:Do:Pr+Ti:Do:Pr+Ba:Ti:Do:Pr+(Ba:Ti:Do)/Ar, data=df, > singular=T); > > > difftime(Sys.time(), x) > > > > > > Time difference of 15.33333 mins > > > > > > > anova(wlm) > > > Analysis of Variance Table > > > > > > Response: Va > > > Df Sum Sq Mean Sq F value Pr(>F) > > > Ba 2 0.1 0.1 0.4262 0.653133 > > > Ti 1 2.6 2.6 16.5055 5.306e-05 *** > > > Do 4 6.8 1.7 10.5468 2.431e-08 *** > > > Pr 15 5007.4 333.8 2081.8439 < 2.2e-16 *** > > > Ba:Ti 2 3.2 1.6 9.8510 5.904e-05 *** > > > Ba:Do 7 2.8 0.4 2.5054 0.014943 * > > > Ba:Pr 30 80.6 2.7 16.7585 < 2.2e-16 *** > > > Ti:Do 4 8.7 2.2 13.5982 9.537e-11 *** > > > Ti:Pr 15 2.4 0.2 1.0017 0.450876 > > > Do:Pr 60 10.2 0.2 1.0594 0.358551 > > > Ba:Ti:Do 7 1.4 0.2 1.2064 0.296415 > > > Ba:Ti:Pr 30 5.6 0.2 1.1563 0.259184 > > > Ba:Do:Pr 105 14.2 0.1 0.8445 0.862262 > > > Ti:Do:Pr 60 14.8 0.2 1.5367 0.006713 ** > > > Ba:Ti:Do:Pr 105 15.8 0.2 0.9382 0.653134 > > > Ba:Ti:Do:Ar 56 26.4 0.5 2.9434 2.904e-11 *** > > > Residuals 840 134.7 0.2 > > > > > > The corresponding SAS program from my collegue is: > > > > > > proc glm data = "the name of the data set"; > > > > > > class B T D A P; > > > > > > model V = B T D P B*T B*D B*P T*D T*P D*P B*T*D B*T*P B*D*P > > T*D*P B*T*D*P A(B*T*D); > > > > > > run; > > > > > > Note, V = Va, B = Ba, T = Ti, D = Do, P = Pr, A = Ar of the > > R-example > > > > -- > > Brian D. Ripley, ripley at stats.ox.ac.uk > > Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ > > University of Oxford, Tel: +44 1865 272861 (self) > > 1 South Parks Road, +44 1865 272866 (PA) > > Oxford OX1 3TG, UK Fax: +44 1865 272595 > > > > > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > >
Thanks All, for your help. There seems to be a lot I can try to speed up the fits. However, I'd like to go for a much simpler model which I think is justified by the experiment itself, e.g; I may think about removing the nestinh "(Ba:Ti:Do)/Ar". The model matrix has 1344 rows and 2970 columns, and the rank of the matrix is 504. Therefore I think I should reformulate the model. I was just stroke my the massive difference in performance when my collegue told me about the difference between SAS and S+. kind regards, Arne -- Arne Muller, Ph.D. Toxicogenomics, Aventis Pharma arne dot muller domain=aventis com> -----Original Message----- > From: Liaw, Andy [mailto:andy_liaw at merck.com] > Sent: 11 May 2004 14:20 > To: Muller, Arne PH/FR; ripley at stats.ox.ac.uk > Cc: r-help at stat.math.ethz.ch > Subject: RE: [R] R versus SAS: lm performance > > > I tried the following on an Opteron 248, R-1.9.0 w/Goto's BLAS: > > > y <- matrix(rnorm(14000*1344), 1344) > > x <- matrix(runif(1344*503),1344) > > system.time(fit <- lm(y~x)) > [1] 106.00 55.60 265.32 0.00 0.00 > > The resulting fit object is over 600MB. (The coefficient > compoent is a 504 > x 14000 matrix.) > > If I'm not mistaken, SAS sweeps on the extended cross product > matrix to fit > regression models. That, I believe, in usually faster than doing QR > decomposition on the model matrix itself, but there are > trade-offs. You > could try what Prof. Bates suggested. > > Andy > > > From: Arne.Muller at aventis.com > > > > Hello, > > > > thanks for your reply. I've now done the profiling, and I > > interpret that the most time is spend in the fortran routine(s): > > > > Each sample represents 0.02 seconds. > > Total run time: 920.219999999453 seconds. > > > > Total seconds: time spent in function and callees. > > Self seconds: time spent in function alone. > > > > % total % self > > total seconds self seconds name > > 100.00 920.22 0.02 0.16 "lm" > > 99.96 919.88 0.10 0.88 "lm.fit" > > 99.74 917.84 99.74 917.84 ".Fortran" > > 0.07 0.66 0.02 0.14 "storage.mode<-" > > 0.06 0.52 0.00 0.00 "eval" > > 0.06 0.52 0.04 0.34 "as.double" > > 0.02 0.22 0.02 0.22 "colnames<-" > > 0.02 0.20 0.02 0.20 "structure" > > 0.02 0.18 0.02 0.18 "model.matrix.default" > > 0.02 0.18 0.02 0.18 "as.double.default" > > 0.02 0.18 0.00 0.00 "model.matrix" > > 0.01 0.08 0.01 0.08 "list" > > > > % self % total > > self seconds total seconds name > > 99.74 917.84 99.74 917.84 ".Fortran" > > 0.10 0.88 99.96 919.88 "lm.fit" > > 0.04 0.34 0.06 0.52 "as.double" > > 0.02 0.22 0.02 0.22 "colnames<-" > > 0.02 0.20 0.02 0.20 "structure" > > 0.02 0.18 0.02 0.18 "as.double.default" > > 0.02 0.18 0.02 0.18 "model.matrix.default" > > 0.02 0.16 100.00 920.22 "lm" > > 0.02 0.14 0.07 0.66 "storage.mode<-" > > 0.01 0.08 0.01 0.08 "list" > > > > I guess this actually means I cannot do anything about it ... > > other than maybe splitting the problem into different > > (independaent parts - which I actually may be able to). > > > > Regarding the usage of lm.fit instead of lm, this might be a > > good idea, since I am using the same model.matrix for all > > fits! However, I'd need to recreate an lm object from the > > output, because I'd like to run the anova function on this. > > I'll first do some profiling on lm versus lm.fit for the > > 12,000 models ... > > > > kind regards + thanks again for your help, > > > > Arne > > > > -- > > Arne Muller, Ph.D. > > Toxicogenomics, Aventis Pharma > > arne dot muller domain=aventis com > > > > > -----Original Message----- > > > From: Prof Brian Ripley [mailto:ripley at stats.ox.ac.uk] > > > Sent: 11 May 2004 09:08 > > > To: Muller, Arne PH/FR > > > Cc: r-help at stat.math.ethz.ch > > > Subject: Re: [R] R versus SAS: lm performance > > > > > > > > > The way to time things in R is system.time(). > > > > > > Without knowing much more about your problem we can only > > > guess where R is > > > spending the time. But you can find out by profiling -- see > > > `Writing R > > > Extensions'. > > > > > > If you want multiple fits with the same design matrix (do > you?) you > > > could look at the code of lm and call lm.fit repeatedly yourself. > > > > > > On Mon, 10 May 2004 Arne.Muller at aventis.com wrote: > > > > > > > Hello, > > > > > > > > A collegue of mine has compared the runtime of a linear > > > model + anova in SAS and S+. He got the same results, but SAS > > > took a bit more than a minute whereas S+ took 17 minutes. > > > I've tried it in R (1.9.0) and it took 15 min. Neither > > > machine run out of memory, and I assume that all machines > > > have similar hardware, but the S+ and SAS machines are on > > > windows whereas the R machine is Redhat Linux 7.2. > > > > > > > > My question is if I'm doing something wrong (technically) > > > calling the lm routine, or (if not), how I can optimize the > > > call to lm or even using an alternative to lm. I'd like to > > > run about 12,000 of these models in R (for a gene expression > > > experiment - one model per gene, which would take far too long). > > > > > > > > I've run the follwong code in R (and S+): > > > > > > > > > options(contrasts=c('contr.helmert', 'contr.poly')) > > > > > > > > The 1st colum is the value to be modeled, and the others > > > are factors. > > > > > > > > > names(df.gene1data) <- c("Va", "Ba", "Ti", "Do", "Ar", "Pr") > > > > > df[c(1:2,1343:1344),] > > > > Va Do Ti Ba Ar Pr > > > > 1 2.317804 000mM 24h NEW 1 1 > > > > 2 2.495390 000mM 24h NEW 2 1 > > > > 8315 2.979641 025mM 04h PRG 83 16 > > > > 8415 4.505787 000mM 04h PRG 84 16 > > > > > > > > this is a dataframe with 1344 rows. > > > > > > > > x <- Sys.time(); > > > > wlm <- lm(Va ~ > > > > > > > Ba+Ti+Do+Pr+Ba:Ti+Ba:Do+Ba:Pr+Ti:Do+Ti:Pr+Do:Pr+Ba:Ti:Do+Ba:Ti > > > :Pr+Ba:Do:Pr+Ti:Do:Pr+Ba:Ti:Do:Pr+(Ba:Ti:Do)/Ar, data=df, > > singular=T); > > > > difftime(Sys.time(), x) > > > > > > > > Time difference of 15.33333 mins > > > > > > > > > anova(wlm) > > > > Analysis of Variance Table > > > > > > > > Response: Va > > > > Df Sum Sq Mean Sq F value Pr(>F) > > > > Ba 2 0.1 0.1 0.4262 0.653133 > > > > Ti 1 2.6 2.6 16.5055 5.306e-05 *** > > > > Do 4 6.8 1.7 10.5468 2.431e-08 *** > > > > Pr 15 5007.4 333.8 2081.8439 < 2.2e-16 *** > > > > Ba:Ti 2 3.2 1.6 9.8510 5.904e-05 *** > > > > Ba:Do 7 2.8 0.4 2.5054 0.014943 * > > > > Ba:Pr 30 80.6 2.7 16.7585 < 2.2e-16 *** > > > > Ti:Do 4 8.7 2.2 13.5982 9.537e-11 *** > > > > Ti:Pr 15 2.4 0.2 1.0017 0.450876 > > > > Do:Pr 60 10.2 0.2 1.0594 0.358551 > > > > Ba:Ti:Do 7 1.4 0.2 1.2064 0.296415 > > > > Ba:Ti:Pr 30 5.6 0.2 1.1563 0.259184 > > > > Ba:Do:Pr 105 14.2 0.1 0.8445 0.862262 > > > > Ti:Do:Pr 60 14.8 0.2 1.5367 0.006713 ** > > > > Ba:Ti:Do:Pr 105 15.8 0.2 0.9382 0.653134 > > > > Ba:Ti:Do:Ar 56 26.4 0.5 2.9434 2.904e-11 *** > > > > Residuals 840 134.7 0.2 > > > > > > > > The corresponding SAS program from my collegue is: > > > > > > > > proc glm data = "the name of the data set"; > > > > > > > > class B T D A P; > > > > > > > > model V = B T D P B*T B*D B*P T*D T*P D*P B*T*D B*T*P B*D*P > > > T*D*P B*T*D*P A(B*T*D); > > > > > > > > run; > > > > > > > > Note, V = Va, B = Ba, T = Ti, D = Do, P = Pr, A = Ar of the > > > R-example > > > > > > -- > > > Brian D. Ripley, ripley at stats.ox.ac.uk > > > Professor of Applied Statistics, > http://www.stats.ox.ac.uk/~ripley/ > > > University of Oxford, Tel: +44 1865 272861 (self) > > > 1 South Parks Road, +44 1865 272866 (PA) > > > Oxford OX1 3TG, UK Fax: +44 1865 272595 > > > > > > > > > > ______________________________________________ > > R-help at stat.math.ethz.ch mailing list > > https://www.stat.math.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide! > > http://www.R-project.org/posting-guide.html > > > > > > > -------------------------------------------------------------- > ---------------- > Notice: This e-mail message, together with any attachments, > contains information of Merck & Co., Inc. (One Merck Drive, > Whitehouse Station, New Jersey, USA 08889), and/or its > affiliates (which may be known outside the United States as > Merck Frosst, Merck Sharp & Dohme or MSD and in Japan, as > Banyu) that may be confidential, proprietary copyrighted > and/or legally privileged. It is intended solely for the use > of the individual or entity named on this message. If you > are not the intended recipient, and have received this > message in error, please notify us immediately by reply > e-mail and then delete it from your system. > -------------------------------------------------------------- > ---------------- >