How should surface elevation table data be analyzed? A comparison of several commonly used analysis methods and one newly proposed approach

Abstract

The use of surface elevation table (SET) instruments to monitor elevation changes at low elevation coastal locations has steadily increased in recent years. A primary focus in the analysis of SET data is the estimation of the overall rate of elevation change, and numerous approaches have been used for this purpose. In this work, we compare and contrast several methods used for estimating the true rate of elevation change at SET station locations, including a novel approach proposed in this work that incorporates spatial dependence. We also discuss theoretical properties of one class of estimators, and undertake a comprehensive simulation study. Additionally, we present two case studies where we illustrate these differences using real SET data. All methods considered here tend to produce similar point estimates, but some confidence interval procedures can generate intervals with empirical coverage rates lower than specified. However, the best analysis approach is likely dependent upon selecting the method that best coincides with the true underlying process. Thus, we do not uniformly recommend one approach for all situations. Instead, we suggest carefully weighing potential advantages and disadvantages of each method before conducting analysis, while keeping in mind the ways in which modeling assumptions may impact this decision.

Appendices

Appendix A: Properties of the least-squares based slope estimators

When comparing the relative strengths of different estimators of the same parameter, statisticians may consider the mean squared error (MSE) criterion. For arbitrary parameter $\theta$ with estimator ${\hat{\theta }}$, its MSE is given by

$$\begin{aligned} MSE({\hat{\theta }}) = E\left( ({\hat{\theta }} - \theta )^2\right) . \end{aligned}$$

Interestingly, MSE may be expressed as the sum of two components: the variance of the estimator and the square of the estimator’s bias. We note that the bias of an estimator is defined to be

$$\begin{aligned} \text {Bias}({\hat{\theta }}) = E({\hat{\theta }}) - \theta . \end{aligned}$$

An estimator is said to be unbiased if its bias is exactly equal to zero. When comparing unbiased estimators, comparison based on the MSE criterion can be based solely on an estimator’s variance, $Var({\hat{\theta }})$. Since $Var({\hat{\theta }})$ may involve unknown parameters, we use the sample analog ${\widehat{Var}}({\hat{\theta }})$ in practice, noting that its corresponding SE is denoted

$$\begin{aligned} SE({\hat{\theta }}) = \sqrt{{\widehat{Var}}({\hat{\theta }})}. \end{aligned}$$

In this appendix, we discuss these properties for the estimators ${\hat{\beta }}^{ATR}$ and ${\hat{\beta }}^{RTA}$. We establish that they are both unbiased, and provide an expression for the variance for each estimator. As these variances are functions of unknown parameters, we also propose estimated variances (and standard errors) which may be calculated based on sample data.

1.1 Properties of the ATR slope estimator

Assuming that no data are missing and $E(y_{t,a,p}) = \alpha + \beta t$ holds, as in the model given in Equation (4), ${\hat{\beta }}^{ATR}$ will be an unbiased estimator of $\beta$. This can be shown by

$$\begin{aligned} E({\hat{\beta }}^{ATR})= & {} E \left( \frac{1}{AP S_{tt}} \sum _{t \in {\mathcal {T}}} \sum _{a=1}^A \sum _{p=1}^P (t - {\bar{t}}) y_{t,a,p} \right) \\= & {} \frac{1}{AP S_{tt}} \sum _{t \in {\mathcal {T}}} \sum _{a=1}^A \sum _{p=1}^P (t - {\bar{t}}) E(y_{t,a,p}) \\= & {} \frac{1}{AP S_{tt}} \sum _{t \in {\mathcal {T}}} \sum _{a=1}^A \sum _{p=1}^P (t - {\bar{t}}) (\alpha + \beta t) \\= & {} \frac{1}{S_{tt}} \sum _{t \in {\mathcal {T}}} (t - {\bar{t}}) (\alpha + \beta t) \\= & {} \frac{1}{S_{tt}} \left( \alpha \sum _{t \in {\mathcal {T}}} (t - {\bar{t}}) + \beta \sum _{t \in {\mathcal {T}}} t (t - {\bar{t}}) \right) \\= & {} \frac{\beta }{S_{tt}} \sum _{t \in {\mathcal {T}}} t(t - {\bar{t}}) \\= & {} \beta . \end{aligned}$$

Also assuming that the pin elevations are independent and have the same error variance, we can derive an expression for the variance of ${\hat{\beta }}^{ATR}$ via

$$\begin{aligned} Var({\hat{\beta }}^{ATR})= & {} Var \left( \frac{1}{AP S_{tt}} \sum _{t \in {\mathcal {T}}} \sum _{a=1}^A \sum _{p=1}^P (t - {\bar{t}}) y_{t,a,p} \right) \\= & {} \left( \frac{1}{AP S_{tt}}\right) ^2 \sum _{t \in {\mathcal {T}}} \sum _{a=1}^A \sum _{p=1}^P (t - {\bar{t}})^2 Var(y_{t,a,p}) \\= & {} \left( \frac{1}{AP S_{tt}}\right) ^2 \sum _{t \in {\mathcal {T}}} \sum _{a=1}^A \sum _{p=1}^P (t - {\bar{t}})^2 \sigma _{a,p}^2 \\= & {} \left( \frac{1}{AP S_{tt}}\right) ^2 \sum _{t \in {\mathcal {T}}} \sum _{a=1}^A \sum _{p=1}^P (t - {\bar{t}})^2 \sigma ^2 \\= & {} \frac{\sigma ^2}{AP S_{tt}^2} \sum _{t \in {\mathcal {T}}} (t - {\bar{t}})^2 \\= & {} \frac{\sigma ^2}{AP S_{tt}}. \end{aligned}$$

Based on the pseudo observations $\{(t,y_t^*)\}_{t \in {\mathcal {T}}}$ and the model in (4), define $MSE^* = (N_t - 2)^{-1} \sum _{t \in {\mathcal {T}}} (y_t^* - {\hat{y}}_t^*)^2$, where ${\hat{y}}_t^* = {\hat{\alpha }}^{ATR} + {\hat{\beta }}^{ATR} t$. As is the case in SLR, $MSE^*$ will be an unbiased estimator of $\sigma ^{*2} = \sigma ^2/(AP)$; therefore, a plug-in estimator of $Var({\hat{\beta }}^{ATR})$ is

$$\begin{aligned} {\widehat{Var}}({\hat{\beta }}^{ATR}) = \frac{MSE^*}{S_{tt}}. \end{aligned}$$

(8)

This estimated variance can be used to constuct hypothesis tests and confidence intervals for ${\hat{\beta }}^{ATR}$ in the usual manner. We also point out that we define the standard error for ${\hat{\beta }}^{ATR}$, denoted $SE({\hat{\beta }}^{ATR})$, to be the square root of the estimated variance in Equation (8).

1.2 Properties of the RTA slope estimator

We can show similar properties for the RTA slope estimator. Assuming that $E({\hat{\beta }}_{a,p}) = \beta ~ \forall ~ a \in {\mathcal {A}}, p \in {\mathcal {P}}$, which will occur if $\beta _{a,p} = \beta ~ \forall ~ a \in {\mathcal {A}}, p \in {\mathcal {P}}$,

$$\begin{aligned} E({\hat{\beta }}^{RTA})= & {} E\left( \frac{1}{AP} \sum _{a=1}^A \sum _{p=1}^P {\hat{\beta }}_{a,p} \right) \\= & {} \frac{1}{AP} \sum _{a=1}^A \sum _{p=1}^P E({\hat{\beta }}_{a,p}) \\= & {} \frac{1}{AP} \sum _{a=1}^A \sum _{p=1}^P \beta \\= & {} \beta . \end{aligned}$$

This verifies that ${\hat{\beta }}^{RTA}$ is an unbiased estimator of $\beta$ if no data are missing. Further assuming independence of the pin elevation measurements and common error variance for the underlying true regression models, an expression for the variance of ${\hat{\beta }}^{RTA}$ is given by

$$\begin{aligned} Var({\hat{\beta }}^{RTA})= & {} Var \left( \frac{1}{AP} \sum _{a=1}^A \sum _{p=1}^P {\hat{\beta }}_{a,p} \right) \\= & {} \left( \frac{1}{AP}\right) ^2 \sum _{a=1}^A \sum _{p=1}^P Var({\hat{\beta }}_{a,p}) \\= & {} \left( \frac{1}{AP}\right) ^2 \sum _{a=1}^A \sum _{p=1}^P \left( \frac{\sigma _{a,p}^2}{S_{tt}} \right) \\= & {} \left( \frac{1}{AP}\right) ^2 \sum _{a=1}^A \sum _{p=1}^P \left( \frac{\sigma ^2}{S_{tt}} \right) \\= & {} \frac{\sigma ^2}{AP S_{tt}}. \end{aligned}$$

Interestingly, we observe that $Var({\hat{\beta }}^{ATR})$ is equal to $Var({\hat{\beta }}^{RTA})$. Furthermore, we consider two approaches for estimating $Var({\hat{\beta }}^{RTA})$. In the first approach, analysis is conducted by assuming that $\{ {\hat{\beta }}_{a,p} \}_{a \in {\mathcal {A}}, p \in {\mathcal {P}}}$ constitutes a random sample to be used for the purpose of estimating the true underlying mean rate of change, $\beta$. To this end, we estimate $Var({\hat{\beta }}(\varvec{s}))$ using the sample variance of the slope estimates from each of the AP pins

$$\begin{aligned} {\widehat{Var}}_1({\hat{\beta }}^{RTA}) = \frac{1}{AP - 1} \sum _{a=1}^A \sum _{p=1}^P \left( {\hat{\beta }}_{a,p} - {\hat{\beta }}^{RTA} \right) ^2. \end{aligned}$$

(9)

Inference is performed using Student’s t based procedures with $AP - 1$ degrees of freedom. This strategy is analogous to estimating the mean of a Normal population with unknown variance based on a simple random sample.

A second approach for quantifying uncertainty is to observe that ${\hat{\sigma }}^2 = \frac{1}{AP} \sum _{a=1}^A \sum _{p=1}^P MSE_{a,p}$ is an unbiased estimator of $\sigma ^2$, as

$$\begin{aligned} E({\hat{\sigma }}^2)= & {} E \left( \frac{1}{AP} \sum _{a=1}^A \sum _{p=1}^P MSE_{a,p} \right) \\= & {} \frac{1}{AP} \sum _{a=1}^A \sum _{p=1}^P E(MSE_{a,p}) \\= & {} \frac{1}{AP} \sum _{a=1}^A \sum _{p=1}^P \sigma _{a,p}^2 \\= & {} \frac{1}{AP} \sum _{a=1}^A \sum _{p=1}^P \sigma ^2 \\= & {} \sigma ^2. \end{aligned}$$

This leads to the plug-in estimator of $Var({\hat{\beta }}^{RTA})$,

$$\begin{aligned} {\widehat{Var}}_2({\hat{\beta }}^{RTA}) = \frac{\sum _{a=1}^A \sum _{p=1}^P MSE_{a,p}}{(AP)^2 S_{tt}}. \end{aligned}$$

(10)

Appendix B: Highlighting differences among least-squares based methods

Differences between the least squares based methods are tied to the ways in which they estimate uncertainty; therefore, it is useful to highlight these dissimilarities through four carefully constructed data sets. Figure 11 gives scatterplots of these four scenarios that are designed to compare and contrast the least squares based methods’ abilities to estimate variability in the overall slope estimator.

For each data set, two pins are measured at three time points. The measurements and regression line associated with the first pin are given in red, while the measurements and regression line associated with pin two are given in blue. The ATR pseudo observations are represented with black stars, and the corresponding regression line is given by a black dotted line.

Applying the estimators defined in Sect. 2.2 to the four data sets plotted in Fig. 11 (given in Table 5) highlights key differences among these methods. In this table, rows one and two give the estimated slope for pin one and two (respectively). Rows three and four give the overall rate of change estimate based on the RTA and ATR methods (respectively). Rows five and six give MSE in the SLR models for pin one and two (respectively). Row seven gives the MSE in the SLR model for the psuedo observations. Rows eight, nine, and ten give each of the standard errors for the overall rate of change estimate.

Interestingly, none of the methods for constructing standard errors produces values that are consistently higher or lower than the others. In the data set plotted in top left panel, we note that $MSE_1 = MSE_2 = MSE^* = 0$, implying that ${\widehat{Var}}_2({\hat{\beta }}^{RTA}) = {\widehat{Var}}({\hat{\beta }}^{ATR}) = 0$. However, the slope estimates for the two pins differ, suggesting that ${\widehat{Var}}_1({\hat{\beta }}^{RTA}) > 0$. In contrast, for the data set plotted in the top right panel, $MSE_1 = MSE_2 > 0$, but $MSE^* = 0$. This indicates that ${\widehat{Var}}_2({\hat{\beta }}^{RTA}) > {\widehat{Var}}({\hat{\beta }}^{ATR})$.

In the bottom left panel of Fig. 11, $MSE_1 = 0$; however, $MSE_2> MSE^* > 0$. This implies that ${\widehat{Var}}_2({\hat{\beta }}^{RTA}) > {\widehat{Var}}({\hat{\beta }}^{ATR})$. In the bottom right panel of Fig. 11, $MSE_1 = MSE_2 = MSE^* > 0$, resulting in ${\widehat{Var}}_2({\hat{\beta }}^{RTA}) = {\widehat{Var}}({\hat{\beta }}^{ATR}) > 0$. However, ${\widehat{Var}}_1({\hat{\beta }}^{RTA}) = 0$ since both slope estimates are identical.

Fig. 11

For illustrative purposes, we present four relatively simple scenarios that compare and contrast the ATR and RTA methods’ abilities to quantify uncertainty in the overall slope estimator. In all four data sets, two pins are measured at three time points. The measurements and regression line associated with pin one (two) are given in red (blue). The ATR pseudo observations are denoted with black stars and the corresponding regression line is given by a black dotted line

Table 5 For the four data sets plotted in Fig. 11, we present several statistics in order to highlight key differences among methods

Appendix C: Choice of data generating process and parameter settings in the simulation study

In this appendix, we provide the reader with intuition regarding the choice of data generating process and parameter settings.

The data generating process described in Eq. (7) is a complicated spatio-temporal process; however, it is useful to think about it as an extension of the model

$$\begin{aligned} Y_t(\varvec{s}) = \alpha + \beta t + \varepsilon _{t}(\varvec{s}), \end{aligned}$$

where $\varepsilon _{t}(\varvec{s}) \overset{iid}{\sim } N(0,\sigma ^2)$. This model with iid errors can simply be thought of as a SLR model, where the intercept at all locations is $\alpha$, and the rate of change (slope) at all locations is $\beta$. If we believe that initial elevations at all locations are not identical, we can add a term to account for this additional source of variability. After including the spatially correlated random effect $\theta _{\alpha }$, we obtain the model

$$\begin{aligned} Y_t(\varvec{s}) = \underbrace{\alpha + \theta _{\alpha }(\varvec{s})}_{Intercept} + \beta t + \varepsilon _{t}(\varvec{s}). \end{aligned}$$

Realizations of $\theta _{\alpha }$ will tend to be more similar for locations that are closer together; therefore, it can be thought of as a random surface that exhibits a degree of spatial smoothness. In the context of SET analysis, this could be interpreted as attempting to model the undulating surface in the immediate vicinity of the SET station at the initial time of measurement.

The models discussed to this point in this appendix all include iid errors, which may not be reasonable for actual SET station data. More specifically, errors for observations that are closer together in space and time may tend to show a higher degree of similarity. This can be modeled by including the spatio-temporal random process W via

$$\begin{aligned} Y_t(\varvec{s}) = \underbrace{\alpha + \theta _{\alpha }(\varvec{s})}_{intercept} + \beta t + \underbrace{W_t(\varvec{s}) + \varepsilon _{t}(\varvec{s})}_{error}. \end{aligned}$$

Under this model specification, the error term is composed of the sum of two random effects, W and $\varepsilon$. The process W allows for these errors to exhibit spatio-temporal dependence, and $\varepsilon$ is the iid error that is typical in linear regression. Including this type of iid error in a spatial or spatio-temporal process allows for what is sometimes referred to as a nugget effect.

We further assume that W is separable through space and time; that is, the spatial dependence structure does not change through time. One way to model separable spatio-temporal dependence is via the product between a spatial covariance function and a temporal covariance function. We employ this approach, and assume that both covariance functions are exponential with range parameters $\rho _W^{sp} = 0.032$ and $\rho _W^{tmp} = 8.9$. Using the traditional rule of thumb described in Banerjee et al. (2014), these parameter values imply effective ranges (the separation in space or time after which dependence is negligible) of $3/0.032 \approx 93.75$cm and $3/8.9 \approx 0.34$ years (respectively). These values were chosen after performing exploratory data analysis of several real SET data sets, but do not necessarily reflect any one SET station perfectly.