If we have a $J$ groups of predictor, outcome (univariate) variable pairs, $$ \{(y_{j1}, x_{j1}) \ldots (y_{jn_j}, x_{jn_j})\}, \quad\text{for $j \in 1\cdots J$}, $$ a hiearchical linear regression model might be
$$ y_{ji} \sim \mathrm{N}(\alpha_j + \beta_j x_{ji}, \sigma^2),\\ \alpha_j \sim \mathrm{N}(\alpha_0, \tau^2_\alpha), \quad\beta_j \sim \mathrm{N}(\beta_0, \tau^2_\beta). $$
How would we make a hierarchical Gaussian process regression model that is analogous to this?
If we rewrite the linear model as $$ y_{ji} = (\alpha_0 + \beta_0 x_{ji}) + (a_j + b_j x_{ji}) + \epsilon_{ji},\\ a_j \sim \mathrm{N}(0, \tau^2_{\alpha}),\quad b_j \sim \mathrm{N}(0, \tau^2_{\beta}), \quad \epsilon_{ji} \sim \mathrm{N}(0, \sigma^2), $$ I guess that means that the Gaussian process analogy would be $$ y_{ji} = f^0_{ji} + f^j_{ji} + \epsilon_{ji},\\ f^0 \sim \mathrm{GP}(0, \kappa^0),\quad f^j \sim \mathrm{GP}(0, \kappa^1), \quad \epsilon_{ji} \sim \mathrm{N}(0, \sigma^2) $$ where $\kappa^0$, $\kappa^1$ are covariance functions. So, in other words, we have one general GP and one random-effects GP (as per comment by @Placidia). The general and group specific GPs are summed for each group, and this sum $f^0 + f^j$ would itself be a GP. The most recent edition of Gelman et al's Bayesian Data Analysis uses examples of additive GPs just like this.
Does that look right?
If this was being done using MCMC, I assume there would be nothing especially difficult about it.
By the way, there is this paper http://jmlr.csail.mit.edu/proceedings/papers/v13/park10a/park10a.pdf where the focus is hierarchical GPs. I have not properly studied it and it is not obvious if they are looking at exactly the same problem that I have described. But in any case, it does not look like it they think it is as easy as I just described it and also, it looks like a main concern is getting analytical approximations for the posterior predictive distributions.
