-
Notifications
You must be signed in to change notification settings - Fork 1
/
wrtdstrnd_sk.R
210 lines (180 loc) · 7.11 KB
/
wrtdstrnd_sk.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
######
#' Get WRTDS trends using seasonal Kendall tests
#'
#' Get WRTDS trends using seasonal Kendall tests
#'
#' @param dat_in input tidal or tidalmean object which must already have fitted model data
#' @param mobrks list of month groupings where each month is an integer from 1 to 12, see examples
#' @param yrbrks numeric vector of breaks for years, see examples
#' @param molabs character vector of names for month breaks, see examples
#' @param yrlabs character vector of names for year breaks, see examples
#' @param tau numeric vector of quantile for estimating trends
#' @param trndvar chr string of variable for trend evaluation, usually back-transformed, flow-normalized results, see details
#' @param ... methods passed to or from other methods
#'
#' @export
#'
#' @details Trends are based on \code{\link{kendallSeasonalTrendTest}} for user-specified time periods. In general, the seasonal Kendall test evaluates monotonic trends using a non-parametric approach that accounts for seasonal variation in the time series.
#'
#' All trends are based on back-transformed, flow-normalized results by default. The variable for evaluating trends can be changed with \code{'trndvar'} as \code{'res'}, \code{'norm'}, or \code{'fit'} for \code{tidal} objects and as \code{'res'}, \code{'bt_norm'}, or \code{'bt_fits'} for \code{tidalmean} objects. In all cases, back-transformed variables are evaluated.
#'
#' The user must supply the annual and monthly aggregation periods to the appropriate arguments. These are passed to \code{\link[base]{cut}} and are left-open, right-closed along the interval.
#'
#' @import dplyr
#' @importFrom lubridate month year
#' @importFrom purrr map
#' @importFrom tidyr nest unnest
#'
#' @return A \code{\link[base]{data.frame}} with summary trends for each grouping, including \code{med} as the median value for the period of observation, \code{tau} as the magnitude and direction of the trend, \code{slope} as the Thiel-Sen slope for change per year, \code{chitest} as the significance test evaluating heterogeneity between seasons, \code{ztest} indicating significance of the overall trend, and \code{perchg} as 100 multiplied by the ratio of the annual slope to the median estimate of the time period (percent change per year).
#'
#' As noted in \code{\link{kendallSeasonalTrendTest}}, the overall test is not appropriate if \code{chitest} indicates a small p-value.
#'
#' @references Hirsch, R.M., Slack, J.R., Smith, R.A. 1982. Techniques of trend analysis for monthly water quality data. Water Resources Research, 18:107-121.
#'
#' Millard, S. P. 2013. EnvStats: An R Package for Environmental Statistics. Springer, New York.
#'
#' @examples
#' ## load a fitted model object
#' data(tidfit)
#' data(tidfitmean)
#'
#' ## get trends
#'
#' # setup month, year categories
#' mobrks <- list(c(1, 2, 3), c(4, 5, 6), c(7, 8, 9), c(10, 11, 12))
#' yrbrks <- c(1973, 1985, 1994, 2003, 2012)
#' molabs <- c('JFM', 'AMJ', 'JAS', 'OND')
#' yrlabs <- c('1974-1985', '1986-1994', '1995-2003', '2004-2012')
#'
#' wrtdstrnd_sk(tidfit, mobrks, yrbrks, molabs, yrlabs)
#' wrtdstrnd_sk(tidfitmean, mobrks, yrbrks, molabs, yrlabs)
#'
wrtdstrnd_sk <- function(dat_in, ...) UseMethod('wrtdstrnd_sk')
#' @rdname wrtdstrnd_sk
#'
#' @export
#'
#' @method wrtdstrnd_sk default
wrtdstrnd_sk.default <- function(dat_in, mobrks, yrbrks, molabs, yrlabs, ...){
# create year mo cats
dat_in <- select(dat_in, date, trndvar) %>%
mutate(
yr = year(date),
mo = month(date)
)
# annual aggs
yrdat <- mutate(dat_in,
yrcat = cut(yr, breaks = yrbrks, labels = yrlabs)
) %>%
na.omit %>%
rename(cat = yrcat) %>%
group_by(cat) %>%
nest %>%
mutate(
sk_res = map(data, ~ kendallSeasonalTrendTest(trndvar ~ mo + yr, .)),
med = map(data, ~ median(.$trndvar, na.rm = TRUE)),
tau = map(sk_res, ~ .$estimate[1]),
slope = map(sk_res, ~ .$estimate[2]),
chitest = map(sk_res, ~ .$p.value[1]),
ztest = map(sk_res, ~ .$p.value[2])
) %>%
select(-data, -sk_res) %>%
unnest
# setup monthly categories
if(length(molabs) != length(mobrks))
stop('molabs are not matched to mobrks')
names(mobrks) <- molabs
mobrks <- as.data.frame(mobrks) %>%
tidyr::gather(.) %>%
arrange(value)
# monthly aggs
modat <- mutate(dat_in,
mocat = factor(mo)
)
levels(modat$mocat) <- mobrks$key
modat <- group_by(modat, yr, mo, mocat) %>%
summarise(trndvar = mean(trndvar, na.rm = T)) %>%
arrange(yr, mo) %>%
rename(cat = mocat) %>%
group_by(cat) %>%
nest %>%
mutate(
sk_res = map(data, ~ kendallSeasonalTrendTest(trndvar ~ mo + yr, .)),
med = map(data, ~ median(.$trndvar, na.rm = TRUE)),
tau = map(sk_res, ~ .$estimate[1]),
slope = map(sk_res, ~ .$estimate[2]),
chitest = map(sk_res, ~ .$p.value[1]),
ztest = map(sk_res, ~ .$p.value[2])
) %>%
select(-data, -sk_res) %>%
unnest
# combine results
out <- rbind(yrdat, modat) %>%
ungroup %>%
mutate(
cat = factor(cat,
levels = c(yrlabs, molabs),
labels = c(yrlabs, molabs)
),
perchg = 100 * slope/med
) %>%
arrange(cat) %>%
data.frame
return(out)
}
#' @rdname wrtdstrnd_sk
#'
#' @export
#'
#' @method wrtdstrnd_sk tidal
wrtdstrnd_sk.tidal <- function(dat_in, mobrks, yrbrks, molabs, yrlabs, tau = NULL, trndvar = 'norm', ...){
# stop if trndvar not one of these
if(!trndvar %in% c('res', 'norm', 'fit'))
stop('var must be one of res, norm, or fit')
# get names if trndvar isn't res
if(trndvar != 'res'){
# get tau if null, otherwise run checks
if(is.null(tau)){
fstr <- paste0('^', trndvar)
tau <- grep(fstr, names(dat_in))
tau<- gsub(fstr, '', names(dat_in)[floor(median(tau))])
toget <- c('^date$', paste0(fstr, tau, '$'))
} else {
if(length(grep(paste0(tau, '$'), names(dat_in))) == 0)
stop('Specified tau not in object')
}
} else {
toget <- c('^date$', paste0('^', trndvar, '$'))
}
# columns to get with regex
toget <- paste(toget, collapse = '|')
# select columns, format date as year, month
# back transform
toeval <- select(dat_in, matches(toget))
names(toeval) <- c('date', 'trndvar')
toeval <- mutate(toeval,
trndvar = exp(trndvar)
)
wrtdstrnd_sk.default(toeval, mobrks, yrbrks, molabs, yrlabs, ...)
}
#' @rdname wrtdstrnd_sk
#'
#' @export
#'
#' @method wrtdstrnd_sk tidalmean
wrtdstrnd_sk.tidalmean <- function(dat_in, mobrks, yrbrks, molabs, yrlabs, trndvar = 'bt_norm', ...){
# stop if trndvar not one of these
if(!trndvar %in% c('res', 'bt_norm', 'bt_fits'))
stop('var must be one of res, bt_norm, or bt_fits')
# columns to get with regex
toget <- paste0('^date$|^', trndvar, '$')
# select columns, format date as year, month
toeval <- select(dat_in, matches(toget))
names(toeval) <- c('date', 'trndvar')
# need to back-transform res if used
if(trndvar == 'res')
toeval <- mutate(toeval,
trndvar = exp(trndvar)
)
wrtdstrnd_sk.default(toeval, mobrks, yrbrks, molabs, yrlabs, ...)
}