2
$\begingroup$

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.

$\endgroup$

closed as off-topic by Martijn Weterings, mkt, Michael Chernick, whuber Feb 13 at 22:46

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – Martijn Weterings, mkt, Michael Chernick, whuber
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ (1) How many states do you have? Any more than 5 will make the visualization hard to read, and you may need something different. (2) How do you want to treat your 150 individuals? Do you want separate plots per individual, or should they all be thrown together, or something else? (3) Do you want something in ggplot, or would base graphics be fine? (4) This may be migrated to SO for being off-topic here. $\endgroup$ – Stephan Kolassa Feb 6 at 9:08
  • $\begingroup$ I have edited my answer. Sorry for not being clear in the first post. $\endgroup$ – asher Feb 6 at 9:26
  • $\begingroup$ Thank you for editing, and no need to apologize. Asking for clarification is what comments are for. $\endgroup$ – Stephan Kolassa Feb 6 at 9:38
  • $\begingroup$ Hm. Can you (1) edit your code to show the precise read.table() function you use to read your csv file, and (2) upload the csv file somewhere we can access it? There may be an issue in whether states etc. are coded as numerical or as factors. The code should only count transitions between adjacent times if (!) your dataset is sorted by time. Is yours? $\endgroup$ – Stephan Kolassa Feb 6 at 17:35
  • $\begingroup$ Thanks for your reply. I have edited the code. Actually, I have used the "Import Dataset" function in RStudio. I have uploaded the file as well. The dataset is sorted by ID, but within one individual, time is in chronological order. $\endgroup$ – asher Feb 6 at 18:00
0
$\begingroup$

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?


enter image description here

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)
  }
}
$\endgroup$
  • $\begingroup$ Dear Martijn, this works great, thanks a lot! Well, I am modeling transitions between these states. Since I had problems with the inclusion of between-subject variability I wanted to check what the distribution of transition numbers looks like. In fact, this plot supports my assumption that the data is not the best for my current modeling approach since the numbers are very low. In later stages of modeling, I would use this kind of plot to compare the distribution of my simulations with the observed number of transitions as some kind of internal validation. $\endgroup$ – asher Feb 13 at 10:02
  • $\begingroup$ @asher I find this a bit strange way to compare the cases. Those transition frequencies will eventually be binomial distributions multiplied with the probability to be in a particular state (if only the ID's would have equal amounts of transitions, so they are actually binomial distributions compounded with a variable parameter n for the amount of times the ID is in a particular state). Why don't you just compute an estimate for the 16 transition probabilities and compare those. Or do you want to test whether the observed transitions are actually corresponding to a Markov model or not? $\endgroup$ – Martijn Weterings Feb 13 at 11:01
  • $\begingroup$ That's an interesting point, thanks. But yes, I just want to look if the model fits to the observed transitions (see also the link to the figure in my first post). I have also planned other validation methods for my model (e.g. different kinds of visual predictive checks). $\endgroup$ – asher Feb 13 at 12:33
3
$\begingroup$

Is this what you had in mind?

histograms

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)
    }
}
$\endgroup$
  • $\begingroup$ Dear Stephan, this looks really great. I will take some time to have a closer look into this in order to understand the different steps. The solid lines in the example are the number of observed transitions whereas the bars represent the model-predicted number of transitions. Thanks a lot! $\endgroup$ – asher Feb 6 at 11:08
  • $\begingroup$ I have further edited my first post. It seems to work, there are however some questions I would like to ask. $\endgroup$ – asher Feb 6 at 15:10

Not the answer you're looking for? Browse other questions tagged or ask your own question.