knitr::opts_chunk$set( collapse = TRUE, fig.ext = 'png', comment = "#>" ) library(waterlevel) library(earthtide) library(data.table) library(plotly) library(rlang) library(purrr) library(dplyr) library(recipes)
To test the efficiency of different methods we will create a response function, apply it to some barometric pressure data to generate a synthetic water level response.
This is the function we will try to recreate! It is the combination of three exponentials of different lengths:
60 second (negative)
600 second (negative)
43201 second (positive)
exp1 <- -exp(-seq(0.01, 8, length.out = 60)) exp2 <- -exp(-seq(0.01, 8, length.out = 43201)) exp3 <- exp(-seq(0.01, 8, length.out = 600)) exp1 <- exp1 / sum(exp1) exp2 <- exp2 / sum(exp2) exp3 <- exp3 / -sum(exp3) exp1 <- c(exp1, rep(0.0, length(exp2) - length(exp1))) exp3 <- c(exp3, rep(0.0, length(exp2) - length(exp3))) kern <- rev(exp1 + exp2 + exp3) resp_fun <- cumsum(rev(kern))*0.5
par(mar = c(5,5,0.5,0.5)) lag_resp_fun <- 0:43200 lag_resp_fun[1] <- 0.1 # for log axis plot(resp_fun~lag_resp_fun, type = 'l', las = 1, log = 'x', col = '#96BC56', lwd = 3, xlab = 'Time lag (s)', ylab = 'Cumulative response')
We can get a barometric response with the synthetic function.
Convolution of barometric pressure data and response function. Water levels have a dampened appearance and are lagged slightly in relation to the barometric response. This is commonly observed at field sites. There is no noise in this example so we hope to perfectly recover the response function.
dat <- synthetic(sd_noise = 0e-2, sd_noise_trend = 0e-12, n = 4 * 86400, linear_trend = 0e-12, seed = 1, baro_kernel = kern)[, list(datetime, baro, wl)] dat <- dat[1:(2.5 * 86400)] dat[, baro := baro - mean(baro, na.rm = TRUE)] dat[, wl := wl - mean(wl, na.rm = TRUE)]
par(mar = c(5,5,0.5,0.5)) plot(baro~datetime, dat, type = 'l', las = 1, lwd = 2, col = '#BC5696', xlab = '', ylab = 'pressure (dbar)') points(wl~datetime, dat, type = 'l', lwd = 2, col = '#5696BC')
Regular lag methods can generate large regression matrices, and therefore, we will subset our data prior to analysis. We start with one second data, but will subset to 60 seconds. The matrix size will be
dt <- 60 dat_sub <- dat[seq(1, nrow(dat), dt)] ba_lags <- seq(0, (43200+dt) / dt, 1) rec <- recipe(wl~., dat_sub) %>% step_lag_matrix(baro, lag = ba_lags) %>% prep() %>% portion() dim(rec)
What if we only had data every 30 minutes?
fit <- lm(outcome~lag_matrix, rec, x = FALSE, y = FALSE, tol = 1e-50, na.action = na.exclude) summarize_lm(fit)[, -c('adj_r_squared'), with = FALSE] resp <- response_from_fit(fit) resp$x <- resp$x * dt resp$x[1] <- 0.1 par(mar = c(5, 5, 0.5, 0.5)) plot(resp_fun~lag_resp_fun, type = 'l', las = 1, log = 'x', col = '#96BC56', lwd = 3, xlab = 'Time lag (s)', ylab = 'Cumulative response') points(value~x, resp, type = 'l', lwd = 2, pch = 20, col = '#5696BC')
We won't subset in this case. We will use fewer regressors but have many more observations. These results are even better.
ba_lags <- log_lags(251, max_time_lag = 43201) rec <- recipe(wl~., dat) %>% step_lag_matrix(baro, lag = ba_lags) %>% prep() %>% portion() dim(rec)
fit <- lm(outcome~lag_matrix, rec, x = FALSE, y = FALSE, tol = 1e-50, na.action = na.exclude) summarize_lm(fit)[, -c('adj_r_squared'), with = FALSE] resp <- response_from_fit(fit) resp$x[1] <- 0.1 par(mar = c(5, 5, 0.5, 0.5)) plot(resp_fun~lag_resp_fun, type = 'l', las = 1, log = 'x', col = '#96BC56', lwd = 3, xlab = 'Time lag (s)', ylab = 'Cumulative response') points(value~x, resp, type = 'l', lwd = 2, pch = 20, col = '#5696BC')
We won't subset in this case and we will only use 15 regressors.
ba_lags <- log_lags(18, max_time_lag = 43201) rec <- recipe(wl~., dat) %>% step_distributed_lag(baro, knots = ba_lags) %>% prep() %>% portion() dim(rec)
fit <- lm(outcome~distributed_lag, rec, x = FALSE, y = FALSE, tol = 1e-50, na.action = na.exclude) summarize_lm(fit)[, -c('adj_r_squared'), with = FALSE] resp <- response_from_fit(fit) resp$x[1] <- 0.1 par(mar = c(5, 5, 0.5, 0.5)) plot(resp_fun~lag_resp_fun, type = 'l', las = 1, log = 'x', col = '#96BC56', lwd = 3, xlab = 'Time lag (s)', ylab = 'Cumulative response') points(value~x, resp, type = 'l', lwd = 2, pch = 20, col = '#5696BC')
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.