In PGWG/analysis:

Using index_by_window()

Environmental data is often sparse and noisy. The non-paramatric Mann-Kendall test is a commonly employed technique for trend detection. For details see Chapter 16 of Gilbert (1987) Statistical Methods for Environmental Pollution Monitoring The Mann-Kendall test can be confounded by non-monotonic trends or negative trend approaching an asymptote. Because of the sensitivity of the test it is prudent to consider whether a trend exists over discrete intervals. We do not discuss it here, but one should pay attention to the issue of multiple hypothesis testing when many discrete intervals are considered. Environmental data is often sparse and noisy. The non-paramatric Mann-Kendall test is a commonly employed technique for trend detection. For details see Chapter 16 of Gilbert (1987) Statistical Methods for Environmental Pollution Monitoring The Mann-Kendall test can be confounded by non-monotonic trends or negative trend approaching an asymptote. Because of the sensitivity of the test it is prudent to consider whether a trend exists over discrete intervals. We do not discuss it here, but one should pay attention to the issue of multiple hypothesis testing when many discrete intervals are considered.

devtools::load_all()

# Toy data with a non-monotonic trend
set.seed(3)
dates <- seq(as.Date("1980/1/1"), as.Date("2005/1/1"), "quarters")
artificial_trend <- sin(seq(1,length(dates))/12)
artificial_trend_plus_noise <- artificial_trend  + rnorm(length(dates), 0, .5)
df <- data.frame(dates = dates, values = artificial_trend_plus_noise)
artificial_trend_plus_noise2 <- artificial_trend_plus_noise
artificial_trend_plus_noise2[60:101] <- -2
artificial_trend_plus_noise2 <- artificial_trend_plus_noise2 + rnorm(length(dates), 0, .1)
df <- data.frame(dates = dates, values = artificial_trend_plus_noise2)
# create a window_indexes
ind <- index_by_window(v= seq(1,length(dates)),by = 1, n = 32)
pvalues <- sapply(ind, function(x) Kendall::MannKendall(df[x,]$values)$s)
period_start <- lapply(ind, function(x) dates[x][1])
period_mid <- lapply(ind, function(x) dates[x][length(x)/2])
period_end<- lapply(ind, function(x) dates[x][length(x)])
period_start <-  do.call("c", period_start)
period_mid   <- do.call("c", period_mid)
period_end   <- do.call("c", period_end)

Here we vizualized the output of the sliding window Kendall::MannKendall(). Statistical significance decays after a long period of stabilitity

# Vizualization 
par(mfrow = c(2,1), mar = c(5,5,1,1) )
plot(dates, artificial_trend , type = "n", xlab = "", ylab = "", bty = "n", ylim = c(-2, 2))
points(dates, artificial_trend_plus_noise2 , type = "l", pch =20)
points(dates, artificial_trend_plus_noise2 , type = "p", pch =20)
plot(period_mid, pvalues, ylab = "Mann-Kendall p-value", xlab = "date", log = "y", bty="n",yaxt = "n", pch=2, main = "")
#arrows(period_mid, pvalues, period_end, pvalues, angle = 90, length = 0.01)
#arrows(period_mid, pvalues, period_start, pvalues, angle = 90, length = 0.01)
pow = c(-1,-2,-3,-4,-5,-6,-7)
ticksat <- as.vector(sapply(pow, function(p) (1:10)*10^p))
axis(2, ticksat, tcl= -.4, label = FALSE)
axis(2, 10^pow, tcl= 0.4, las = 2, labels = c(
  expression(10^-1),
  expression(10^-2),
  expression(10^-3),
  expression(10^-4),
  expression(10^-5),
  expression(10^-6),
  expression(10^-7)

))
abline(h = .05 , lty = 3)
text(as.Date( "1989-07-01"), 0.1, "p = 0.05 " )


# create a smaller window_indexes and vizualize teh difference
ind <- index_by_window(v= seq(1,length(dates)),by = 1, n = 16)
pvalues <- sapply(ind, function(x) Kendall::MannKendall(df[x,]$values)$s)
period_start <- lapply(ind, function(x) dates[x][1])
period_mid <- lapply(ind, function(x) dates[x][length(x)/2])
period_end<- lapply(ind, function(x) dates[x][length(x)])
period_start <-  do.call("c", period_start)
period_mid   <- do.call("c", period_mid)
period_end   <- do.call("c", period_end)

points(period_mid, pvalues, ylab = "Mann-Kendall p-value", xlab = "date", pch = 3)
legend("bottomright", pch = c(2,3), c("8 year window","4 year window"))

PGWG/analysis documentation built on May 7, 2019, 11:49 p.m.