3
$\begingroup$

I would like to conduct a meta analysis in order to collate the information from a number of studies. The parameter of interest is a probability $\theta$.

In each of the studies, the observed data $y_i$ is seen as the result of a Poisson process:

$y_i|x_i\sim Poisson(f_i(x_i,\theta_i))$

where $x_i$ represents a set of observable variables.

Note importantly, that I have written $f_i$ to denote a function that depends on the particular characteristics (data length, data completeness) of study $i$. The exact relationship between $f$ and the data characteristics (not $x_i$ which represents the actual values of the data) is too long and complex to be able to write down in a particular form. I would like it to be determined at run-time via an auxillary function. As a simple example, imagine that the function is $\theta$ for a study of length 2, whereas $\theta^2$ for one of length 3 or 4, and $\theta^3$ for lengths longer than 4.

In reality, the function depends on a number of relatively complex conditional statements, and cannot be written in algebraic form.

In each study, I imagine that $\theta_i$ is drawn from a hyper-distribution:

$Log(\frac{\theta_i}{1-\theta_i}) \sim N(\mu,\sigma^2)$

where the hyperprior on $\mu$ is uniform on $[0,1]$, and is flat on $log(\sigma)$ space.

The main issue I have in implementing the above, is that the function $f_i$ depends on the particular data for that individual study. There are too many studies for me to code by hand a separate function for each of them, this will need to be done dynamically. The function is relatively complex, but is determined completely by the data set characteristics.

Does anyone know how I could do this in BUGS, JAGS or STAN? Is it possible to call a function which will do this from within BUGS, JAGS or STAN?

Many thanks,

Ben

Improve this question
$\endgroup$
5
  • $\begingroup$ The statement "The function...is determined completely by the data set characteristics." makes me think that $f$ should not have an $i$ subscript, but instead $f$ only depends on $x_i$ (the data set characteristics). $\endgroup$
    – jaradniemi
    Apr 16, 2015 at 14:43
  • $\begingroup$ @jaradniemi - thanks for your comment. However, I like to disaggregate a couple of different effects. $x_i$ represents the values of the data, whereas I want $f$ to vary according to the data dimensions. This is why I use $f_i$ to account for the latter. I will edit the above to make this clear. Best, Ben $\endgroup$
    – ben18785
    Apr 16, 2015 at 15:09
  • $\begingroup$ The example function you give is just $(1-\theta)^{length-1}$. In any case, if you can parameterize your function in terms of $z_i$ (the data set characteristics), then you can program it. $\endgroup$
    – jaradniemi
    Apr 16, 2015 at 15:27
  • $\begingroup$ @jaradniemi - sorry if my example was not taken for what it was. In this simple case, it would be very easy to write. Suppose that instead the function is a complex-control-structure based function, can I still call it from within BUGS to find the exact functional form? $\endgroup$
    – ben18785
    Apr 16, 2015 at 15:38
  • $\begingroup$ @jaradniemi - I have now changed the example I give, to yield one that yields more similar results to the type of function I would like to call. Hope that clarifies it. Best, Ben $\endgroup$
    – ben18785
    Apr 16, 2015 at 15:47

1 Answer 1

Reset to default
3
+50
$\begingroup$

I am assuming that you have a function that can be written in a programming language such as R or C, but is not translatable directly to JAGS/BUGS. In this case you have two options:

  1. If your function can be discretised into n-dimensional space, it is possible to pre-calculate a grid of evaluated points in R and then pass these evaluated values AS DATA into JAGS/BUGS. You can then simply use the sampled parameter values to index the array and find the corresponding value of the evaluated function. This has the advantage that computation is moved out of JAGS/BUGS, so the model itself will be much faster. However, if you have a lot of dimensions (i.e. parameters) and/or don't want to discretise your parameter space this is a non-starter.

  2. You can write your own function in C++ and load this into JAGS as an external module (it is easiest in JAGS if you know C++). A very useful tutorial on how to do this is given at http://www.ncbi.nlm.nih.gov/pubmed/23959766 - but you can avoid a lot of the configure setup by putting the JAGS module inside an R package. For an example of the same module both inside an R package and as a standalone JAGS module you can look at http://runjags.sourceforge.net

I think there were plans to add customised functions to Stan but I'm not sure if this has been done yet - it is a while since I last played with Stan!

Hope that helps,

Matt

Improve this answer
$\endgroup$
3
  • $\begingroup$ Hi Matt, thanks for your message. Unfortunately, the function can't really be discretised into n-dimensional space. The issue is that it depends not on the $x_i$ values, but on its dimensions. I also don't think that the second option will work either. Basically, I think I want to do a Poisson regression each time, with different independent variables. Best, $\endgroup$
    – ben18785
    Apr 20, 2015 at 13:14
  • $\begingroup$ It is not clear to me why you think the second option won't work? If you can define your function of interest in C++ as a function of your data and parameters then it should work. A lambda value calculated from this function can then be used for the Poisson regression (or whatever) in JAGS. Unless you mean that the dimensions of the parameters need to be sampled as well as the parameter values, in which case you need something like reversible jump MCMC. Matt $\endgroup$
    – Matt Denwood
    Apr 21, 2015 at 8:20
  • $\begingroup$ Hi Matt, thanks for your help here. I ended up using Stan in the end. It seems by far the most flexible language (as well as its speed advantage) for MCMC. Best, Ben $\endgroup$
    – ben18785
    Apr 25, 2015 at 15:54

Your Answer

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

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