Skip to content
GitLab

Included the modified-Mann-Kendall significance test for Trend

Open Carlos Delgado Torres requested to merge develop-Trend-sign into master

Hi @aho and @erifarov,

As discussed in today's meeting, I have included the modified Mann-Kendall test to assess the statistical significance of the trend (described in Hamed and Rao 1998; doi:10.1016/S0022-1694(97)00125-X). For it, I have added the parameter sign.test in Trend(). The options for this new parameter are "ANOVA" (the test that was already included) and "modified-Mann-Kendall" (the new test). The default value is "ANOVA" for reproducibility.

@vtorralba is using the function and kindly agreed on testing it.

I have already checked that the results are identical to the previous version if sign.test="ANOVA", and plotted a figure for comparing both methods:

library(startR)
library(s2dv)
library(multiApply)
source('/esarchive/scratch/cdelgado/gitlab/s2dv/R/Trend.R')

data <- startR::Start(dat = '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h/$var$_$year$$month$.nc',
                      var = 'tas',
                      year = paste0(1981:2015),
                      lon = values(list(-85,45)),
                      lat = values(list(27,81)),
                      month = '01',
                      return_vars = list(lat = NULL, lon = NULL), 
                      lon_var = 'lon', 
                      lat_var = 'lat',
                      lat_reorder = Sort(decreasing = FALSE),
                      lon_reorder = CircularSort(-180,180),
                      num_procs = 1, 
                      retrieve = TRUE)
lat <- attr(data,'Variables')$common$lat
lon <- attr(data,'Variables')$common$lon
data <- drop(data)

trend_anova_old <- s2dv::Trend(data = data, time_dim = 'year', ncores = 4)
trend_anova_new <- Trend(data = data, time_dim = 'year', sign.test = "ANOVA")
trend_mk_new <- Trend(data = data, time_dim = 'year', sign.test = "modified-Mann-Kendall")

identical(trend_anova_old,trend_anova_new) # TRUE

s2dv::PlotLayout(fun = PlotEquiMap, plot_dims = c('lon','lat'), lon = lon, lat = lat, ncol = 2, 
                 var = list(trend_anova_new$trend[2,,], trend_mk_new$trend[2,,]),
                 special_args = list(list(dots = trend_anova_new$p.val <= 0.05), 
                                     list(dots = trend_mk_new$p.val <= 0.05)),
                 # toptitle = title_metrics, title_margin_scale = 3, title_scale = 1.5, 
                 titles = c('Student-T test', 'Modified Mann-Kendall'), sizetit = 1.5,
                 brks = seq(-0.24, 0.24, by = 0.04), bar_extra_labels = seq(-0.24, 0.24, by = 0.04), triangle_ends = c(T,T),
                 bar_scale = 0.8, bar_label_scale = 3, axelab = FALSE, draw_separators = FALSE, colNA = 'grey', 
                 country.borders = TRUE, coast_width = 2, dot_symbol = 4, dot_size = 4, filled.continents = FALSE, filled.oceans = FALSE, 
                 res = 50, width = 30, height = 10, fileout = '/esarchive/scratch/cdelgado/test.png')

image

Thank you,
Carlos