Linear regression analysis of coastal data
Linear regression analysis of beach level data is demonstrated here using a set of beach profile measurements carried out at locations along the Lincolnshire coast (UK) by the National Rivers Authority (now the Environment Agency) and its predecessors between 1959 and 1991, as described in Sutherland et al. 2007[1]. A short introduction to linear regression theory is given in Appendix A.
Contents
Use of trend line for prediction
In order to better understand and predict beach level lowering in front of a seawall, beach levels were measured in front of seawalls along the Lincolnshire coast (UK). Figure 1 shows the bed level change measured in front of the seawall at Mablethorpe Convalescent Home.
Straight lines fitted to beach level time series give an indication of the rate of change of elevation and hence of erosion or accretion. The measured rates of change are often used to predict future beach levels by assuming that the best-fit rate from one period will be continued into the future. Similarly, long-term shoreline change rates can be determined using linear regression on cross-shore position versus time data.
Genz et al. (2007[2]) reviewed methods of fitting trend lines, including using end point rates, the average of rates, ordinary least squares (including variations such as jackknifing, weighted least squares and least absolute deviation, with and without weighting functions). Genz et al. recommended that weighted methods should be used if uncertainties are understood, but not otherwise. The ordinary least squares (OLS), jackknifing and least absolute deviation (LAD) methods were preferred (with weighting, if appropriate). If the uncertainties are unknown or not quantified then the LAD method is preferred.
The following question then arises: how useful is a best-fit linear trend as a predictor of future beach levels? In order to examine this, the thirty years of Lincolnshire data have been divided into sections: from 1960 to 1970, from 1970 to 1980, from 1980 to 1990 and from 1960 to 1990, for most of the stations. In each case a least-squares best-fit straight line was fitted to the data and the rates of change in elevation from the different periods are shown below:
- From 1960 to 1970 the rate of change was -17 mm/year;
- From 1970 to 1980 the rate of change was -63 mm/year;
- From 1980 to 1990 the rate of change was +47 mm/year;
- From 1960 to 1990 the rate of change was -25 mm/year.
The data above indicate that 10-year rates provide little predictive capability for estimating the change in elevation for the next 10 years, let alone for the planning horizon that might need to be considered for a coastal engineering scheme. The large differences between the decadal rates, including a reversal from erosion to accretion, illustrate the limited reliability of short-term trends for predicting longer-term beach evolution.
The prediction horizon is the lead time up to which the extrapolated trend outperforms the baseline prediction. Sutherland et al. (2007)[1] devised a method of determining the prediction horizon for an extrapolated trend using the Brier Skill Score, here defined as 1 minus the mean squared error of the prediction divided by the mean squared error of the baseline prediction (Sutherland et al., 2004[3]). The prediction horizon is then the lead time at which the score ceases to be positive. The baseline prediction was that future beach levels would be the same as the average of the measured levels used to define the trend. A 10 year trend was found to have a prediction horizon of only 4 years at Mablethorpe Convalescent Home (Fig. 2). Similar values have been found at other sites in Lincolnshire.
Uncertainty in extrapolated beach levels increases with the prediction horizon. A distinction should be made between confidence intervals for the estimated mean trend and prediction intervals for future beach levels, which additionally include variability around the trend. Both rely on the assumption that the statistical properties of the underlying process remain unchanged. They therefore do not account for changes in coastal processes that alter the rate of erosion or accretion. For beach level time series, errors are often correlated because of seasonal cycles, storm erosion and recovery, and multi-year fluctuations; the large differences between the decadal rates at Mablethorpe suggest such fluctuations. Confidence limits from ordinary least squares would then be too narrow. A fitted linear trend describes the observed change but does not identify the processes responsible for it; extrapolation is therefore particularly uncertain when the governing coastal conditions change.
For a discussion of beach lowering in front of seawalls see also the article Seawalls and revetments.
Conditions for application
A linear trend provides an estimate of the average rate of change over the observation period. Its use for prediction requires the additional assumption that this average rate remains approximately representative of future changes.
For ordinary least squares to provide the best linear unbiased estimates and for the conventional uncertainty formulas to apply, the regression errors are assumed to have zero mean, constant variance and no mutual correlation. A Gaussian distribution of errors is additionally assumed for conventional statistical inference based on the fitted coefficients.
The histogram in Fig. 3 suggests that the residual beach levels at Mablethorpe are approximately Gaussian-distributed. This does not imply that successive residuals are statistically independent. The important requirement of uncorrelated observation errors can be checked with an autocorrelation plot.
If observation errors are significantly correlated in time or space, ordinary least-squares uncertainty estimates are generally not appropriate and the correlation structure should be taken into account. Methods include generalized least squares and, for spatial data, Kriging.
Appendix A: Introduction to linear regression theory
Although regression analysis is explained in many textbooks, a short mathematical introduction is given below. The purpose of linear regression is to analyze the relationship between a measurable property [math]h[/math] (the so-called 'dependent variable' or 'target') and a number of [math]K-1[/math] measurable factors [math]y_k, \; k=2, .., K[/math] (the 'independent variables' or 'regressors'). Linear regression is based on the assumption that [math]h[/math] can be estimated through a linear relationship:
[math]h = \hat{h} +\hat{\epsilon} \; , \quad \hat{h} = a_1 + \sum_{k=2}^K a_k y_k = \sum_{k=1}^K a_k y_k \; , \qquad (A1)[/math]
where [math]\hat{h}[/math] is the estimator, [math]a_k, \; k=2, .., K[/math] are the regression coefficients, [math]a_1[/math] is the intercept and [math]y_1=1[/math]. The term [math]\hat{\epsilon}[/math] is the 'residual', i.e., the difference between the observed value [math]h[/math] and the estimator [math]\hat{h}[/math]. (Note: Linear regression means linear in the regression coefficients; the regressors [math]y_k[/math] can be nonlinear quantities, for example [math]y_3=y_2^2[/math].) Now assume that we have observations of [math]N[/math] different instances of the [math]K-1[/math] regressors, [math]y_{ki}, \; i=1, .., N[/math], and corresponding observations of the dependent variable, [math]h_i[/math]. If [math]N\gt K[/math], these observations can be used to estimate the values of the regression coefficients [math]a_k[/math] by looking for the best solution to the linear system
[math]h_i = \sum_{k=1}^K a_k y_{ki} + \hat{\epsilon}_i \; . \qquad (A2)[/math]
The deviations from the assumed regression relationship are represented by stochastic errors [math]\epsilon_i[/math], which account for unresolved variability, model approximation and measurement errors in the dependent variable. The errors are assumed to have zero mean, [math]E[\epsilon_i]=0[/math], and variance [math]\sigma_i^2 \equiv E[\epsilon_i^2][/math] (we use the notation [math]E[x][/math]= the mean value from a large number of trials of the stochastic variable [math]x[/math]). After fitting the regression coefficients, the corresponding deviations between observations and fitted values are called residuals, denoted by [math]\hat{\epsilon}_i[/math].
What is 'the best solution'? Different criteria can be used. Least-squares optimization is the most often used criterion, i.e., the solution for which half the sum of the squared residuals [math]\Phi[/math] is minimum,
[math]\Phi = \tfrac{1}{2} \sum_{i=1}^N \hat{\epsilon}_i^2 = \tfrac{1}{2} \sum_{i=1}^N \Big(h_i - \sum_{k=1}^K a_k y_{ki}\Big)^2 . \qquad (A3)[/math]
At minimum, [math]\Phi[/math] increases for any change in one of the coefficients [math]a_k[/math], which implies that the partial derivatives are zero:
[math]\Large\frac{\partial \Phi}{\partial a_k}\normalsize = 0 , \; k=1, ..,K . \qquad (A4)[/math]
This condition yields the set of [math]K[/math] linear equations
[math] - \sum_{i=1}^N y_{ki} h_i + \sum_{i=1}^N \sum_{k'=1}^K y_{ki} y_{k'i} a_{k'} = 0 , \qquad (A5)[/math]
from which the regression coefficients [math]a_k[/math] can be solved.
In the more compact matrix notation we may write: [math]H[/math] is the [math]N[/math]-vector with elements [math]h_i[/math], [math]\; A[/math] is the [math]K[/math]-vector with elements [math]a_k[/math] and [math]Y[/math] is the [math]N\times K[/math] matrix with elements [math]Y_{ik}=y_{ki}[/math] (row i contains the regressor values of observation i). The linear equations (A5) can then be written
[math]Y^T \, H = Y^T \, Y \, A . \qquad (A6)[/math]
The solution is found by inversion of the [math]K \times K[/math] matrix [math]Y^T \, Y[/math], giving
[math]A = (Y^T \, Y)^{-1} Y^T \, H . \qquad (A7)[/math]
This solution requires the regressors to be linearly independent, so that [math]Y^TY[/math] is invertible.
The estimated covariance matrix of the regression coefficients is [math]\text{Cov}(A) = s^2 (Y^T Y)^{-1}[/math], where [math]s^2 = \sum_{i=1}^N \hat{\epsilon}_i^2 / (N-K)[/math] is the estimated error variance and [math]\hat{\epsilon}_i[/math] are the residuals of the fit. The standard error of coefficient [math]a_k[/math] is the square root of the k-th diagonal element. This estimate is only valid if the residuals are uncorrelated and of equal variance; with positively correlated residuals it underestimates the true uncertainty.
Under the usual linear-model assumptions, least squares provides the best linear unbiased estimator when the observation errors [math]\epsilon_i[/math] have zero mean, equal variance [math]\sigma_i^2=\sigma^2[/math] and are uncorrelated. For conventional exact statistical inference it is also assumed that errors are Gaussian-distributed. This is often not the case in practice. A few other cases are discussed below.
Case 1. If positive data have approximately multiplicative errors, so that error magnitude scales with the measured (positive) value, a logarithmic transformation can often stabilize the variance. Linear regression can then be applied to the log-transformed data. Back-transformed predictions estimate the median rather than the mean.
Case 2. The error variance differs among observations ('heteroscedasticity'), [math]\sigma_i \ne \sigma_j , \, i \ne j[/math]. If estimates of the variances are known, one can use weighted least squares instead of least squares, i.e., replace (A3) with [math]\Phi = \tfrac{1}{2} \sum_{i=1}^N \big( \hat{\epsilon}_i / \sigma_i \big)^2[/math].
Case 3. Multicollinearity occurs when two or more regressors are strongly linearly related. Perfect linear dependence makes the regression coefficients indeterminate and requires removal or reformulation of redundant regressors. Strong but imperfect multicollinearity does not prevent fitting, but increases the uncertainty of individual regression coefficients.
Case 4. The errors of the [math]N[/math] observations are correlated, [math]c_{ij} \equiv E[\epsilon_i \, \epsilon_j] \ne 0[/math]. If the error covariance matrix [math]c_{ij} [/math] can be estimated, the generalized least squares method can be used. The solution (A7) is then replaced by
[math]A = (Y^T \, C^{-1} \, Y)^{-1} Y^T \, C^{-1} H , \qquad (A8)[/math]
where [math]C[/math] is the covariance matrix with elements [math]c_{ij}[/math]. For spatial prediction, methods such as Kriging can also be used to account for spatial correlation; see Data interpolation with Kriging.
Application to Mablethorpe example
Estimating the trend lines in the Mablethorpe example is the simplest case of Eq. (A1): one regressor, the time [math]t[/math], so that [math]K=2[/math], [math]y_2=t[/math] and [math]\hat{h} = a_1 + a_2 t[/math]. The slope [math]a_2[/math] is the rate of change of the beach level. Solving Eqs. (A5) for this case gives
[math]a_2 = \dfrac{\sum_{i=1}^N (t_i-\bar{t})(h_i-\bar{h})}{\sum_{i=1}^N (t_i-\bar{t})^2} , \quad a_1 = \bar{h} - a_2 \bar{t} , \qquad (A9)[/math]
where [math]\bar{t}[/math] and [math]\bar{h}[/math] are the mean values of [math]t_i[/math] and [math]h_i[/math]. The standard error of the rate [math]a_2[/math] follows from the covariance matrix:
[math]s_{a_2} = \dfrac{s}{\sqrt{\sum_{i=1}^N (t_i-\bar{t})^2}} , \quad s^2 = \dfrac{1}{N-2} \sum_{i=1}^N \hat{\epsilon}_i^2 . \qquad (A10)[/math]
The standard error decreases as the record becomes longer, because [math]\sum (t_i-\bar{t})^2[/math] grows rapidly with the length of the record. Like the general covariance estimate, Eq. (A10) assumes uncorrelated errors of equal variance; for positively correlated residuals it underestimates the uncertainty of the rate.
Related articles
- Data analysis techniques for the coastal zone
- Analysis of coastal processes with Empirical Orthogonal Functions
- Wavelet analysis of coastal processes
- Artificial Neural Networks and coastal applications
- Data interpolation with Kriging
- Random Forest Regression
- Support Vector Regression
References
- ↑ 1.0 1.1 1.2 Sutherland, J., Brampton, A.H., Obhrai, C., Motyka, G.M., Vun, P.-L. and Dunn, S.L. 2007. Understanding the lowering of beaches in front of coastal defence structures, Stage 2. Defra/EA Joint Flood and Coastal Erosion Risk Management R&D programme Technical Report FD1927/TR
- ↑ Genz, A.S., Fletcher, C.H., Dunn, R.A., Frazer, L.N. and Rooney, J.J. 2007. The predictive accuracy of shoreline change rate methods and alongshore beach variation on Maui, Hawaii. Journal of Coastal Research 23(1): 87 – 105
- ↑ Sutherland, J., Peet, A.H. and Soulsby, R.L. 2004. Evaluating the performance of morphological models. Coastal Engineering 51, pp. 917-939.
|
Please note that others may also have edited the contents of this article.
|
