Model formulation and parameter estimation

Model-based Geostatistics for Global Public Health

Emanuele Giorgi

Lancaster University

Overview of topics

• An overview of staionary Gaussian processes in 2 dimensions
• How to formulate and fit a geostatistical model
• How to interpret the results from a model fit

Stationary and Isotropic Gaussian Processes

• A Gaussian Process (GP) is a collection of random variables, any finite subset of which follows a multivariate normal distribution.
• A GP is stationary if its statistical properties do not change with location shifts.
• It is isotropic if its properties do not change with locations rotations.
• We denote a GP as \(S(x)\), where \(x\) represents location.
• For a stationary and isotropic GP, the correlation function is purely a function of the distance between locations.

The Matérn Process

• The Matérn covariance function is the most widely used in geostatistics.
• It is controlled by three parameters:
  • Variance (\(\sigma^2\)) – controls process magnitude
  • Smoothness (\(\kappa\)) – determines differentiability
  • Scale parameter (\(\phi\)) – governs spatial correlation decay
• The covariance function is given by: \[ C(h) = \sigma^2 \frac{2^{1-\kappa}}{\Gamma(\kappa)} \left( h/\phi \right)^\kappa K_\kappa \left(h/\phi \right), h = ||x-x'|| \] where \(K_\kappa\) is the modified Bessel function of the second kind.

matern_simulator.R

Special Cases of the Matérn Covariance

The Matérn family includes two important special cases:

• Exponential Covariance (\(\kappa = 1/2\))
  \[ C(h) = \sigma^2 \exp\left(-\frac{h}{\phi}\right) \]
  Properties:
  • Relatively rough process, non-differentiable
    
  • Continuous but not differentiable
    
  • Markovian (memoryless) property

• Gaussian (Squared Exponential) Covariance (\(\kappa \to \infty\))
  \[ C(h) = \sigma^2 \exp\left(-\frac{h^2}{2\phi^2}\right) \]
  Properties:
  • Infinitely differentiable
    
  • Produces very smooth realizations
    
  • Sometimes too smooth for physical processes

Practical Range of Spatial Correlation

• Definition: Distance (\(h^*\)) where correlation drops to 0.05 (5%)

  • Represents the effective spatial dependence range
  • Commonly used cutoff in geostatistics

• General Solution: \[ \rho(h^*) = 0.05 \] Solve for \(h^*\) to find practical range

• Exponential Covariance Example: \[ \rho(h^*) = \exp\left(-\frac{h^*}{\phi}\right) = 0.05 \] \[ \Rightarrow -\frac{h^*}{\phi} = \log(0.05) \approx -3 \] \[ \Rightarrow h^* \approx 3\phi \]

  • Practical range: \(3\phi\)
  • At \(h^* = 3\phi\), correlation drops to ~5%

Generalized linear geostatistical models

• Let \(S(x)\) be a stationary and isotropic Gaussian process, with Matérn correlation function.

• Let \(Z_i\) be independent identically distributed random variables.

• Let \(d(x_i)\) be a vector of covariates with associated regression coefficients \(\beta\).

• Let \(Y_i\) be a random variable that, conditionally on \(S(x_i)\) and \(Z_i\), follows an exponential family distribution (e.g. Gaussian, Binomial, Poisson).

• Let \(g(\cdot)\) be a link function.
  • \(\mathbb{E} \left[ Y_i \: | \: S(x_i), Z_i \right] = m_i \mu_i\)
    
  • \({\rm Var}\left[Y_i \: | \: S(x_i), Z_i\right] = m_i V(\mu_i)\)
  • \(g(\mu_i) = d(x_i)^\top \beta + S(x_i) + Z_i\)

Binomial and Poisson Geostatistical Models

• Binomial Geostatistical Model
  • \(Y_i \: | \: S(x_i), Z_i \sim \text{Binomial}(m_i, \mu_i)\)
  • Link function: Logit link \(g(\mu_i) = \log \left( \frac{\mu_i}{1 - \mu_i} \right) = d(x_i)^\top \beta + S(x_i) + Z_i\)
  • Variance function: \(V(\mu_i) = m_i \mu_i (1 - \mu_i)\)

• Poisson Geostatistical Model
  • \(Y_i \: | \: S(x_i), Z_i \sim \text{Poisson}(m_i \mu_i)\)
  • Link function: Log link \(g(\mu_i) = \log (\mu_i) = d(x_i)^\top \beta + S(x_i) + Z_i\)
  • Variance function: \(V(\mu_i) = m_i \mu_i\)

Monte Carlo Maximum Likelihood (MCML)

• Let \(W_i = S(x_i) + Z_i\) and \(W = (W_1, \ldots, W_n)\)

• The likelihood function for parameters \(\theta = (\beta, \sigma^2, \phi, \tau^2)\) is given by:

  \[ L(\theta) = \int N(W; D\beta, \Omega) f(y| W) dW \]

  where \(\Omega = \sigma^2 R(\phi) + \tau^2I_n\).

• We approximate this using Monte Carlo integration:

  \[ L_m(\theta) = \frac{1}{B} \sum_{j=1}^{B} \frac{N\left(W^{(j)}; D\beta, \Omega \right)}{N\left(W^{(j)}; D\beta_0, \Omega_0\right)} \] where \(W^{(i)}\) are sampled from the distribution of \(W\) given \(y\) using an MCMC algorithm.

MCML: Iterative Estimation

• Choose initial values \(\theta_0\) and generate samples \(W^{(j)}\), for \(j = 1,\ldots, B\).

• Estimate \(\theta\) by maximizing the Monte Carlo approximation, denoted by \(L_m(\theta)\), to the likelihood. Denote with \(\hat{\theta}_m\) this estimate.

• Update \(\theta_0 = \hat{\theta}_m\) and iterate until convergence.

• Standard errors are computed using the inverse Hessian of \(L_m(\theta)\) at \(\hat{\theta}_m\).

glgm_fit.R

The Nugget Effect as Small-Scale Spatial Correlation

• Interpretation of \(Z_i\):
  • Represents microscale variation (below measurement distance)
  • Captures measurement error or sub-resolution variability
• Example Process with Threshold Correlation: \[ \rho(h) = \begin{cases} \delta & \text{if } h \leq h_0 \\ 0 & \text{if } h > h_0 \end{cases} \]
  • Where \(\delta\) is small-scale correlation
  • \(h_0\) is resolution threshold
• Practical Implications:
  • If minimum distance between samples \(> h_0\):
    • Process appears as \(Z_i \sim N(0, \tau^2)\) (pure noise)
    • Spatial structure in \(Z_i\) becomes undetectable
  • Cannot disentangle \(Z_i\) from measurement error in the lineara geostatistical model