slr_slope: Calculate Long Term Estimate of SLR
In ccb60/SLRSIM: Tools to Analyse Sea Level Rise

Description Usage Arguments Details Value See Also Examples

Given an series of (typically autocorrelated) observations of sea level or mean sea level, and associated time coordinates (usually Dates), return a named vector containing a summary of a generalized least squares linear trend, optionally scaled to annual values. This function mimics the calculation of long term estimates of sea level rise and associated uncertainties conducted by NOAA. Where we have checked, results of this function, when based on NOAA's seasonally adjusted mean trend data (See Details), agree with NOAA's trend estimates within rounding error.

slr_slope(
  .data,
  .sl,
  .dt,
  .ci = 0.95,
  t_fit = FALSE,
  .mode = c("year", "unscaled")
)

`.data`	Source data frame for data. Use NULL if no data frame is used and all data is passed from the enclosing environment.
`.sl`	The data variable (usually found in the data frame) showing sea level or mean sea level. Must be a named variable, not an expression.
`.dt`	Data variable containing corresponding midpoint dates for the period of averaging used to calculate .sl. Must be a named variable, not an expression. Midpoint dates for a given month can be approximated with `as.Date(paste(year, month, 15, sep = '-')`.
`.ci`	P value for two sided confidence intervals for the estimated slope, based on a (possibly naive) normal approximation.
`t_fit`	Should the underlying generalized linear model's correlation structure be fit based on the time coordinate (.dt), or only on the sequence of observations in the data? Setting this to TRUE is safer if you are uncertain of the sequence of observations in the source data, or significant missing values in the data, but it significantly slows model fitting. Results tend to be very similar for complete or near complete data.
`.mode`	one of c('year', 'unscaled') indicating whether the results should be scaled to years (the default) or left unscaled. Other sea level rate functions accept other `.mode` values. Scaling requires the time coordinate to inherit from R `Date` or `POSIXt` objects. Rates based on `Dates` are converted to annual values by multiplying by 365.2422 days / year. Rates based on `POSIXt` objects are converted to annual values by multiplying slope estimates by multiplying by 31556926 seconds / year (with a warning about rounding error).

Sea level data tends to be highly autocorrelated over multiple time scales. This function uses generalized least squares (GLS) to fit a linear model with a correlation structure defined by fitting an AR1 process to model residuals.

Slope estimates from the AR1 models tend (in this context) to be indistinguishable from estimates generated by ordinary least squares, but the uncertainty of the estimates is substantially greater. As observations are not independent, ordinary least squares (which assumes independence) will overstate model precision.

The AR1 process, an "autoregressive process of order 1", is a simple assumption about correlation structures that reduces the risk of overstating the precision of trend estimates. In an AR1 model, each observation is regressed against the prior observation. Informally, An AR1 process generates a correlation structure based on the lag 1 autocorrelation. Modeled autocorelations between observations are thus proportional to the lag 1 autocorrelation raised to the power of the number of timesteps between the two observations.

The AR1 structure, however, may not be ideal for all tidal time series. Actual tidal data series tend to have more complex autocorrelation structure. Where regular seasonal patterns are evident, explicit modeling is probably called for. Thus the function is principally intended for use with data averaged over monthly or longer periods.

This function is a convenience wrapper around nlme::gls(..., correlation = corAR1()). For more control, or to examine model diagnostics, the user should fit the generalized least squares model directly.

The function does not fit periodic terms. Thus it is generally not suitable for use with "raw" sea level observations, especially those collected at high resolution, like NOAA hourly data. Tidal data contains many periodic components (~ 24:50 tidal period, spring tide / neap tide cycles, annual cycles, etc.). Those periodic structures will affect estimates of autocorrelation, ultimately making estimates of model uncertainty less reliable. If working from high frequency data, short-term structure should be removed before using this function, typically by calculating average sea level over convenient and meaningful time periods, like months or years.

The function's original use case was based on secondary analysis of NOAA "mean trend" data. NOAA calculates long-term sea level rise estimates based on a modified monthly water level data series. The source data is not the same as what one gets if one downloads monthly mean water level through the standard NOAA API. NOAA's "mean trend" data contains "monthly mean sea level without the regular seasonal fluctuations due to coastal ocean temperatures, salinities, winds, atmospheric pressures, and ocean currents. The seasonally detrended data is available from NOAA for manual download from water level station web pages, but apparently not through the API. The prov_meantrend data included in SLRSIM is an example, for Providence, Rhode Island.

SLRSIM provides a a function, get_slr_trend(), to access NOAA's long term trend analysis result directly, including error terms and period of record.

If t_fit = FALSE, the function fits an autocorrelation structure to the sequential observations, ordered by .dt, without regard to how much time has passed between subsequent observations. This approach works well with evenly spaced dates and complete or nearly complete records, but is less successful if there are many periods with missing observations.

If t_fit = TRUE, the function fits an autocorrelation structure based on the time coordinate (.dt). In this case, the time coordinate must be integer valued. R's Date class is an integer under the hood. It is easy to build Dates, either by converting POSIX times with as.Date() or by building them up from month, day, and year information with as.Date(paste(year, month, 15, sep = '-').

By default, analysis is conducted by years (.mode == 'year'). The slope estimate and standard error are scaled from days (based on R Dates) or seconds (based on POSXIt objects) to (approximate) annual values by multiplying by 365.242 2 (days / Dates) or 31556926 (seconds / POSIXt).

A named vector, with the following components and attributes.

Vector Components

Estimate: The trend. Units depend on value of the .mode argument. For .mode == 'year, units are per year, otherwise per unit of the .dt variable.
Std_Err: Standard Error of the trend. Units as above.
P_Val: P value of the trend.
Lower_CI: Lower confidence interval for the estimate.
Upper_CI: Upper confidence interval for teh estimate.

Attributes

CI_P: Nominal p value used to generate the confidence intervals.
span: How many calander years are represented in the data?
start: First year included in the data.
stop: last year included in the data.
cor_struct: Were parameters fit based on order or time

Other sea level rate functions: get_sl_trend(), slr_change_comp(), slr_change()

prov_meantrend$MSL_mm <- prov_meantrend$MSL * 1000 (convert to mm)

# Basic Usage
slr_slope(prov_meantrend, MSL_mm, MidDate)

# if not passed dates, the result is unscaled
 slr_slope(prov_meantrend, MS_mm, unclass(MidDate))

# One can also fit the model to the time coordinate explicitly
slr_slope(prov_meantrend, MSL_mm, MidDate, t_fit = TRUE)