如果时间序列和使用非 NA 数据，我如何限制 Mann Kendall 测试使用开始日期和结束日期？答案

【问题标题】：How do I limit a Mann Kendall test to use the start and end dates if the time series and to use non-NA data?如果时间序列和使用非 NA 数据，我如何限制 Mann Kendall 测试使用开始日期和结束日期？
【发布时间】：2020-09-22 15:48:05
【问题描述】：

我需要对不同长度的时间序列数据运行 Mann-Kendall 测试（R 中的包趋势，https://cran.r-project.org/web/packages/trend/index.html）。目前，时间序列分析是使用我手动指定的开始年份运行的，但这可能不是实际的开始日期。我的很多网站包含不同的开始年份，有些可能有不同的结束年份。我将我的数据压缩成以下内容。这是水质数据，因此存在数据缺失和开始/结束日期不同的问题。

我还在时间序列的中间和开头处理 NA。我想在时间序列中间消除丢失的 NA。如果 NA 在开头，我想从第一个实际值开始时间序列。


+---------+------------+------+--------------+-------------+-------------+---------------+--------------+
| SITE_ID | PROGRAM_ID | YEAR |  ANC_UEQ_L   |  NO3_UEQ_L  |  SO4_UEQ_L  | SBC_ALL_UEQ_L | SBC_NA_UEQ_L |
+---------+------------+------+--------------+-------------+-------------+---------------+--------------+
|    1234 | Alpha      | 1992 | 36.12        | 0.8786      | 91.90628571 | 185.5595714   | 156.2281429  |
|    1234 | Alpha      | 1993 | 22.30416667  | 2.671258333 | 86.85733333 | 180.5109167   | 154.1934167  |
|    1234 | Alpha      | 1994 | 25.25166667  | 3.296475    | 92.00533333 | 184.3589167   | 157.3889167  |
|    1234 | Alpha      | 1995 | 23.39166667  | 1.753436364 | 97.58981818 | 184.5251818   | 160.2047273  |
|    5678 | Beta       | 1983 | 4.133333333  | 20          | 134.4333333 | 182.1         | 157.4        |
|    5678 | Beta       | 1984 | 2.6          | 21.85       |      137.78 | 170.67        | 150.64       |
|    5678 | Beta       | 1985 | 3.58         | 20.85555556 | 133.7444444 | 168.82        | 150.09       |
|    5678 | Beta       | 1986 | -5.428571429 | 40.27142857 |       124.9 | 152.4         | 136.2142857  |
|    5678 | Beta       | 1987 | NA           | 13.75       |      122.75 | 137.4         | 126.3        |
|    5678 | Beta       | 1988 | 4.666666667  | 26.13333333 | 123.7666667 | 174.9166667   | 155.4166667  |
|    5678 | Beta       | 1989 | 6.58         | 31.91       |      124.63 | 167.39        | 148.68       |
|    5678 | Beta       | 1990 | 2.354545455  | 39.49090909 | 121.6363636 | 161.6454545   | 144.5545455  |
|    5678 | Beta       | 1991 | 5.973846154  | 30.54307692 | 119.8138462 | 165.4661185   | 147.0807338  |
|    5678 | Beta       | 1992 | 4.174359     | 16.99051285 | 124.1753846 | 148.5505115   | 131.8894862  |
|    5678 | Beta       | 1993 | 6.05         | 19.76125    |    117.3525 | 148.3025      | 131.3275     |
|    5678 | Beta       | 1994 | -2.51666     | 17.47167    |   117.93266 | 129.64167     | 114.64501    |
|    5678 | Beta       | 1995 | 8.00936875   | 22.66188125 |    112.3575 | 166.1220813   | 148.7095813  |
|    9101 | Victor     | 1980 | NA           | NA          |      94.075 | NA            | NA           |
|    9101 | Victor     | 1981 | NA           | NA          |       124.7 | NA            | NA           |
|    9101 | Victor     | 1982 | 33.26666667  | NA          | 73.53333333 | 142.75        | 117.15       |
|    9101 | Victor     | 1983 | 26.02        | NA          |        94.9 | 147.96        | 120.44       |
|    9101 | Victor     | 1984 | 20.96        | NA          |       82.98 | 137.4         | 110.46       |
|    9101 | Victor     | 1985 | 29.325       | 0.157843137 |      84.975 | 144.45        | 118.45       |
|    9101 | Victor     | 1986 | 28.6         | 0.88504902  |      81.675 | 139.7         | 114.45       |
|    9101 | Victor     | 1987 | 25.925       | 1.065441176 |       74.15 | 131.875       | 108.7        |
|    9101 | Victor     | 1988 | 29.4         | 1.048529412 |      80.625 | 148.15        | 122.5        |
|    9101 | Victor     | 1989 | 27.7         | 0.907598039 |      81.025 | 143.1         | 119.275      |
|    9101 | Victor     | 1990 | 27.4         | 0.642647059 |       77.65 | 126.825       | 104.775      |
|    9101 | Victor     | 1991 | 24.95        | 1.228921569 |        74.1 | 138.55        | 115.7        |
|    9101 | Victor     | 1992 | 29.425       | 0.591911765 |       73.85 | 130.675       | 106.65       |
|    9101 | Victor     | 1993 | 22.53333333  | 0.308169935 | 64.93333333 | 117.3666667   | 96.2         |
|    9101 | Victor     | 1994 | 29.93333333  | 0.428431373 | 67.23333333 | 124.0666667   | 101.2333333  |
|    9101 | Victor     | 1995 | 39.33333333  | 0.57875817  | 65.36666667 | 128.8333333   | 105.0666667  |
|    1121 | Charlie    | 1987 | 12.39        | 0.65        |       99.48 | 136.37        | 107.75       |
|    1121 | Charlie    | 1988 | 10.87333333  | 0.69        | 104.6133333 | 131.9         | 105.2        |
|    1121 | Charlie    | 1989 | 5.57         | 1.09        |      105.46 | 136.125       | 109.5225     |
|    1121 | Charlie    | 1990 | 13.4725      | 0.8975      |      99.905 | 134.45        | 108.9875     |
|    1121 | Charlie    | 1991 | 11.3         | 0.805       |     100.605 | 134.3775      | 108.9725     |
|    1121 | Charlie    | 1992 | 9.0025       | 7.145       |      99.915 | 136.8625      | 111.945      |
|    1121 | Charlie    | 1993 | 7.7925       | 6.6         |      95.865 | 133.0975      | 107.4625     |
|    1121 | Charlie    | 1994 | 7.59         | 3.7625      |     97.3575 | 129.635       | 104.465      |
|    1121 | Charlie    | 1995 | 7.7925       | 1.21        |      100.93 | 133.9875      | 109.5025     |
|    3812 | Charlie    | 1988 | 18.84390244  | 17.21142857 | 228.8684211 | 282.6540541   | 260.5648649  |
|    3812 | Charlie    | 1989 | 11.7248      | 21.21363636 | 216.5973451 | 261.3711712   | 237.4929204  |
|    3812 | Charlie    | 1990 | 2.368571429  | 35.23448276 | 216.7827586 | 286.0034483   | 264.3137931  |
|    3812 | Charlie    | 1991 | 33.695       | 40.733      |      231.92 | 350.91075     | 328.443      |
|    3812 | Charlie    | 1992 | 18.49111111  | 26.14818889 |    219.1488 | 301.3785889   | 281.8809222  |
|    3812 | Charlie    | 1993 | 17.28181818  | 27.65394545 | 210.6605091 | 290.064       | 271.9205455  |
+---------+------------+------+--------------+-------------+-------------+---------------+--------------+

如果我将开始年份更改为错过早期数据中的 NA，则这是当前将为我的实际数据运行时间序列的代码。它适用于在整个时间段内都有价值的网站，但在考虑不同的开始/结束年份时会给出奇怪的结果。

Mann_Kendall_Values_Trimmed <- filter(LTM_Data_StackOverflow_9_22_2020, YEAR >1984) %>% #I manually trimmed the data here to prevent some errors
  group_by(SITE_ID) %>% 
  filter(n() > 2) %>% #filter sites with more than 10 years of data
  gather(parameter, value, SO4_UEQ_L, ANC_UEQ_L, NO3_UEQ_L, SBC_ALL_UEQ_L, SBC_NA_UEQ_L ) %>% 
  #, DOC_MG_L) 
  group_by(parameter, SITE_ID, PROGRAM_ID) %>% nest() %>% 
  mutate(ts_out = map(data, ~ts(.x$value, start=c(1985, 1), end=c(1995, 1), frequency=1))) %>% 
#this is where I would like to specify the first year in the actual time series with data. End year would also be tied to the last year of data.
  mutate(mk_res = map(ts_out, ~mk.test(.x, alternative = c("two.sided", "greater", "less"),continuity = TRUE)),
         sens = map(ts_out, ~sens.slope(.x, conf.level = 0.95))) %>% 
  #run the Mann-Kendall Test
  mutate(mk_stat = map_dbl(mk_res, ~.x$statistic), 
         p_val = map_dbl(mk_res, ~.x$p.value)
         , sens_slope = map_dbl(sens, ~.x$estimates)
  ) %>% 
  #Pull the parameters we need
  select(SITE_ID, PROGRAM_ID, parameter, sens_slope, p_val, mk_stat) %>% 
  mutate(output = case_when(
  sens_slope == 0 ~ "NC",
  sens_slope > 0 & p_val < 0.05 ~ "INS",
  sens_slope > 0 & p_val > 0.05 ~ "INNS",
  sens_slope < 0 & p_val < 0.05 ~ "DES",
  sens_slope < 0 & p_val > 0.05 ~ "DENS"))

如何处理数据中间的 NA？
如何让时间序列在实际数据的日期自动开始和结束？作为参考，每个 site_id 都有以下日期范围（不包括 NA）：

+-----------+-----------+-------------------+-----------+-----------+
|   1234    |   5678    |       9101        |   1121    |   3812    |
+-----------+-----------+-------------------+-----------+-----------+
| 1992-1995 | 1983-1995 | 1982 OR 1985-1995 | 1987-1995 | 1988-1993 |
+-----------+-----------+-------------------+-----------+-----------+

【问题讨论】：

由于测试是基于等级的，平滑/插值中间值 a) 是否有帮助，b) 有必要吗？似乎平滑内部 NA 只是增加了您稍后要测试的趋势。并且由于测试是非参数的且基于等级的，因此省略一个点会有所不同吗？如果系列是单调的（ish），无论您是否跳过一年，排名都会给出正确的答案。
简短版：我认为省略 NA 值是最明智的做法。
@Jason 我对数据进行了更深入的研究，似乎所有缺失值都发生在 1990 年之前，所以我删除了 1990 年之前的数据（这类似于我们在 Excel Mann-Kendall 中所做的测试）。然后我可以参考代码（参见下面的答案）来参考 .x$YEAR 的最小值和最大值，由于参数的分组，我以前无法做到这一点。

标签： r time-series trend

【解决方案1】：

为了使数据更加一致，我决定在导入 R 之前在 Oracle 中将数据组织为单独的时间序列（按参数、年份、站点 ID、程序分组）。

+---------+------------+------+--------------+-----------+
| SITE_ID | PROGRAM_ID | YEAR |    Value     | Parameter |
+---------+------------+------+--------------+-----------+
|    1234 | Alpha      | 1992 | 36.12        | ANC       |
|    1234 | Alpha      | 1993 | 22.30416667  | ANC       |
|    1234 | Alpha      | 1994 | 25.25166667  | ANC       |
|    1234 | Alpha      | 1995 | 23.39166667  | ANC       |
|    5678 | Beta       | 1990 | 2.354545455  | ANC       |
|    5678 | Beta       | 1991 | 5.973846154  | ANC       |
|    5678 | Beta       | 1992 | 4.174359     | ANC       |
|    5678 | Beta       | 1993 | 6.05         | ANC       |
|    5678 | Beta       | 1994 | -2.51666     | ANC       |
|    5678 | Beta       | 1995 | 8.00936875   | ANC       |
|    9101 | Victor     | 1990 | 27.4         | ANC       |
|    9101 | Victor     | 1991 | 24.95        | ANC       |
|    9101 | Victor     | 1992 | 29.425       | ANC       |
|    9101 | Victor     | 1993 | 22.53333333  | ANC       |
|    9101 | Victor     | 1994 | 29.93333333  | ANC       |
|    9101 | Victor     | 1995 | 39.33333333  | ANC       |
|    1121 | Charlie    | 1990 | 13.4725      | ANC       |
|    1121 | Charlie    | 1991 | 11.3         | ANC       |
|    1121 | Charlie    | 1992 | 9.0025       | ANC       |
|    1121 | Charlie    | 1993 | 7.7925       | ANC       |
|    1121 | Charlie    | 1994 | 7.59         | ANC       |
|    1121 | Charlie    | 1995 | 7.7925       | ANC       |
|    3812 | Charlie    | 1990 | 2.368571429  | ANC       |
|    3812 | Charlie    | 1991 | 33.695       | ANC       |
|    3812 | Charlie    | 1992 | 18.49111111  | ANC       |
|    3812 | Charlie    | 1993 | 17.28181818  | ANC       |
|    1234 | Alpha      | 1992 | 0.8786       | NO3       |
|    1234 | Alpha      | 1993 | 2.671258333  | NO3       |
|    1234 | Alpha      | 1994 | 3.296475     | NO3       |
|    1234 | Alpha      | 1995 | 1.753436364  | NO3       |
|    5678 | Beta       | 1990 | 39.49090909  | NO3       |
|    5678 | Beta       | 1991 | 30.54307692  | NO3       |
|    5678 | Beta       | 1992 | 16.99051285  | NO3       |
|    5678 | Beta       | 1993 | 19.76125     | NO3       |
|    5678 | Beta       | 1994 | 17.47167     | NO3       |
|    5678 | Beta       | 1995 | 22.66188125  | NO3       |
|    9101 | Victor     | 1990 | 0.642647059  | NO3       |
|    9101 | Victor     | 1991 | 1.228921569  | NO3       |
|    9101 | Victor     | 1992 | 0.591911765  | NO3       |
|    9101 | Victor     | 1993 | 0.308169935  | NO3       |
|    9101 | Victor     | 1994 | 0.428431373  | NO3       |
|    9101 | Victor     | 1995 | 0.57875817   | NO3       |
|    1121 | Charlie    | 1990 | 0.8975       | NO3       |
|    1121 | Charlie    | 1991 | 0.805        | NO3       |
|    1121 | Charlie    | 1992 | 7.145        | NO3       |
|    1121 | Charlie    | 1993 | 6.6          | NO3       |
|    1121 | Charlie    | 1994 | 3.7625       | NO3       |
|    1121 | Charlie    | 1995 | 1.21         | NO3       |
|    3812 | Charlie    | 1990 | 35.23448276  | NO3       |
|    3812 | Charlie    | 1991 | 40.733       | NO3       |
|    3812 | Charlie    | 1992 | 26.14818889  | NO3       |
|    3812 | Charlie    | 1993 | 27.65394545  | NO3       |
+---------+------------+------+--------------+-----------+

在 R 中，我能够将代码编辑为以下代码开头的代码。其余代码相同。

Mann_Kendall_Values_Trimmed <- filter(LTM_Data_StackOverflow_9_22_2020, YEAR >1989, PARAMETER != 'doc') %>% 
  #filter data to start in 1990 as this removes nulls from pre-1990 sampling
  group_by(SITE_ID) %>% 
  filter(n() > 10) %>% #filter sites with more than 10 years of data
  #gather(SITE_ID, PARAMETER, VALUE) #I believe this is now redundant %>% 
  group_by(PARAMETER, SITE_ID, PROGRAM_ID) %>% nest() %>% 
  mutate(ts_out = map(data, ~ts(.x$VALUE, start=c(min(.x$YEAR), 1), c(max(.x$YEAR), 1), frequency=1)))

这达到了我需要的所有时间序列的结果，这些时间序列有足够的长度（我相信大于 2）来运行 mann-kendall 测试。有这些问题的参数将在单独的 R 代码中处理。

【讨论】：