A wraper function to convert coda files into a BUGS object
R graph with two y-axes
I’ve been asked how to do this several times, so I thought it might help to put an example online.
x <- 1:5
y1 <- rnorm(5)
y2 <- rnorm(5,20)
par(mar=c(5,4,4,5)+.1)
plot(x,y1,type="l",col="red")
par(new=TRUE)
plot(x, y2,,type="l",col="blue",xaxt="n",yaxt="n",xlab="",ylab="")
axis(4)
mtext("y2",side=4,line=3)
legend("topleft",col=c("red","blue"),lty=1,legend=c("y1","y2"))
R graph with two y-axes
I’ve been asked how to do this several times, so I thought it might help to put an example online.
x <- 1:5
y1 <- rnorm(5)
y2 <- rnorm(5,20)
par(mar=c(5,4,4,5)+.1)
plot(x,y1,type="l",col="red")
par(new=TRUE)
plot(x, y2,,type="l",col="blue",xaxt="n",yaxt="n",xlab="",ylab="")
axis(4)
mtext("y2",side=4,line=3)
legend("topleft",col=c("red","blue"),lty=1,legend=c("y1","y2"))
Drop unused factor levels
> a <- factor(letters)
> a
[1] a b c d e f g h i j k l m n o p q r s t u v w x y z
Levels: a b c d e f g h i j k l m n o p q r s t u v w x y z
## Now, even though b only includes five letters,
## all 26 are listed in the levels
> b <- a[1:5]
> b
[1] a b c d e
Levels: a b c d e f g h i j k l m n o p q r s t u v w x y z
## This behavior can be changed using the following syntax:
> b <- a[1:5,drop = TRUE]
> b
[1] a b c d e
Levels: a b c d e
Another way to deal with this is to use the dropUnusedLevels() command in the Hmisc library. The only issue here is that behavior is changed globally which may have undesired consequences (see the post listed above).
****UPDATE****
As Jeff Hollister mentions in the comments, there is another way to do this:
a<-factor(letters)
b<-factor(a[1:5])
Yet another way, if you are working with data frames that by default convert characters into factors, was suggested on r-sig-ecology by Hadley Wickham:
options(stringsAsFactors = FALSE)
a <-data.frame("alpha"=letters)
b<-a[1:5]
R package: codetools
Paper on the Gillespie Stochastic Simulation Algorithm in press
Just got news that my revisions to the reviewer’s comments on my paper GillespieSSA: Implementing the Gillespie Stochastic Simulation Algorithm in R were accepted. Hence, this paper is not officially in press in the Journal of Statistical Software. 
Here’s the abstract:
The deterministic dynamics of populations in continuous time are traditionally described using coupled, first-order ordinary differential equations. While this approach is accurate for large systems, it is often inadequate for small systems where key species may be present in small numbers or where key reactions occur at a low rate. The Gillespie stochastic simulation algorithm (SSA) is a procedure for generating time-evolution trajectories of finite populations in continuous time and has become the standard algorithm for these types of stochastic models. This article presents a simple-to-use and flexible framework for implementing the SSA using the high-level statistical computing language R and the package GillespieSSA. Using three ecological models as examples (logistic growth, Rosenzweig-MacArthur predator-prey model, and Kermack-McKendrick SIRS metapopulation model), this paper shows how a deterministic model can be formulated as a finite-population stochastic model within the framework of SSA theory and how it can be implemented in R. Simulations of the stochastic models are performed using four different SSA Monte Carlo methods: one exact method (Gillespie’s Direct method); and three approximate methods (Explicit, Binomial, and Optimized tau-leap methods). Comparison of simulation results confirms that while the time-evolution trajectories obtained from the different SSA methods are indistinguishable, the approximate methods are up to four orders of magnitude faster than the exact methods.
A preprint is available on my Publication page.
E-mail updates
