Building a table of probability transition matrices with plyr

Question

I'm trying to model a system of continuous time Markov chains where in different time intervals I have different rates.

I build a rate matrix for each time period like this

make.rate.matrix <- function(c1, c2, m12, m21) {
  matrix(
    c(# State 1: lineages in different populations
      -(m12+m21), m21, m12, 0,
      # State 2: both lineages in population 1
      2*m12, -(2*m12+c1), 0, c1,
      # State 3: both lineages in population 2
      2*m21, 0, -(2*m21+c2), c2,
      # State 4: coalesced (catches both populations; absorbing)
      0, 0, 0, 0),
    byrow = TRUE, nrow=4)
}

(if you are interested it is modelling the coalescence density in a two-deme system with migration)

The rates, the cs and ms, differs in different time periods, so I want to build a rate matrix for each time period and then a transition probability matrix for each period.

With two periods I can specify the rates like this

rates <- data.frame(c1 = c(1,2), c2 = c(2,1), m12 = c(0.2, 0.3), m21 = c(0.4, 0.2))

and I want to use the first rates from time 0 to t and the second set of rates from time t to s, say.

So I want to have a table of rate matrices for the first and second period, and probability transition matrices for moving from state a to b through the first and second period.

mlply(rates, make.rate.matrix)

gives me a list of the two rate matrices, and if I want a table where I can easily look up the rate matrices, I can do something like

> xx <- array(unlist(mlply(rates, make.rate.matrix)), dim=c(4,4,2))
> xx[,,1]
     [,1] [,2] [,3] [,4]
[1,] -0.6  0.4  0.2    0
[2,]  0.4 -1.4  0.0    1
[3,]  0.8  0.0 -2.8    2
[4,]  0.0  0.0  0.0    0
> xx[,,2]
     [,1] [,2] [,3] [,4]
[1,] -0.5  0.2  0.3    0
[2,]  0.6 -2.6  0.0    2
[3,]  0.4  0.0 -1.4    1
[4,]  0.0  0.0  0.0    0

I can then get the probability transition matrices like

> library(Matrix)
> t <- 1; s <- 2
> P1 <- expm(xx[,,1] * t)
> P2 <- expm(xx[,,2] * (s - t))

but I somehow cannot figure out how to get a table of these like I can get for the rate matrices.

I feel that I should be able to get there with aaply, but I am stomped as to how to get there.

How do I get a table P, where P[,,1] is P1 and P[,,2] is P2?

Building a table of probability transition matrices with plyr

0 Answers0