7
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • $\begingroup$ This looks like a mark-recapture problem, where you're basically guessing $P_{00}$, either through ML or Bayesian models. $\endgroup$
    – Roman Luštrik
    Commented Jan 16, 2013 at 11:05
  • $\begingroup$ @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? $\endgroup$
    – agstudy
    Commented Jan 16, 2013 at 12:07
  • $\begingroup$ @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. $\endgroup$
    – Gago-Silva
    Commented Jan 16, 2013 at 13:15

1 Answer 1

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
7
  • 1
    $\begingroup$ 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] $\endgroup$
    – David J. Harris
    Commented Jan 18, 2013 at 4:45
  • 2
    $\begingroup$ @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... $\endgroup$
    – Tomas
    Commented Jan 18, 2013 at 9:38
  • $\begingroup$ (+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. $\endgroup$
    – David J. Harris
    Commented Jan 18, 2013 at 18:23
  • $\begingroup$ 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)) $\endgroup$
    – Gago-Silva
    Commented Jan 24, 2013 at 12:41
  • 1
    $\begingroup$ @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.. $\endgroup$
    – Tomas
    Commented Jan 24, 2013 at 13:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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