Markov Chain process - Missing values

I have a vector of sequences with presence (1) and absence (0), from were I have calculated the first order Markov Process.

This is how the data looks like:

dataset=c(NA,NA,0,0,0,0,1,0,1,1,0,1,1,NA,NA,NA,NA,NA,0,
1,1,1,1,1,1,0,NA,1,1,1,1,1,1,0,0,NA,NA,0,1,1,
NA,NA,0,NA,0,0,0,1,1,1,0,1,1,1,1,0,1,NA,0,1)


For thsi vector I have calculated:

• Probability of Presence (P1) and Absence (P0)
• Probability of having Presence followed by Presence (P11), Presence followed by Absence (P10), Absence followed by Presence (P01) and Absence followed by Absence (P00)

The transition probabilities were obtained by using a loop the checks for sequence of 2 values: to calculate P_00 I am using

   P_00 : dataset[j]==0 & dataset[j+1]==0 , etc.


These are the results:

P_0=0.3913043
P_1=0.6086957
P_0+P_1=1

P_00=0.1538462
P_01=0.2307692
P_11=0.4615385
P_10=0.1538462
P_00+P_01+P_11+P_10=1


The idea is to populate the NA with presence/absence sequences according to the Probability values obtained for the sequence. The problem is that for now I am not being able to find the best and more adequate process for this problem and I am a newbie with this type of problems. Even if there is already a package that can do this for me, I would prefer to understand how can I do it.

-

migrated from stackoverflow.comJan 17 '13 at 17:05

This question came from our site for professional and enthusiast programmers.

This looks like a mark-recapture problem, where you're basically guessing $P_{00}$, either through ML or Bayesian models. –  Roman Luštrik Jan 16 '13 at 11:05
@A.R which package do already what do you want to do ? and if I undesrtand you try to find a (best)method to replace you NA according to probs? –  agstudy Jan 16 '13 at 12:07
@agstudy for example with the HiddenMarkov package I can predict the states (using the probabilities obtained with the "observed states") and compare the result with the observed states. However I don't want to loose the "observed states" and just want to "predict" the NA states. –  A.R Jan 16 '13 at 13:15

This is best and most transparently done in bayesian modelling. Bayesian inference works using MCMC (Markov Chain Monte Carlo) simulations, and that's exactly what you need. Example model code in bugs (not tested!):

dataset[1] ~ dbern(P1) # you must somehow solve the first element
for (i in 2:n) { # markov chain - define how each value depends on the previous one
dataset[i] ~ dbern(p.presence[dataset[i - 1]])
}


dbern(p) stands for Bernoulli distribution with probability p. The array p.presence can be defined as (using P* variables as you defined them):

p.presence[0] = P01 / P0   # probability of presence, given absence in previous step
p.presence[1] = P11 / P1   # probability of presence, given presence in previous step


This way you can pass p.presence as input data (as in your example), but you could also let bugs to estimate it from the data!! (this is much more common and reasonable). You can of course get posterior distributions and MCMC simulation samples for the missing values (NA) in the dataset, and compute various other statistics on it.

-
I agree that MCMC is the way to go for this. But don't you also need to worry about transitions from the NAs to what was observed? I don't think your method "cares" if the third element of A.R.'s data set is a 1 or a 0, which seems wrong. In short, I think you need to add dependence on dataset[i + 1], not just dataset[i - 1] –  David J. Harris Jan 18 '13 at 4:45
@David, no! This is why bayesian modelling is so great. The NAs will be replaced by MCMC simulation samples using the model. Look at the rules I specified in the model (for loop) - this is exactly the markov chain principle, and that's all you have to care about! Don't be bothered with some details like missing values :-). We are in Bayesian world! You don't have to understand how it works at all, but if you want, read something about MCMC sampling or gibbs sampling... –  Curious Jan 18 '13 at 9:38
(+1) Ah, right, you're using BUGS. Yeah, it should take care of that for you. I was thinking of a different approach that only looked at the NAs. Nice answer. –  David J. Harris Jan 18 '13 at 18:23
Can't I do something like this: model { for (i in 1:N) { dataset[i] ~ dbern(p.presence[dataset1[i]]) dataset1[i]<-dataset[i]+1 } } list(N=10,dataset=c(0,0,1,1,0,NA,NA,0,NA,NA),p.presence=c(0.60,0.40))  –  A.R Jan 24 '13 at 12:41
@A.R. 1) and have you tried p.presence[dataset[i-1]+1]? And you mistakenly use i instead of i-1! 2) I don't know, I was at course with Marc Kery... maybe his book Introduction to WinBUGS for Ecologists could be of interest for you, mostly if you are ecologist.. –  Curious Jan 24 '13 at 13:41