How to estimate uncertainty in Markov chain simulations

Question

Consider how I fit a Markov chain to my data with R:

library(markovchain)
library(dplyr)
library(ggplot2)
library(data.table)

#Data
A<-structure(c("sunny", "sunny", "sunny", "sunny", "sunny", "sunny", 
"rain", "cloudy", "rain", "cloudy", "sunny", "cloudy", "cloudy", 
"cloudy", "cloudy", "sunny", "sunny", "sunny", "sunny", "rain", 
"sunny", "rain", "sunny", "sunny", "rain", "cloudy", "rain", 
"sunny", "sunny", "cloudy", "rain", "cloudy", "rain", "sunny", 
"rain", "rain", "rain", "sunny", "cloudy", "cloudy", "cloudy", 
"cloudy", "cloudy", "cloudy", "sunny", "cloudy", "rain", "rain", 
"cloudy", "cloudy", "sunny", "sunny", "cloudy", "cloudy", "cloudy"
), .Dim = c(5L, 11L), .Dimnames = list(NULL, c("time1", "time2", 
"time3", "time4", "time5", "time6", "time7", "time8", "time9", 
"time10", "time11")))

#estimate transition matrix
B<-markovchainFit(data=A,name="weather")
mcWeather<-B$estimate

##### Do the forecasting over time and find uncertainty due to small sampling size

KKK<-list()
for(j in 1:10000){
LL<-list()
for(i in 1:5){
  LL[[i]]<-data.frame(cat=rmarkovchain(n = 10, object = mcWeather, t0 = "sunny",include.t0 = TRUE),index=i,time=1:11)  
}

LLL<-rbindlist(LL)
KKK[[j]]<-LLL %>% group_by(time,cat) %>% summarize(freq=n()/i)

KKK[[j]]$perm=j
}

KOO<-rbindlist(KKK)

KKX<-KOO %>% group_by(time,cat) %>% summarize(mean=mean(freq),lq=quantile(freq,0.025),up=quantile(freq,0.975))

# Plot results
ggplot(KKX,aes(x=time,y=mean,color=cat))+geom_line()+ geom_ribbon(aes(ymin=lq, ymax=up),color="grey",alpha=0.3)+facet_wrap(~cat)  

Some more detail to the code: I have 5 individuals that show a sequence of states over time, which can be expressed as a Markov chain.

I fit a Markov chain model to my data to obtain my transition matrix. With this I can now forecast the expected probabilities or expected distribution of my states over time. E.g., via:

W0<-t(as.matrix(c("cloudy"=0,"rainy"=0,"sunny"=1))) #start category sunny
for (time in 1:10){
W0 * (B$estimate ^ time)
}

But if I repeated my experiment with another 5 individuals I would not necessarily observe my expected distribution of the states over time, because this can be seen as 5 random draws of my Markov chain. These are not enough samples to hit the expected distribution perfectly. With this simulation I try to account for that by 10000 times draw sequences for 5 individuals and calculate the uncertainty of the distribution of my states over time. With this I can account for the uncertainty due to small sampling size (low number of individuals) and better compare different experiments all based on 5 individuals.

So with this code I have to some extent accounted for uncertainty of the small sampling size of 5. (See how confidence increases when changing to for(i in 1:5000){.. right? Or is this way wrong already?)

Now my question is: Does my estimated transition matrix—the one I used for the simulation above—not already have some uncertainty? Each entry of the transition matrix is estimated from very few observations (sequences of 5 individuals) as well. I saw the function markovchainFit() includes confidence interval estimates for the matrix entries, but I don't know how to link and combine this to the uncertainty estimation I already have done, so that in the end I get a 'global' estimation of the uncertainty in the forecast including (1) the uncertainty in the probability estimates of the transition matrix I fitted to my data and (2) the uncertainty I simulate above due to small sample size.

Answer 1

(1) the uncertainty in the probability estimates of the transition matrix I fitted to my data and (2) the uncertainty I simulate above due to small sample size.

(1) I believe markovchainFit assumes that the parameters of estimated transition matrix follows normal distribution. Because the returned confidence intervals match the values of mean +- 1.96 * SE. Where 1.96 is the critical value for normal distribution at 95% confidence level.

# upper CI provided by the model 
B$upperEndpointMatrix
#        cloudy    rain   sunny
#cloudy 0.67698 0.60005 0.52126
#rain   1.00000 0.49190 0.36706
#sunny  0.56847 0.41684 0.78315

# calculated as for normal distribution
B$estimate@transitionMatrix + 1.96 * B$standardError
#        cloudy    rain   sunny
#cloudy 0.67698 0.60006 0.52126
#rain   1.04182 0.49191 0.36707
#sunny  0.56847 0.41684 0.78316

EDIT

As you pointed out, to sample individual entries of the transition matrix independently from normal distribution is problematic as row sums can be larger than one.

To stochastically generate a transition matrix, we can use bootstrap method provided in markovchainFit.

B <- markovchainFit(data=c(A[1,], NA, A[2,], NA, A[3,], NA, A[4,], NA, A[5,]),
                 name="weather",
                 method='bootstrap',
                 nboot=1000) #number of models generated

#Check the row sums are indeed 1s
rowSums(B$bootStrapSamples[[1]])
#cloudy   rain  sunny 
#     1      1      1 

(2) The uncertainty due to small sample size could be viewed as binomial distribution with small sample size n. In your example, the expected distribution for "sunny" at step 10 with 5 observation is B(n=5, p=0.30408) / 5.

W0 * (B$estimate ^ 10) # prob of sunny at step 10
#  cloudy    rain   sunny
# 0.42746 0.26846 0.30408

The mean of B(n, p)/n is p, and the variance is p(1-p)/n. As n becomes larger, the variance reduces.

So, to simulate with both (1) and (2) uncertainties, we may first sample a transition matrix from bootstrapping normal distribution with mean and SE provided by the model. And followed by your method for random Markov chain sampling.

for(j in 1:10000){
LL<-list()
for(i in 1:5){
  mcWeather2 <- mcWeather
  #mcWeather2@transitionMatrix <- mcWeather2@transitionMatrix + 
  #                               rnorm(9, sd=B$standardError) # add parameter uncertainty
  #mcWeather2@transitionMatrix[mcWeather2@transitionMatrix>1] = 1
  #mcWeather2@transitionMatrix[mcWeather2@transitionMatrix<0] = 0

  #EDIT: using bootstrap method for parameter uncertainty
  mcWeather2@transitionMatrix <- B$bootStrapSamples[[sample(length(B$bootStrapSamples), 1)]]
  
  LL[[i]]<-data.frame(cat=rmarkovchain(n = 10, object = mcWeather2, t0 = "sunny",include.t0 = TRUE),index=i,time=1:11)  
}