Step 1

Load libraries.

libraries <- c('Kendall', 'trend', 'changepoint')
lapply(libraries, require, character.only = TRUE)
## [[1]]
## [1] TRUE
## 
## [[2]]
## [1] TRUE
## 
## [[3]]
## [1] TRUE


Step 2

rain <- read.csv("D:/OneDrive - University of Vermont/Curriculum/19_ Personal webpage/TropicalFreshwaterEcology/_pages/Lectures/data/StreamEcology_TS_precipitation.csv")
rainEVFS <- rain[,3] # Select colunm 3.


Step 3

rain.ts <- ts(rainEVFS, # Convert "Precipitation" to a time series object.
                    frequency=12, start=c(1975,1), end=c(2016,12))
head(rain.ts, n=24) 
##  [1] 301.6  76.1  63.8 189.4 188.9  75.4 201.5 293.9 620.7 431.9 367.0 589.9
## [13] 158.3 272.8 192.2 134.0 215.9 394.4 179.0 278.7 444.2 368.6 190.1 205.4


Step 4

plot(rain.ts, main="Precipitation El Verde")


Mann-Kendall Trend Test

Step 5

We can use the Mann-Kendall Trend Test to see if there is a pattern in the data.

MannKendall(rain.ts)
## tau = 0.0301, 2-sided pvalue =0.31262
# tau is a measure of the strength of the trend
# p-value for trend is also shown

The test statistic is 0.0301, and the two-sided p-value associated with it is not less than 0.05. We cannot reject the null hypothesis of the test and conclude that no trend exists in the data because the p-value is greater than 0.05.


Step 6

Plot and trend Now we can add a smooth line to visualize the trend.

plot(rain.ts, main="Precipitation El Verde")
lines(lowess(time(rain.ts), rain.ts), col='red')


Sen’s slope

Step 7

The sens.slope function in the trend package is used with a time series object.
The sea.sens.slope function performs the test while taking into account the seasonality of the data.

sens.slope(rain.ts)
## 
##  Sen's slope
## 
## data:  rain.ts
## z = 1.0097, n = 504, p-value = 0.3126
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.04469697  0.13496333
## sample estimates:
## Sen's slope 
##  0.04663441


Pettitt’s test for change in value

Step 8

The pettitt.test function in the trend package identifies a point at which the values in the data change. The results here suggest that the stage values were higher after May 2009 than before.

pettitt.test(rain.ts)
## 
##  Pettitt's test for single change-point detection
## 
## data:  rain.ts
## U* = 7802, p-value = 0.116
## alternative hypothesis: two.sided
## sample estimates:
## probable change point at time K 
##                             244
rain[244,]
##     year month Rainfall
## 244 1995   abr   136.87


Multiple changes points

Step 9

The next few commands will be used to find a potential change point in the mean.

mvalueS = cpt.mean(rain.ts, method="BinSeg", Q=5) # BinSeg Binary Segmentation
cpts(mvalueS)
## [1]   7 189 255 420 468
rain[c(7,189,255,420,468),]
##     year month Rainfall
## 7   1975   jul   201.50
## 189 1990   sep   346.70
## 255 1996   mar   146.30
## 420 2009   dic   232.51
## 468 2013   dic   425.72
plot(mvalueS, cpt.width = 4)