\(\underbrace{P(parameters|Data)}_\text{Posterior}\propto \underbrace{P(Data|parameters)}_\text{Likelihood}\underbrace{P(parameters)}_\text{Prior}\)

• For modelling puposes we consider a latent species response variable

• Each latent species variable has a functional relationship with elevation:

\[{\bf s}_j = f_{\theta_j}(\text{elevation}) + {\bf\epsilon}_j\]

• f is a penalised-spline function and \({\bf \theta}_j\) is a vector of parameters controlling the shape of the spline for species \(j\).

• \(\epsilon_{ij} \sim N(0, \sigma_j^2)\) and is added to account for overdispersion (i.e, the data here are likely to exhibit more variation than may be capturped with the data model).

• The observed species abundances are multinomial

  \[(y_{i1},....,y_{iM}) \sim Multi(p_{i1},....,p_{iM}, N_i)\]

• The probabilites of the multinomial distribution are estimated as a function of the latent species variable via a softmax transformation where

\[{p}_{ij} = \frac{e^{s_{ij}}}{\sum_{j=1}^M e^{s_{ij}}}\]

\(\underbrace{P(parameters|Data)}_\text{Posterior}\propto \underbrace{P(Data|parameters)}_\text{Likelihood}\underbrace{P(parameters)}_\text{Prior}\)

A Case Study from New Jersey, USA

Consider the multivariate normal distribution for some k-dimension random vector y \[y \sim MVN(0, \Sigma)\] The matrix Sigma currently looks like this: \[\Sigma = \left[ \begin{array}{cccc} \sigma^2 & 0 & \ldots & 0 \\ 0 & \sigma^2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \sigma^2 \end{array} \right]\]

• The key idea in Gaussian Processes is to change the off-diagonals of \(\Sigma\) so that the you get higher correlations between \(y_i\) and \(y_j\) when \(t_i\) and \(t_j\) are close.

• A zero-mean Gaussian process is defined as follows:

\[y \sim MVN(0,K)\]

• one commonly used option is the squared exponential covariance function:

\[K_{i,j} = \sigma_g^2exp\bigg(-\rho^2(t_i-t_j)^2\bigg)\]

• \(\sigma_g^2\) should be high for functions that cover a broad range on the y-axis

• if \(t_i\) is distant from \(t_j\) then \(K_{i,j} \approx 0\)

• the effect of \(t_i-t_j\) on \(K_{i,j}\) will depend on \(\rho\)

• We assume a GP model for the rate of sea-level change

\[\omega(t) \sim MVN(0, K)\]

Where \(K_{i,j} = \sigma_g^2exp\bigg(-\rho^2(t_i-t_j)^2\bigg)\)

• Then the sea-level process is obtained by integrating the rate process

\[s(t) = \int_0^t \omega(u) du\]

• The observed sea-level values are assumed to be normally distributed

\[y_i \sim N(s(t_i), \sigma_i^2)\] where \(\sigma_i^2\) will capture the variation of the observed sea-level values around the mean.

N Cahill, A C. Kemp, B P. Horton and A C. Parnell. A Bayesian Hierarchical Model for Reconstructing Relative Sea Level: From Raw Data to Rates of Change.Climate of the Past} , 12(2):525-542, 2016.

N Cahill, A C. Kemp, B P. Horton and A C. Parnell. Modeling sea-level change using errors-in-variables integrated Gaussian processes. Annals of Applied Statistics, 9(2): 547-571, 2015.