replaced mvtnorm::rmvnorm with mvnfast::rmvn

[singmann]

This pull request replaces mvtnorm::rmvnorm with mvnfast::rmvn and adds a new argument, ncores, to the corresponding functions which is passed through to mvnfast::rmvn. This allows generating random variables from the multivariate normal distribution with multiple cores (via openMP), which can speed up the simulation quite considerably for large models. I have also removed all import statements for the old package.

Let me know if this is something you might consider. I am happy to discuss this should things be unclear.

[singmann]

I just noticed that there is a copy and paste error in this pull request in line 166 of derivatives.R. I am happy to submit a new pull request which corrects this.

I probably should have figured this out by running the unit tests after rebuilding the documentation, but I was not sure how this work here. Does this work with devtools::document() and devtools::test()`? If not, can you tell me how to do that in case I find more things?

[gavinsimpson]

@singmann Thanks for this. Ideally we'd have discussed this in an Issue first, but as I haven't exactly outlined how I'd like to receive contributions, that's on me :-)

Could you outline what the perceived advantages of mvnfast are beyond just speed? I'm assuming that there is a memory consideration for big models/data sets spread over n cores? What's the accuracy like? I'll take a look at the package myself but it would be good to get input on this.

I don't want to make a habit of making changes to the package that will affect backwards compatibility. However, the package is currently early/young enough in the dev cycle that now is a good time to consider this and once we make a decision we can stick with it going forward.

The tests are all standard testthat stuff. I don't use devtools myself much, but I assume that it will be able to run the tests in whatever way it interacts with testthat tests. I just run R CMD check as that tests much more than simply the unit test — it'll check that docs match the code etc, as well as running all the tests.

If you are running the tests, note that there are a couple of quite long-running tests that relate to a paper that we have submitted using gratia. You can turn this off by setting an environment variable NOT_CRAN to true, or perhaps this is something that devtools will manage for you? Otherwise the tests can take tens of minutes to complete. These aren't run normally on Travis either as it times out whilst waiting for the model to fit and it seems sub-optimal to delay all CI tests just for these couple of models.

[singmann]

Sorry for the delay, but I have done a bit more testing.

First of all, my main interest is speed. To this end, I have now compared both functions, mvtnorm::rmvnorm and mvnfast::rmvn, with a rather large gam model (50,000 data points, random effect with ~1000 levels, adaptive smooth across 20 factor levels, total model size is 3 GB). The difference in the confint method is quite large, even when not using multiple cores.

mvtnorm::rmvnorm needs over 1000 seconds for a simulation with 1000 simulations:

> system.time(
+ pred1 <- confint(mr_2, "timen2", type = "simultaneous", nsim = 1000)
+ )
   user  system elapsed 
1234.70    2.03 1237.04

mvnfast::rmvn is considerably faster (by at least a factor of 2 to 3), even when only using one core:

> system.time(
+ pred2 <- confint(mr_2, "timen2", type = "simultaneous", nsim = 1000, shift = TRUE)
+ )
   user  system elapsed 
 462.11    0.77  463.23 
> system.time(
+ pred3 <- confint(mr_2, "timen2", type = "simultaneous", 
+                  ncores = 8, nsim = 1000, shift = TRUE)
+ )
   user  system elapsed 
 672.30    0.95  315.14 

Note, I have rerun each call at least once and got essentially the same results to the one shown here each time. Also, the difference in the effect of using only one core to multiple cores for mvnfast increases the larger nsim gets. I do not show the results here, but that seems to be a quite clear pattern.

I also checked that all three methods produce equivalent results:

> all.equal(as.data.frame(pred1), 
+           as.data.frame(pred3), tolerance = 0.025)
[1] TRUE
> all.equal(as.data.frame(pred1), 
+           as.data.frame(pred2), tolerance = 0.025)
[1] TRUE
> all.equal(as.data.frame(pred2), 
+           as.data.frame(pred3), tolerance = 0.025)
[1] TRUE
> nrow(pred1)
[1] 4000

I cannot really say much about memory efficiency, as I am haven't yet figured out to profile the memory easily. But from just looking at the graph of the OS, mvnfast also seems to require less memory, and the memory footprint does not seem to increase with cores. I guess this is a consequence of using OpenMP.

From my perspective this generally suggests superiority of mvnfast. However, I suspect that changing the method in the call to mvtnorm::rmvnorm (e.g., to "chol") should also come with some improvements. If it will be similar large, I did not test.

Unfortunately, mvnfast does not come with a paper yet, but the author, Matteo Fasiolo, seems trustworthy. He works on GAMs and actually has papers with Simon Wood. We could @-him here as he is on github as well if you want to hear his opinion.

However, I am fine with not adding it to the package. It is of course your call. I know that backwards compatibility is important and the changes where really small and helped for my use case already. So not much work is lost if you reject it. In the end it was also easier to directly submit the pull request compared to opening an issue first. I had already done the changes.

One more thing. about the package. Which command do you use to go from the roxygen documentation to get the R documentation? You use markdown in the roxygen2 so I am afraid that testthat::document() will do more harm than the right thing (I mainly needed this to rewrite the NAMESPACE file after the changes to which functions are imported).

[singmann]

@gavinsimpson Did you have time to consider this pull request? I still feel like it would be a nice addition that can increase speed with little downsides. I also removed the small bug and made sure it works with the current version.

[gavinsimpson]

@singmann Yep; I've thought about this on an off whilst I've been busy with work and I agree with your evaluation and am happy to make this change. I'll take a look at merging this within in the next few days.

Thanks for doing this. I'll do a review and see what needs doing to merge this. Could you add yourself as a contributor (ctb) in the DESCRIPTION - if you're not familiar with how Authors@R works I can also do this when I merge the PR.

[singmann]

Great. I have made a few more changes to this pull request to remove all remaining issues and make it easy to merge.

1. I have removed a few more copy-and-paste errors, all code I changed should now work in all functions.
2. I have rebuild the documentation and NAMESPACE (simply using devtools::document()). This allows to run R CMD check, which throws no warnings or errors. Even Travis CI runs through now (no idea what is happening with AppVeyor though).
3. As asked, I have added myself as contributor to the DESCRIPTION

Let me know if there is anything else I can do to make this pull request easy for you.

[gavinsimpson]

Thanks for this PR @singmann. This is now merged with a few tweaks.

[singmann]

Thanks for accepting this. Should be quite helpful for larger models.