Estimating confidence intervals of a Markov transition matrix

Question

I have a series of n=400 sequences of varying length containing the letters ACGTE. For example, the probability of having C after A is:

and which can be calculated from the set of empirical sequences, thus

Assuming:

Then I get a transition matrix:

But I'm interested in calculating the confidence intervals for Phat, any thoughts on how I could I go about it?

Heinz

this is better suited on stats.stackexchange.com – Mcgaw

I think the reason for it being here is to attract Matlab specialists, programmers and enthusiasts. the stats forum question (stats.stackexchange.com/questions/64309/…) is not attacting any useful answers – Heinz

ok good point. I just figured you'd get better explanations there, I am not a statistician myself :) – Mcgaw

Answer 1

You could use bootstrapping to estimate confidence intervals. MATLAB provides bootci function in the Statistics toolbox. Here is an example:

%# generate a random cell array of 400 sequences of varying length
%# each containing indices from 1 to 5 corresponding to ACGTE
sequences = arrayfun(@(~) randi([1 5], [1 randi([500 1000])]), 1:400, ...
    'UniformOutput',false)';

%# compute transition matrix from all sequences
trans = countFcn(sequences);

%# number of bootstrap samples to draw
Nboot = 1000;

%# estimate 95% confidence interval using bootstrapping
ci = bootci(Nboot, {@countFcn, sequences}, 'alpha',0.05);
ci = permute(ci, [2 3 1]);

We get:

>> trans         %# 5x5 transition matrix: P_hat
trans =
      0.19747       0.2019      0.19849       0.2049      0.19724
      0.20068      0.19959      0.19811      0.20233      0.19928
      0.19841      0.19798       0.2021       0.2012      0.20031
      0.20077      0.19926      0.20084      0.19988      0.19926
      0.19895      0.19915      0.19963      0.20139      0.20088

and two other similar matrices containing the lower and upper bounds of confidence intervals:

>> ci(:,:,1)     %# CI lower bound
>> ci(:,:,2)     %# CI upper bound

I am using the following function to compute the transition matrix from a set of sequences:

function trans = countFcn(seqs)
    %# accumulate transition matrix from all sequences
    trans = zeros(5,5);
    for i=1:numel(seqs)
        trans = trans + sparse(seqs{i}(1:end-1), seqs{i}(2:end), 1, 5,5);
    end

    %# normalize into proper probabilities
    trans = bsxfun(@rdivide, trans, sum(trans,2));
end

As a bonus, we can use bootstrp function to get the statistic computed from each bootstrap sample, which we use to show a histogram for each of the entries in the transition matrix:

%# compute multiple transition matrices using bootstrapping
stat = bootstrp(Nboot, @countFcn, sequences);

%# display histogram for each entry in the transition matrix
sub = reshape(1:5*5,5,5);
figure
for i=1:size(stat,2)
    subplot(5,5,sub(i))
    hist(stat(:,i))
end

Mcgaw

This is fantastic! Thank you so much! Exactly why it's not in the stats forum! Giving bounty as soon as available! However, I do have a concern that actually row entries to Phat are not univariate normal distributions because they cannot vary individually. In fact each row forms a multivariate normal distribution don't you think? e.g. as one entry increases, at least one other in the same row must necessarily decrease to maintain sum(P,2)=1 – Heinz

Wow, lots of effort given here! – Arjan

Glad I could help. As for P_hat the rows should indeed sum to one (up to a certain precision). In the calculation function countFcn, I first accumulate the counts of co-occurrences from each sequence (using a less-known syntax of sparse function but we could have also used accumarray), then I divide by the row sums (the bsxfun call). One caveat: if a symbol doesn't occur in a sequence (say we only had [1 2 4 5]), then we could be dividing by zero, which would cause NaN values in the transition matrix. A common solution is to add +1 to the entire matrix before normalizing the rows – Mcgaw

@HughNolan Indeed this is the most helpful answer I could have hoped for! And worth every bit of rep! :) – Heinz

@Mcgaw Lovely! Is it possible to calculate the covariance of entries in the transition matrix P_hat? – Heinz

@HCAI: you can compute that from the output of bootstrp above as: C = cov(stat). Each row of stat is a flattened 5x5 transition matrix from one of the drawn samples. cov works on the columns of the input matrix, the result would be a 25x25 matrix corresponding to the covariances between each pair of entries of P_hat – Mcgaw

on another note, if you have the PCT toolbox, you can parallelize the computation of bootci and bootstrp (see the 'UseParallel' option in the docs) – Mcgaw

I'll have a look at parallelising this too. I'm considering the fact that actually my P_hat is sparse and as such does not well approximate P even with a bootstrap method. Do you think I should open a new question about Laplace smoothing to help this? E.g. people.cs.umass.edu/~dasmith/inlp/lect5-cs585.pdf – Heinz

@HCAI: I think Laplace smoothing is just that, adding a small lambda to smooth the probabilities (the +1 I mentioned above in the comments). This is a sort of regularization if you think about it, which should also prevent overfitting. – Mcgaw

... Now I dont know the about the domain background, but perhaps you could also try a more restricted MC structure which is not fully-connected, and not all transitions between states are permitted. For example a left-to-right model where state sequences are only allowed to move right (or stay the same) as time increases, but never go back.. As I said, I don't know if this is applicable for DNA sequences, just something to consider :) – Mcgaw

@Mcgaw Hmmm, this chain is actually attempting to represent nurses touching surfaces when in a patient's room (ACGTE attracts attention better). So technically all transitions are physically possible. However the number of observed sequences is clearly not big enough to account for this. So the $\lambda$ idea is worth pursuing... But technically this Markov chain's transition matrix can be represented by a multivariate normal distribution with mean(trans,3) and cov(stat)? – Heinz

@HACI: oops, of course (after all there is no "E" in DNA) :) Anyway, since we have discrete categorical variables here ("A,C,G,T,E" as opposed to continuous like "1.2, -3.4"), then the transition matrix is actually representing multinomial distributions one for each state, rather than a MVN distribution... – Mcgaw

Answer 2

Not sure whether it is statistically sound, but an easy way to get an indicative upper and lower bound:

Cut your sample in n equal pieces (for example 1:40,41:80,...,361:400) and calculate the probability matrix for each of these subsamples.

By looking at the distribution of probabilities amongst subsamples you should get a pretty good idea of what the variance is.

The disadvantage of this method is that it may not be possible actually calculate an interval with a desired given probability. The advantage is that it should give you good feeling for how the series behaves, and that it may capture some information that could be lost in other methods due to the assumptions that other methods (for example bootstrapping) are based on.

Chronicle

If you actually have multiple series of length n, you are probably better off calculating one probability matrix per series (rather than chopping up each series in smaller pieces) – Chronicle