Our new R package: R2jags

I have got emails occasionally from JAGS users, asking about our new R package: R2jags. Basically, R2jags runs JAGS via R and makes postanalysis easier to be done in R. Taking advantage of the functions provided by JAGS, rjags and R2WinBUGS, R2jags allows users to run BUGS in the same way as they would do it in R2WinBUGS. Nonetheless, R2jags has some powerful features that facilitate the BUGS model fitting:

  1. If your model does not converge, update it. If you used to be a R2WinBUGS user, you must feel frustrated that your model does not converge. The best thing you can do is to use your current states of parameters as the starting values for the next MCMC. On the other hand, in R2jags, you just type in the R console:

    fit <- update(fit) 
    

    This will update the model.

  2. If you have to shutdown your machine but the model is still not converged.... In R2jags, you can go ahead and save the R workspace and shutdown your machine. When you are ready to run the model again, load the workspace in R and type:

    recompile(fit)
    fit.upd <- update(fit)
    

    This will recompile the model, which means you can update the model again!!

If you want to explore more features of the R2jags, just type ?jags in the R console. The example code contains all the functions in the R2jags.

Installing JAGS and rjags can be tricky in the Mac OS or in the Linux system. I have a blog entry here that shows how this can be done. If you are not a Window user, this post might help you.

Recent changes in arm (arm > 1.1-8)

There are some changes in arm. Most of them are related to the big changes lme4 has made. If your arm version number is > arm_1.1.8, then you might want to read the followings (or if you are a reader of Data Analysis Using Regression and Multilevel/Hierarchical Models, you should read this):

New functions:

1. multicomp.plot(): mcplot() is the wraper for it. It plots the significant differnece of a simulated array. For example, see Figure 5 in this paper.

2. traceplot(): we made traceplot into a generic function. So now it identifies “mcmc.list” and “bugs” objects. The old traceplot() in coda will be masked. But the functionality of it will be passed into arm as the S4 method.

Changes related to lme4.

The new releass of lme4 (version # > 0.999375-16) has many changes. See the detailed here. Basically, it adopts new alogrithm, which makes lme4 quicker and more robust. Here are some changes that affects arm since some functions in arm depend on lme4:

1. lmer(): now lmer() breaks into lmer(), glmer(), and nlmer() (linear, generalized linear, nonlinear and generalized nonlinear mixed-effects models). So if you are fitting a generalized linear model, instead of using lmer() and specify the link function, you should use glmer() and specify the link function.

2. mcsamp(): mcsamp() depends on mcmcsamp() in lme4. It is a way to simulate the mixed model fitted by lme4. Whereas arm does provide a sim() function that does the same thing, mcsamp() can deal with the sigularity problem of the cov-variance matrix. The most common error message related to the sigularity problem when you use sim() to simulate a “mer” model is “Sigma is not postive definite!” However, the new mcmcsamp() contains some problems. It is still under development. So the new arm takes out mcsamp(). We will put it back when mcmcsamp() works properly. This means if you used to use mcsamp() to do the simulation, you would not consider upgrading to the newest arm (arm > 1.1-8).

A trick to speed up R matrix calculation

This is another example of why defaults matter a lot.

I got an email of Evan Cooch forward by Matt, saying that there exists a trick to speed up R matrix caculation. He found that if we replace the default Rblas.dll in R with the proper one. It can boost R’s speed in doing matrix caculation.

The file is here (This file only works under Windows). For Mac and Linux users, see here.

Here are the steps to replace the Rblas.dll file (for Windows users):

1. Check what kind of processor (CPU) your PC or laptop is using (My computer –> Property). Download Rblas.dll from the corresponding directory under http://cran.r-project.org/bin/windows/contrib/ATLAS/.

2. Go to your R directory tp locate where the Rblas.dll is, for example, c:/program files/R/R.2-7.0/bin. Rename it as Rblasold.dll so that if the new Rblass.dll doesn’t fit, you can use the old one by renaming it back.

3. Copy the new Rblass.dll you just download in this folder.

4. Restart R!

Here is an example to test:

require(Matrix)
set.seed(123)
X <- Matrix(rnorm(1e6), 1000)
print(system.time(for(i in 1:25) X%*%X))
print(system.time(for(i in 1:25) solve(X)))
print(system.time(for(i in 1:10) svd(X)))

Here is a test result on my machine (Intel Pentium (R) M process 1.73GHz with 1 GB RAM).

Default Rblas.dll

> print(system.time(for(i in 1:25) X%*%X))
   user  system elapsed 
 114.19    0.38  121.04 
> print(system.time(for(i in 1:25) solve(X)))
   user  system elapsed 
  87.03    0.28   89.31 
> print(system.time(for(i in 1:10) svd(X)))
   user  system elapsed 
 232.29    1.44  242.64 


New Rblas.dll

> print(system.time(for(i in 1:25) X%*%X))
   user  system elapsed 
  37.18    0.36   37.89 
> print(system.time(for(i in 1:25) solve(X)))
   user  system elapsed 
  30.62    0.56   31.78 
> print(system.time(for(i in 1:10) svd(X)))
   user  system elapsed 
 102.89    2.17  107.17 

Overall, R with the new Rblas.dll is twice faster than R with the old one. Now we wonder why this isn't a default design for R.