Branching process HMM in NIMBLE

[Adam Howes]

Hi all,

Nice to meet you, I'm Adam, and new to using NIMBLE.

At the moment I'm working on model for the amplfication of DNA during a qPCR reaction in collaboration with a wet lab. The model that I would like to specify is a hidden Markov model with a discrete Binomial branching process latent space. Something like:

n_0 ~ Binomial(N_0, r)
n_{c + 1} = n_c + Binomial(n_c, E)
f_c = alpha * n_c + eps_c
eps_c ~ Normal(0, sigma_f^2)

where:

* N_0 is the initial number of plasmids

* n_0 is the initial number of DNA

* r is the proportion of DNA which are succesfully released from the plasmid

* n_c is the number of DNA at qPCR cycle c
* E is the "efficiency" (it would be in [0, 1])
* f_c is the fluorescence
* alpha is the fluoresence produced by a single DNA
* eps is Gaussian noise

Do you think this type of model would be possible to fit in NIMBLE? I had thought it would be a good choice based on my research for being able to handle discrete parameters, and popularity for fitting HMMs in ecology.

I'm a bit rusty on HMM terminology, but I'm most interested in "filtering" (inference on the n_c) and "smoothing" (inference on f_c [this more so than filtering, apart from n_0]) and not at all on predicting. But the most important parameters that I'd like to infer are E, the efficiency, n_0, the initial number of DNA going into the reaction, and N_0, the initial number of plasmids. Experiments can be designed to help if some of these things are not directly identifiable.

Do you think it would be possible to specify and fit this model using NIMBLE? Are there any suggested resources I should read? (I have the Efficient Markov chain Monte Carlo sampling for hierarchical hidden Markov models open now!). Would be very appreciative of any other pointers if you have them. :)

Best wishes,
Adam

[Perry de Valpine]

Hi Adam,

My first thought would be to set this up directly rather than using HMM-specific tools.  I think the model code would be something like this:

N_0_cont ~ dunif(0, 1000) #  prior for N_0, over whatever range seems relevant 

N_0 <- round(N_0_cont)

r ~ dunif(0, 1) # prior for r

n[1]~ Binomial(N_0, r) # indexing starts at 1, not 0, so this is your n_0

for(i in 2:C_max) {

   n_new[i] ~ dbinom(n[i-1], E)

   n[i] <- n[i-1] + n_new[i]

   f[i] ~ dnorm(alpha * n[i], sd = sigma_f)

}

sigma_f ~ dunif(0, 10) # flat prior over whatever range seems relevant

E ~ dunif(0, 1) # flat prior

alpha ~ dunif(0, 100) # or alpha ~ T(dnorm(0, sd = 100), 0, Inf) to get a half-normal

If you want to see more about direct (and discrete) HMM support, take a look at the vignette for nimbleEcology .  The HMMs there are quite general.  But if the N_0 and n_c values are large, you'd need to set up large transition and observation matrices, which you'd have to populate with values of interest, and you'd have to discretize f[c].  But if the above code makes sense, maybe try that first.

Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/0ddcc7fe-7de6-4b5d-b661-33930aba84ben%40googlegroups.com.