Tell me more ×
Cross Validated is a question and answer site for statisticians, data analysts, data miners and data visualization experts. It's 100% free, no registration required.
up vote 6 down vote favorite
1

Last night I started a complex calculation with gamm() and it took me... 

     user        system       elapsed 
    9259.76      326.05     9622.64 (s)

...meaning it took me 160 minutes or 2.67 hours for that calculation. The problem is that I have to do around 50 or even 100 more of these! So I was wondering if there is any way that could speed up these calculations. I compared the 32bit with the 64bit version (4gb) and R 2.12.2 to calculate a less complex gamm().

32bit solution

 User      System        elapsed 
 41.87        0.01       42.01

64 bit solution

  User      System      elapsed
 40.06        2.82       43.05

but it took even longer using 64bit!

My question now:

Would it help to simply buy more ram, for example 8GB DDR3? or would that be a waste of money? Or would the compiler package in R 2.13.0 be able to handle that properly? I do not think that rcpp can handle gamm() functions, or am I wrong?

any comments welcome!

the gamm() model call for the 160min process was: 

  g1 <- gamm(CountPP10M    ~ s(tempsurf,bs="cr") 
                           + s(salsurf,bs="cr") 
                           + s(speedsurf,bs="cr") 
                           + s(Usurf,bs="cr")
                           + s(Vsurf,bs="cr")   
                           + s(Wsurf,bs="cr")
                           + s(water_depth,bs="cr")
                           + s(distance.to.bridge,bs="cr")
                           + s(dist_land2,bs="cr")
                           + s(Dist_sventa,bs="cr"),
                           data=data,
                           random=list(ID_Station=~1),
                           family=poisson,
                           methods="REML",
                           control=lmc)
            )
share|improve this question
I just tried the compiler package in R v2.13.0 in 64bit modus. Now it took "only" 38.39 seconds, so roughly minus 5 seconds.. – Jens Jul 7 '11 at 7:45
...that's 11 percent – Jens Jul 7 '11 at 7:49
2  
The multicore won't speed up the single model as it is doing computations that are not parallel. My point was, if you have 4 cores, you can reduce your compute time by about 4 over the entire set of 50-100 models by doing them 4 at a time at the weekend, not one at a time as R only uses a single core. I fitted several thousand complex gamm() timeseries models with pretty large data sets and left my workstation going for a week on 3 cores to do the fitting. I'd worked out the code needed first on a small data set, then set it to work on the main problem. – Gavin Simpson Jul 7 '11 at 8:21
1  
I doubt the amount of RAM will help the single model timings. Check if you are maxing out on RAM during a model fit. If R still has plenty of head room then RAM is not the issue for the individual model. Where it might be an issue is if you follow the advice to break the job down into 4 chunks to run on separate cores (12-25 models in each). Then I could see you quickly running out of RAM - the models I mentioned fitting on 3 cores were on a workstation with 16GB of RAM and the computations routinely used ~10GB of that in total over the three processes. – Gavin Simpson Jul 7 '11 at 9:58
1  
rcpp is a package that allows easy wrapping/integration of C++ code within R and writing of C++ that exploits R object structures etc. gamm() is already coded in C. Yes it could be made faster (hence the gamm4 package which uses lme4 which is much faster at the mixed model computations than nlme was/is), but you can't make it faster using rcpp without rewriting the entire model in your own C++ using the C++ classes/API that rcpp provides. The help for gamm does warn you that you are stressing the computational ability of the code and your machine when fitting even modestly sized models. – Gavin Simpson Jul 7 '11 at 10:10
show 4 more comments

2 Answers

up vote 6 down vote accepted

You are not going to be able to achieve substantial speed-up here as most of the computation will be being done inside compiled C code.

If you are fitting correlation structures in gamm() then you can either simplify the correlation structure you want to fit (i.e. don't use corARMA(p=1, .....) when corAR1(....) would suffice. Or nest the correlations within years if you have many observations per year, rather than for the whole time interval.

If you aren't fitting correlation structures, gam() can fit simple random effects, and if you need more complex random effects, consider the gamm4 which is by the same author as mgcv but which uses the lme4 package (lmer()) instead of the slower/older nlme package (lme()).

You could try simpler bases for the smooth terms; bs = "cr" rather than the default thin-plate spline bases.

If all else fails, and you are just facing big-data issues, the best you can do is exploit multiple cores (manually split a job into ncores chunks and run them in BATCH mode over night, or by one of the parallel processing packages on R) and run models as the weekend. If you do this, make sure you wrap your gamm() calls in try() so that the whole job doesn't stop because you have a convergence problem part way through the run.

share|improve this answer
I think it really needs temporal autocorrelation see the link for a png imageshare.web.id/images/jgn53ohyvvmkthqj09e4.png so probably only too complex calculations. I could shrink it down by linearizing some smoothers. but I am not sure if they need smoother terms after dropping another factor. – Jens Jul 7 '11 at 8:34
1  
Fit without the correlation structure. That will speed things up immeasurably. Then look at the model and especially acf(resid(mod$lme, type = "normalised")) and look to see if the residuals are autocorrelated or not. Also, look at the smooth functions, see which might be simple and dropped or made linear. For others, perhaps try fixing the smoothness for some smooths to reasonable values; setting k to a suitable value and setting fx = TRUE in the s() terms you want to fix. – Gavin Simpson Jul 7 '11 at 9:52
1  
If that doesn't work, then you can drop the correlation structure, fit using gam() and it's new bs = "re" splines. Fix each smoother at a reasonable level of smoothness (see above). Once you have the fitted model, then look at residuals and decide what time series model you need for them. Then you could manually produce the covariance matrix of the residuals using that time series process and replace the one in the fitted GAM model object with the new one. That way you get proper standard errors on terms etc. This is getting into difficult territory but doable. – Gavin Simpson Jul 7 '11 at 9:55
up vote 0 down vote

If gamm() is in R code rather than C it might be worth using the byte-code compiler that is new in R 2.13. There is a new core package called compiler and you can compile a function using the cmpfun() function.

More details can be found here: http://www.r-bloggers.com/the-new-r-compiler-package-in-r-2-13-0-some-first-experiments/

share|improve this answer
@Jens tried that for small improvement. Most of the heavy lifting in gamm() is eventually done in lme() in the nlme package. That will be exploiting compiled C code for most of the fitting. – Gavin Simpson Jul 7 '11 at 10:03
hi, I am currently trying the compiler package and it seem to work quite well! I did the same 160min gamm() now within 5 minutes ! but I must admit, I changed the Count response to a gaussian response (I know...the beta regression problem again).anyway...5 minutes. Then I did the original one gamm() including (!) temporal autocorrelation i.e. correlation=corAR1(form = ~ 1 | ID_Station) and with compile it took only 42 minutes. so 75% better than the other one. I am still computing the original one with compile. will tell you in about 5-42 minutes :) – Jens Jul 7 '11 at 10:40
2  
@Jens the problem with that comparison is that with Gaussian errors (not response!) the function calls go gamm() -> lme() -> compiled code and iterate between the GAM and LMM bits. With non-Gaussian errors, you are using a PQL approach to fitting GLMMs in a LMM framework. So the calls now go gamm() -> glmmPQL() -> lme() etc and each of those steps need to be iterated to convergence. So you can't compare timing of compiler compiled code versus non-compiler code between Guassian and non-Gaussian runs. – Gavin Simpson Jul 7 '11 at 12:43
Yes, I see, I did not notice that PQL is not involved in fitting the gaussian. many thanks for that comment. and again you are right. compile was not helpful for doing the original gamm() with PQL. 161minutes... – Jens Jul 7 '11 at 13:16

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.