Nothing
#' Test hydrograph changes
#'
#' Use this function to test interannual changes or hydrograph separation variables returned by [grwat::gr_summarize()]. Pettitt test is used to detect the change year — i.e. the year which divides the time series into the statistically most differing samples. Student (Welch) and Fisher tests are used to estimate the significance of mean and variance differences of these samples. Theil-Sen test calculates the trend slope value. Mann-Kendall test is performed to reveal the significance of the trend.
#'
#' Number of observations formally required for various tests: Pettitt > 0, Mann-Kendall > 2, Theil-Sen > 1, Student > 1, Fisher > 1.
#'
#' @param df `data.frame` as produced by [grwat::gr_summarize()] function.
#' @param ... Names of the tested variables (quoted).
#' @param year Integer value of year used to divide series in two samples compared by Student and Fisher tests. Defaults to `NULL` which means that the year is calculated automatically by Pettitt test.
#' @param exclude Integer vector of years to be excluded from tests.
#'
#' @return `list` of testing results with following elements:
#'
#' | __Element__ | __Description__ |
#' | ------------ | --------------- |
#' | `ptt` | Pettitt tests for change year |
#' | `mkt` | Mann-Kendall test for trend significance |
#' | `tst` | Theil-Sen test for slope estimation |
#' | `ts_fit` | Theil-Sen linear model fit |
#' | `tt` | Student (Welch) test for significance of mean differences between two periods |
#' | `ft` | Fisher test for significance of variance differences between two periods |
#' | `year` | Integer value of year used to divide series in two samples compared by Student and Fisher tests |
#' | `maxval` | Maximum value for the variable along the full time series |
#' | `fixed_year` | Boolean `TRUE` or `FALSE` value indicating if the year was fixed |
#' | `pvalues` | p-values of all tests summarized as a single table for all variables |
#'
#' @export
#'
#' @example inst/examples/gr_test_vars.R
#'
gr_test_vars <- function(df, ..., year = NULL, exclude = NULL){
fields = rlang::exprs(...) %>% as.character()
if(length(fields) == 0)
fields = params_out %>%
dplyr::filter(.data$Order != 0) %>%
dplyr::arrange(.data$Order) %>%
dplyr::pull(.data$Name)
prms = params_out %>%
dplyr::filter(.data$Name %in% fields) %>%
dplyr::slice(match(fields, .data$Name))
fixed = !is.null(year)
desc = switch(grenv$loc,
'RU' = prms$Desc,
'UA' = prms$Descua,
'EN' = prms$Descen)
nn = nrow(prms)
ch_year = setNames(rep(NA, nn), prms$Name) # change years
mean1 = setNames(vector(mode = 'list', length = nn), prms$Name) # means for first period
mean2 = setNames(vector(mode = 'list', length = nn), prms$Name) # means for second period
mratio = setNames(rep(NA, nn), prms$Name) # variance for first period
sd1 = setNames(rep(NA, nn), prms$Name) # variance for first period
sd2 = setNames(rep(NA, nn), prms$Name) # variance for first period
maxval = setNames(vector(mode = 'list', length = nn), prms$Name) # maximum values
ptt = setNames(vector(mode = 'list', length = nn), prms$Name) # Pettitt test
mkt = setNames(vector(mode = 'list', length = nn), prms$Name) # Mann-Kendall test
tst = setNames(vector(mode = 'list', length = nn), prms$Name) # Theil-Sen slope estimation
ts_fit = setNames(vector(mode = 'list', length = nn), prms$Name) # Theil-Sen regression
tt = setNames(vector(mode = 'list', length = nn), prms$Name) # Student t test for periods
ft = setNames(vector(mode = 'list', length = nn), prms$Name) # Fisher F test for periods
df = df %>%
dplyr::filter(!(.data$Year1 %in% exclude)) %>%
dplyr::mutate_if(params_out$Winter == 1,
replace_year)
cli::cli_progress_bar('Testing variables', total = nn)
for (i in 1:nn) {
cli::cli_progress_update()
# PETTITT TEST FOR CHANGE DETECTION
vl = dplyr::pull(df, prms$Name[i])
uvals = unique(vl)
if (length(uvals[!is.na(uvals)]) < 3)
next
isdate = FALSE
if(prms$Unitsen[i] %in% c('Date', 'Month')) {
isdate = TRUE
lubridate::year(vl) = 2000
if (prms$Winter[i] == 1) {
vl = replace_year(vl)
}
}
vl_cmp = !is.na(vl)
vl_cmp_sum = cumsum(vl_cmp)
ptt[[i]] = trend::pettitt.test(vl[vl_cmp])
nyear = match(ptt[[i]]$estimate, vl_cmp_sum)
ch_year[i] = ifelse(is.numeric(year),
year,
as.vector(as.matrix(df[nyear, "Year1"]))[1])
maxval[[i]] = max(vl, na.rm = TRUE)
# MANN-KENDALL TEST FOR TREND SIGNIFICANCE
if(length(vl_cmp) > 2) { # mk.test requires at least 3 observations
if(isdate){
mkt[[i]] = trend::mk.test(vl[vl_cmp] %>% as.integer())
} else {
mkt[[i]] = trend::mk.test(vl[vl_cmp])
}
}
# THEIL-SEN SLOPE ESTIMATION
df.theil = df %>%
dplyr::select('Year1', prms$Name[i])
# values = df.theil[[prms$Name[i]]] %>%
# as.matrix() %>%
# as.vector()
if(isdate){
# values = values %>%
# as.Date() %>%
# as.integer()
df.theil[2] = as.integer(vl)
}
frml = substitute(y ~ x,
list(y = as.name(prms$Name[i]),
x = as.name('Year1')))
if (nrow(df.theil) > 1) { # slope testing requires at least two observations
fltr = ! (is.infinite(df.theil[[2]]) | is.nan(df.theil[[2]]) | is.na(df.theil[[2]]))
ts_fit[[i]]= mblm::mblm(eval(frml), data = df.theil[fltr, ], repeated = FALSE)
tst[[i]] = trend::sens.slope(df.theil[[2]][fltr])
}
# PERIOD TESTS
vl_int = vl
if(isdate) vl_int = as.integer(vl)
d1 = vl_int[df$Year1 < ch_year[i]]
d2 = vl_int[df$Year1 >= ch_year[i]]
d1 = d1[!(is.infinite(d1) | is.na(d1) | is.nan(d1))]
d2 = d2[!(is.infinite(d2) | is.na(d2) | is.nan(d2))]
mean1[[i]] = round(mean(d1, na.rm = TRUE), 5)
mean2[[i]] = round(mean(d2, na.rm = TRUE), 5)
sd1[i] = round(sd(d1, na.rm = TRUE), 5)
sd2[i] = round(sd(d2, na.rm = TRUE), 5)
if(isdate) {
mean1[[i]] = mean1[[i]] %>% as.integer() %>% as.Date(origin = '1970-01-01')
mean2[[i]] = mean2[[i]] %>% as.integer() %>% as.Date(origin = '1970-01-01')
mratio[i] = lubridate::yday(mean2[[i]]) - lubridate::yday(mean1[[i]])
mean1[[i]] = mean1[[i]] %>% format("%d-%b")
mean2[[i]] = mean2[[i]] %>% format("%d-%b")
sd1[i] = as.integer(sd1[i])
sd2[i] = as.integer(sd2[i])
} else {
mratio[i] = round(100 * (mean2[[i]] - mean1[[i]]) / mean1[[i]], 1)
}
# Student and Fisher tests requre at least two observations in each sample
if (length(d1) > 1 & length(d2) > 1) {
if ((sum(abs(diff(d1)), na.rm = TRUE) != 0) && (sum(abs(diff(d2)), na.rm = TRUE) != 0)) {
tt[[i]] = t.test(d1, d2)
ft[[i]] = var.test(d1, d2)
}
}
}
cli::cli_progress_done()
pvalues = data.frame(
N = 1:nn,
Variable = desc,
Change.Year = ch_year,
Trend = sapply(ts_fit, function(X) { if(is.null(X)) NA else round(coef(X)[2], 5) }),
M1 = sapply(mean1, function(X) { if(is.null(X)) NA else X }),
M2 = sapply(mean2, function(X) { if(is.null(X)) NA else X }),
MeanRatio = mratio,
sd1 = sd1,
sd2 = sd2,
sdRatio = round(100 * (sd2 - sd1) / sd1, 1),
Mann.Kendall = sapply(mkt, function(X) { if(is.null(X)) NA else round(X$p.value, 5) }),
Pettitt = sapply(ptt, function(X) { if(is.null(X)) NA else round(X$p.value, 5) }),
Student = sapply(tt, function(X) { if(is.null(X)) NA else round(X$p.value, 5) }),
Fisher = sapply(ft, function(X) { if(is.null(X)) NA else round(X$p.value, 5) })
)
row.names(pvalues) = 1:nn
tests = list(ptt = ptt,
mkt = mkt,
tst = tst,
ts_fit = ts_fit,
tt = tt,
ft = ft,
year = ch_year,
maxval = maxval,
fixed_year = fixed,
pvalues = pvalues)
return(tests)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.