Visualization of the number of transitions between states I am currently developing a Markov model for ordinal data. In order to proceed with the modeling, I would like to check the distribution of the number of transitions per individual in my data set.
The data set looks like this (example for one individual, I have 150 individuals in total):
> (markov <- data.frame(ID=1,TIME=0:6,STATE=c(1,1,2,3,4,3,3)))
  ID TIME STATE
1  1    0     1
2  1    1     1
3  1    2     2
4  1    3     3
5  1    4     4
6  1    5     3
7  1    6     3

There are 4 different states (1-4). 
I have thought about a visualization like in this figure. So there should be one plot per possible transition in which the individuals are summarized, a plot per individual is not necessary. 
I would like to use R to do this. However, I'm a total beginner except for simple plots with ggplot. I have searched for this topic in this forum but the transition plots shown here are not suitable for me. How can I count the number of transitions between two states per individual and plot the results? A ggplot example would be fine, but a simple basic visualization is also totally alright.

EDIT: Thanks to Stephan. I have tried using the provided code. It seems to work; however, I don't get a distribution of transitions, just one vertical line per plot. I suppose, they represent the total number of a transition in the dataset. How can I solve this issue? Is it because I have directly imported and used the data from a csv file without further editing it?
Also, it seems that the code counts every transition within an individual. How can I edit the code that only "neighboring" transitions (e.g. transition between TIME=0 and 1, not TIME=0 and 2) are included?
Here's the current code; certainly there is something wrong with it.
nSTATE <- 4
data <- read.csv("~/.../.../transitiondata.csv")
transitions <- by(data,data$ID, # the file "data" is directly imported from a csv file
                  function(x)data.frame(ID=head(data$ID,-1),
                                         TIME=tail(data$TIME,-1),
FROM=head(data$STATE,-1),TO=tail(data$STATE,-1)))
transition_table <- lapply(transitions,function(x)with(x,table(FROM,TO)))
min_n_transitions <- min(unlist(transition_table))
max_n_transitions <- max(unlist(transition_table))
max_freq <- 100  # adapt by hand

par(mfrow=rep(nSTATE,2),mai=c(.4,.4,.4,.1))
for ( from in 1:nSTATE ) {
  for ( to in 1:nSTATE ) {
    foo <- sapply(transition_table,"[",from,to)
    hist(foo,freq=TRUE,
         #breaks=seq(min_n_transitions-.5,max_n_transitions+0.5),
         xlim=c(min_n_transitions,max_n_transitions),
         ylim=c(0,max_freq),xlab="n transitions",ylab="Frequency",
         main=paste("From",from,"to",to),las=1,col="lightgray")
    abline(v=quantile(foo,c(0.05,0.95)),lty=2,lwd=2)
  }
}


Here's the dataset.
 A: Is this what you had in mind?

I didn't quite understand what the solid vertical lines in the plot were, but the dashed ones give 5% and 95% quantiles as in the original.
Let's simulate some dummy data:
n_IDs <- 150
TIME <- 0:1000
n_STATEs <- 4

set.seed(1)
(markov <- data.frame(ID=rep(1:n_IDs,each=length(TIME)),
    TIME=rep(TIME,n_IDs),
    STATE=sample(1:n_STATEs,n_IDs*length(TIME),replace=TRUE)))

We first extract all transitions by judiciously using head() and tail() within each ID (by using by()):
transitions <- by(markov,markov$ID,
 function(xx)data.frame(ID=head(xx$ID,-1),TIME=tail(xx$TIME,-1),FROM=head(xx$STATE,-1),TO=tail(xx$STATE,-1)))

Next, we create transition matrices for each ID:
transition_table <- lapply(transitions,function(xx)with(xx,table(FROM,TO)))

Do take a look at these two data structures to understand what is happening.
We next make sure that all histograms align in terms of common horizontal and vertical axes, for comparability:
min_n_transitions <- min(unlist(transition_table))
max_n_transitions <- max(unlist(transition_table))
max_freq <- 50  # adapt by hand

I didn't find a good way to automate setting max_freq (for the upper end of the vertical axes). You may need to run the loop below twice, not plotting in the first run and just collecting the counts (which hist() returns). Or just set by hand, as I did.
Finally, plot in panels:
par(mfrow=rep(n_STATEs,2),mai=c(.4,.4,.4,.1))
for ( from in 1:n_STATEs ) {
    for ( to in 1:n_STATEs ) {
        foo <- sapply(transition_table,"[",from,to)
        hist(foo,freq=TRUE,
            # breaks=seq(min_n_transitions-.5,max_n_transitions+0.5),
            xlim=c(min_n_transitions,max_n_transitions),
            ylim=c(0,max_freq),xlab="",ylab="",
            main=paste("From",from,"to",to),las=1,col="lightgray")
        abline(v=quantile(foo,c(0.05,0.95)),lty=2,lwd=2)
    }
}

A: It is a bit unclear what you are doing with the lapply functions. You seem to be generating some big matrix for each of the ID's, but they contain the same content (creating the single bar). The example below makes the 4x4x150 transition table in a slightly simpler way.

A more important question... why do you use this graphic, or how can it be improved? (I would critique that the graphic is not working so well because the dynamic range of your data is not so great. You only have very few transitions per ID)
What is the information that you wish to convey? Or is it (just) a programming/debugging question?


nSTATE <- 4
data <- read.csv("~/Downloads/transitiondata.csv")


transitions <-  data.frame(ID=head(data$ID,-1),
                           TIME=tail(data$TIME,-1),
                           FROM=head(data$STATE,-1),
                           TO=tail(data$STATE,-1))
transition_table <- with(transitions,table(FROM,TO,ID))
min_n_transitions <- min(unlist(transition_table))
max_n_transitions <- max(unlist(transition_table))
max_freq <- 100  # adapt by hand

par(mfrow=rep(nSTATE,2),mai=c(.4,.4,.4,.1))
for ( from in 1:nSTATE ) {
  for ( to in 1:nSTATE ) {
    foo <- transition_table[from,to,]
    hist(foo,freq=TRUE,
         #breaks=seq(min_n_transitions-.5,max_n_transitions+0.5),
         xlim=c(min_n_transitions,max_n_transitions),
         ylim=c(0,max_freq),xlab="n transitions",ylab="Frequency",
         main=paste("From",from,"to",to),las=1,col="lightgray")
    abline(v=quantile(foo,c(0.05,0.95)),lty=2,lwd=2)
  }
}

