Bayesian meta analysis: implementation in BUGS/JAGS/STAN 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
 A: 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:


*

*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.

*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
