Weighted model calibration with spatial conditional information

Cost functions which are often used for model calibration can be expressed as:

where $\mathbf{y}=(y_{i},\dots,y_{n})$ is the vector of $n$ observations and $\hat{\mathbf{y}}=(\hat{y}_{i},\dots,\hat{y}_{n})$ denotes their corresponding model estimates. The function $g$ describes the way the individual model residuals (the difference between the modelled and observed values) contribute to the cost function. In this formulation, since the cost function is a simple sum of these contributions from the $n$ observations, the observations are thought to provide equal and independent information. This assumption is implicit for mean square error (MSE) where $g(\hat{y}_{i}-y_{i})=\frac{1}{n}(\hat{y}_{i}-y_{i})^{2}$ and mean absolute error where $g(\hat{y}_{i}-y_{i})=\frac{1}{n}|\hat{y}_{i}-y_{i}|$.

In the case of MSE, the cost function is directly proportional to the negative log-likelihood of independent Gaussian random variables (less the scalar constant $n\log(\sigma\sqrt{2\pi})$):

Hence, minimising MSE is akin to maximising the likelihood. Here, the random variables are $Y_{i}\sim N(\hat{y}_{i},\sigma^{2})$ for some standard deviation $\sigma>0$.
Instead of independent Gaussian likelihoods, a full spatial likelihood can be used to express the spatial dependence in the model residuals and maximum likelihood estimation through for example, kriging can be used to estimate model parameters (Chiles & Delfiner, 1999). The likelihood, however, now takes the form of a multivariate Gaussian density involving the inverse and determinant of a $n\times n$ spatial covariance matrix. For large datasets, exact analytical methods become infeasible, especially when we need to estimate multiple parameters as in the case of process-based models.
One way to approximate a dependent likelihood is to weight independent likelihoods according to some measure of their dependence (Varin et al., 2011). A weighted cost function can be expressed as:

where $w_{i}$ denotes a calibration weight related to spatial dependence. For example, weighted MSE can be defined as:

In kriging, the variogram is often used to describe spatial dependence. Let $e_{i}=\hat{y}_{i}-y_{i}$, the model residual corresponding to observation $i$ at site $s_{i}$. The variogram is defined as half the expected square difference between model residuals at different locations:

There is a one-to-one relationship between the variogram and the spatial covariance:

Under the common assumptions of stationarity and isotropy, the variogram becomes a function of the distance between observations, $\gamma(s_{1},s_{2})=\gamma(h)$ where $h=|s_{1}-s_{2}|$. An example is the exponential variogram:

where $\sigma^{2}$ denotes the partial sill (the amount of spatially correlated variance), $\phi$ is a range parameter (controlling the rate of correlation decay with distance) and $\tau^{2}$ is the nugget variance (the amount of spatially independent variance or measurement error). Figure 1 shows two exponential variograms with the same nugget variance and partial sill but different range parameters. The larger range parameter leads to a slower decay in spatial correlation.

In our analyses, we will use Matérn variograms which commonly used in Spatial Statistics and which extend the exponential variogram towards greater flexibility in terms of the smoothness of the generated processes:

where $\Gamma$ denotes the Gamma function, $\kappa$ is a smoothness parameter and $K_{\kappa}$ is the modified Bessel function of the second kind.

Next, we introduce the concept of spatial conditional information (CI). This is related to the full spatial likelihood which involves spatial covariance matrices related to variograms. Based on transformations of CIs which are chosen to obtain desirable properties, we develop a weighting scheme for model calibration. Inference results from the weighted calibration will be compared with unweighted calibration in later Sections.

We consider both random spatial sampling and clustered sampling. In line with our tephra data, our base setting involves 80 observations for each simulation. Figure 2 shows the four sampling schemes we consider: 80 randomly chosen locations, 40 randomly chosen locations, 80 locations with relative positions fixed to be the same as that in the tephra load data from the 2014 Kelud eruption and 80 locations exhibiting higher clustering than our case study data. To generate the latter, we used k-means clustering to divide the original locations into five clusters (five being the number chosen by the elbow method) and halved the distance of each observation to its associated cluster centre.
To generate Gaussian process values with spatial dependence at the chosen locations, we use a Matérn variogram with nugget $\tau^{2}=1$, sill $\sigma^{2}=4$ and range parameter $\phi=2.5$. These relative values are similar to what we obtain for the case study in Section 4.1. Note that the smoothness parameter $\kappa$ is set to $1$ which is a common choice which allows for slightly smoother surfaces than an exponential and the default value in the R package R-INLA (R-INLA, 2022; Rue et al., 2009). We generate $100$ sets of simulated values per sampling scheme and estimate the process mean for each by optimising the full spatial likelihood (via the ‘likfit’ function of the geoR package) as well as weighted and unweighted MSE.
Figure 3 shows the box plots of the process mean estimates from the 40 and 80 randomly sampled locations for the three calibration methods. With less observations, there is higher bias and spread in the estimates as seen in the larger differences between the bold line denoting the median values from the dashed line denoting the true mean of zero as well as the longer widths of the box plots. While the WMSE calibration has slightly lower bias in the case with 40 observations than the other two methods, its benefit is hard to ascertain with 80 observations. This is likely because, by random sampling, we have a good representation of the study region and process values, and this is even more so with more observations considered.
Holding the number of observations fixed at 80, we examine the effect of clustering our observations. Figure 4 shows the box plots of the process mean estimates from the 80 locations corresponding to those of the Kelud case study and after increasing the extent of clustering by halving the distance of the observations to their cluster centres. Under increased clustering, we see that the bias of the unweighted calibration increases while the estimates from the weighted calibration remain unbiased. The precision of all three methods decrease with increased clustering but based on the width of the boxes, unweighted calibration is affected the most. The same conclusions can be drawn when we compare the box plots corresponding to the highly clustered observations to that for 80 randomly selected locations in Figure 3.

From the previous section, we saw that weighted calibration is important for clustered observations. In this section, we vary the spatial dependence of the Gaussian process to see its effect on the mean estimates. We do this by increasing the proportion of the nugget variance to the sill from $0.2$ to $0.4$ (hence increasing random noise and reducing spatial dependence from the ‘Mid’/base case) and increasing the spatial range parameter $\phi$ from $2.5$ to $3.5$ (hence increasing dependence). In Figure 5, the former case is represented by ‘Low’ and the latter by ‘High’. Generally, we see that when we increase the spatial dependence in the Gaussian process, the precision of the mean estimates decreases (the boxes are wider). Nevertheless, the estimates from the weighted calibration remain relatively unbiased under Mid-High spatial dependence unlike those from the unweighted calibration and exhibit higher precision.

We conduct weighted and unweighted calibration of the Tephra2 model using the tephra load data from the 2014 Kelud eruption. Note that calibration was not conducted using a full spatial likelihood unlike in the simulation experiements for the Gaussian process in Section 3 because of the infeasibility of this approach for large environmental datasets. The fitted variogram for the model residuals from the unweighted calibration is shown in Figure 6. The increase from a nugget value of $\tau^{2}=0.0106$ to the sill of $0.0512$ indicates that there is spatial correlation present since the expected difference in residuals increases with increasing distance between observations.
We use the fitted variogram to construct the spatial covariance matrix required to compute the CI weights and conduct the weighted calibration. The computed weights are shown on a map in Figure 7. We see that clustered observations are given lower weights. Table 1 shows the parameter estimates from the weighted and unweighted calibration. We see that the difference in estimates ranges from -71.587% for fall time threshold to +59.295% for the standard deviation in particle diameter. Hence, the effect of accounting for spatial dependence through weighting is non-negligible. As expected, the difference in total mass ejected is relatively low at +3.364% as it is known to be the parameter which is most independent of the rest and easiest to estimate (ref from Seb?).
The difference in fitted parameters manifests in the estimated tephra load surface as illustrated in Figure 8. The contours from the weighted calibration, denoting areas of similar load values, are less elongated towards the left than those from the unweighted calibration. In addition, there is generally lower loads at distal and medial regions (between the 30-100 kg/m² contours).

With the estimated variogram parameters (nugget, partial sill and range), we can generate 100 sample tephra load datasets from the same distribution as the 2014 Kelud data by adding simulated residual data to the Tephra2 fitted values from weighted inference. This means that we treat the Tephra2 fit as a mean surface of tephra load and use the Gaussian process with the fitted variogram to generate spatially correlated residuals. Examples of two simulated tephra load datasets are shown in Figure 12 of the Appendix.
We perform weighted and unweighted calibration of the Tephra2 model for each of these 100 simulated datasets to examine the differences in the inferred parameters. Figure 9 shows the box plots of the parameter estimates from the two calibration methods and Table 2 summarises the results in terms of percentage median bias (the difference between the median estimate and the true parameter value expressed in terms of percentage of the latter) as well as percentage median absolute deviation, a measure of spread. From these, we see that weighted calibration clearly improves inference of certain parameters. The median bias for beta parameter and median grain size decreases from 23.86% to 11.94% and from 8.97% to 6.56% respectively. In addition, their median absolute deviation decreases from 61.93% to 55.15% and from 41.95% to 35.31% respectively. Apart from the diffusion coefficient (for which accuracy and precision seems unaffected) and the standard deviation in particle diameter (for which estimates slight increase in absolute bias), weighting based on spatial conditional information either improves accuracy or precision for the remaining Tephra2 parameters. The median absolute deviation of total mass ejected decreases from 30.61% to 25.37% while that for fall time threshold decreases from 83.23% to 65.60%. On the other hand, maximum column height sees a reduction in median bias from 7.46% to 5.75%. Note that the procedure of estimating parameters from data simulated from the fitted model can be used to construct parametric bootstrap confidence intervals for the parameters. In the cases where the median absolute deviation is lower for estimates from the weighted calibration, the corresponding confidence intervals will be narrower than those derived from unweighted calibration.
The unequal effects of weighting on the different parameters is likely due to their interactions in the physical model. Figure 10 shows that there are some high correlations between parameters estimates from both unweighted and weighted calibration. In particular, standard deviation in particle diameter (std_dev) is highly correlated to median grain size (median_size) and total mass ejected (tot_mass) and the beta parameter (beta_param) is highly correlated with the maximum column height (max_column_ht).
Figure 11 shows the maximum absolute difference in weights between consecutive calibration iterations. Using a convergence threshold of 0.1 to determine if additional iterations are required, 90 out of the 100 datasets had weights which converged after one round of weighting. Only 6 out of the 100 datasets required two rounds of weighting for the weights to converge while 4 out of the 100 datasets required three rounds.

In the context of system biology, a similar problem is faced where limited experimental data is used to fit complex models and not all data points provide equal amounts of information on the system. To address this, Jeong & Qiu (2017) propose a calibration approach, where more informative data points were given higher weights. In their case, these were data points taken during the dynamic phase of the experiment rather than those taken during the steady state.
In their paper, Jeong & Qiu (2017) consider the weighted sum of squares as a cost function:

where $\theta$ represents the model parameters, $Y_{i}^{obs}$ the $i^{th}$ observation, $Y_{i}^{pred}$ the corresponding prediction and $w_{i}$ its allocated weight. Through a second-order Taylor expansion of $C$ centered around the optimal parameters and some simplifications, parameter uncertainty can be expressed in terms of the Fisher information matrix. In a similar way, the uncertainty of estimating one set of data points using another set of data points can be expressed in terms of their Fisher information matrices.
The iterative algorithm is as follows:

1. 1.

   Assign equal weight to all data points and use a random initial parameter setting.

2. 2.

   Optimise the parameters using an initial unweighted cost function.

3. 3.

   Compute Jacobian $J$ at the optimised parameter values.

4. 4.

   Compute the uncertainty of data point $i$ given all the other data points as:

   $$U(i|S_{-i})=\frac{1}{m}trace(I_{i}I_{S_{-i}}^{-1})$$

   (23)

   where $m$ is the number of parameters, $I_{i}=J_{i}^{T}J_{i}$ where $J_{i}$ is the $i^{th}$ row of $J$ and $I_{S_{-i}}$ is computed using all the other rows of $J$. The term for which we compute trace approximates the derivatives of the other data points with respect to $i$ (Jeong 2019, p.42).

5. 5.

   Compute the weights $w_{i}$ by normalising the uncertainty values so that they sum up to the total number of data points $n$.

6. 6.

   Use the weights in the next round of model calibration via the weighted cost function and use the previously estimated parameters as initial parameter settings.

7. 7.

   Iterate the calibration and and computing of weights until the weights converge.