Markov chain in R

Question

Suppose we have a ten state system where an observation can enter the system in any one of the ten states with equal probability and move from the given state into a new state also with equal probability (the observation's new state isn't conditional on its previous state and the observation can remain in its current state). Additionally, an individual can "die" within any of the given ten states at any time including its current state. How exactly could this be setup in R or is this not even possible in R?

Why isn't `markovchain` package working for you? Alternatively, have tried to implement it step by step: http://alexhwoods.com/markov-chains-in-r/ ? — PSzczesny, Mar 07 '18 at 09:59

score 1 · Accepted Answer · answered Mar 07 '18 at 10:55

This sounds like a compartmental modelling problem. You could solve it with the SimInf package:

library(SimInf)

Define names of your 10 compartments:

compartments <- letters[1:10]

Define rates of enter and exit. In this case, as you suggested all enter events have the same rate to all compartments: k1 enter and k2 exit.

enterexit <- unlist(lapply(compartments, function(x){
    c(paste0("@ -> k1 -> ", x), paste0(x, "-> k2*", x, " -> @"))
}))

Define transitions between compartments. All get the same rate k3. So the individuals can just bounce around all 10 compartments at an equal probability:

transitions <- unlist(lapply(1:10, function(x){
    unlist(lapply(compartments[-x], function(y){
        paste0(compartments[x], "-> k3*", compartments[x], "->", y)
    }))
}))

Now define the initial state of the compartments. I'm going to put 0 in all ten compartments; you could also add some individuals to start if you wish:

u0 <-  data.frame(a = 0,
                  b = 0,
                  c = 0,
                  d = 0,
                  e = 0,
                  f = 0,
                  g = 0,
                  h = 0,
                  i = 0,
                  j = 0)

define the number of time steps for the model to run:

tspan = 1:100

Initialize the model:

model <- mparse(transitions = c(enterexit, transitions),
                compartments = compartments,
                k1 = 0.5,
                k2 = 0.5,
                k3 = 0.5)
model <- init(model, u0, tspan)

Run the model

ob <- run(model)

Plot it

plot(ob, N = TRUE)

Get a dataframe of the number of units in each compartment at each time step

df <- trajectory(ob)

thanks a lot for your example! However, what if I wanted to know the decay associated with each state i.e. the probability of how long the observations will stay within the current state for x amount of days? — menzd53, Mar 07 '18 at 21:10
The example above uses the [Gillespie algorithm](https://en.wikipedia.org/wiki/Gillespie_algorithm). The mean time that the individual spends in each compartment is the inverse of the sum of the rates that cause it to exit the compartment. In this case: 1/(k2+k3) ie. 1/(death_rate + movement_rate). The time to the next event is exponentially distributed. — trosendal, Mar 08 '18 at 15:23

Markov chain in R

1 Answers1