I would like to calculate the sens slope of a time serie organized by stations.
date station value
01/05/2007 1001 12.27
01/06/2007 1001 12.33
01/07/2007 1001 6.26
01/08/2007 1001 2.75
01/09/2007 1001 1.15
01/10/2007 1001 3.68
01/11/2007 1001 4.72
01/12/2007 1001 6.16
01/01/2008 1001 7.80
01/02/2008 1001 8.42
01/03/2008 1001 8.85
01/04/2008 1001 9.55
01/05/2008 1001 11.44
01/06/2008 1001 12.90
01/07/2008 1001 11.67
01/08/2008 1001 2.56
01/09/2008 1001 2.72
01/10/2008 1001 3.30
01/11/2008 1001 4.82
01/12/2008 1001 5.77
01/01/2009 1001 6.65
01/02/2009 1001 7.35
01/03/2009 1001 8.53
01/04/2009 1001 10.75
01/05/2009 1001 12.02
01/06/2009 1001 13.68
01/07/2009 1001 9.36
01/08/2009 1001 9.18
01/09/2009 1001 6.10
01/10/2009 1001 8.36
01/11/2009 1001 8.86
01/12/2009 1001 9.67
01/01/2010 1001 10.49
01/02/2010 1001 11.36
01/03/2010 1001 12.05
01/04/2010 1001 13.74
01/05/2010 1001 14.75
01/06/2010 1001 15.28
01/07/2010 1001 13.87
01/08/2010 1001 2.66
01/09/2010 1001 1.84
01/10/2010 1001 2.49
01/11/2010 1001 3.55
01/12/2010 1001 6.27
01/01/2011 1001 6.73
01/02/2011 1001 8.64
01/03/2011 1001 9.85
01/04/2011 1001 10.91
01/05/2011 1001 12.31
01/06/2011 1001 12.93
01/07/2011 1001 8.89
01/08/2011 1001 3.84
01/05/2016 1511 21.22
01/06/2016 1511 21.88
01/07/2016 1511 21.15
01/08/2016 1511 20.25
01/09/2016 1511 14.70
01/10/2016 1511 7.63
01/11/2016 1511 8.00
01/12/2016 1511 8.16
01/01/2017 1511 8.41
01/02/2017 1511 10.27
01/03/2017 1511 11.28
01/04/2017 1511 13.82
01/05/2017 1511 14.77
01/06/2017 1511 16.35
01/07/2017 1511 17.26
01/08/2017 1511 13.47
01/09/2017 1511 11.06
01/10/2017 1511 5.98
01/11/2017 1511 6.81
01/12/2017 1511 7.82
01/02/2018 1511 9.02
01/03/2018 1511 9.82
01/04/2018 1511 10.24
01/05/2018 1511 9.23
01/06/2018 1511 9.32
01/07/2018 1511 9.06
01/08/2018 1511 7.55
01/09/2018 1511 9.00
01/10/2018 1511 8.35
01/11/2018 1511 8.46
01/10/2007 1512 7.08
01/11/2007 1512 6.71
01/12/2007 1512 7.41
01/01/2008 1512 7.87
01/02/2008 1512 8.27
01/03/2008 1512 8.63
01/04/2008 1512 3.64
01/05/2008 1512 8.91
01/06/2008 1512 8.89
01/07/2008 1512 9.54
01/08/2008 1512 7.82
01/09/2008 1512 6.76
01/10/2008 1512 6.74
01/11/2008 1512 6.66
01/12/2008 1512 6.06
01/01/2009 1512 3.22
01/02/2009 1512 7.12
01/03/2009 1512 7.71
01/04/2009 1512 4.25
01/05/2009 1512 4.85
01/06/2009 1512 8.96
01/07/2009 1512 5.48
01/08/2009 1512 9.87
01/09/2009 1512 9.64
01/10/2009 1512 8.88
01/11/2009 1512 9.12
01/12/2009 1512 9.49
01/01/2010 1512 10.08
01/02/2010 1512 10.26
01/03/2010 1512 10.59
01/04/2010 1512 11.88
01/07/2017 1515 5.92
01/08/2017 1515 5.64
01/09/2017 1515 4.93
01/10/2017 1515 4.05
01/11/2017 1515 4.17
01/12/2017 1515 4.64
01/02/2018 1515 4.80
01/03/2018 1515 5.36
01/04/2018 1515 5.72
01/05/2018 1515 5.72
01/06/2018 1515 6.01
01/07/2018 1515 6.26
01/08/2018 1515 5.33
01/09/2018 1515 5.13
01/10/2018 1515 4.22
01/11/2018 1515 4.19
01/12/2018 1515 4.59
01/01/2019 1515 4.70
01/07/2019 1515 7.23
Hi, I have this dataframe and I m able to calculate Mann Kendall test with this code:
#function to perform Mann Kendall test
MK <- function(x){
MKfun <- possibly(Kendall::MannKendall,
list(tau = NA_real_,
sl = NA_real_,
S = NA_real_,
D = NA_real_,
varS = NA_real_))
out <- MKfun(x)
class(out) <- "list"
data.frame(out)}
#calling functions and store it in a dataframe
TS_MK <-
TS %>%
group_by(station) %>%
summarise(N = n(),
MK(value))
But I wanted to calculate the sen's slope and when I try this:
sens.slope(value)
I have this error:
Erreur : Problem with `summarise()` column `SEN`.
i `SEN = sens.slope(WL_m)`.
x `SEN` must be a vector, not a `htest` object.
i The error occurred in group 1: station = 1001.
Run `rlang::last_error()` to see where the error occurred.
> rlang::last_error()
<error/dplyr_error>
Problem with `summarise()` column `SEN`.
i `SEN = sens.slope(WL_m)`.
x `SEN` must be a vector, not a `htest` object.
i The error occurred in group 1: station = 1001.
How to calculate the sen slope? And in the same time than the Mann Kandall test if it's possible to have everything in the same data frame at the end. Thanks for your help.
------------------------UPDATE-----------------------------
Here is my last code version (thanks @AKRUN) and still it doesnt work:
library(trend)
library(broom)
library(tidyr)
library(dplyr)
library(Kendall)
#Loading data
TS = read.csv('Data/DATA.csv',sep=';')
#Calculate Mann Kendall test values for each station
TS %>%
mutate(value = readr::parse_number(value)) %>%
group_by(station) %>%
summarise(N = n(), out = list(MK(value))) %>%
unnest_wider(c(out))
#calculate Sen slope value for each station
TS %>%
mutate(value = readr::parse_number(value)) %>%
group_by(station) %>%
summarise(out = list(broom::tidy(sens.slope(value)))) %>%
unnest(out)
I have these errors:
Erreur : Problem with `mutate()` column `value`.
i `value = readr::parse_number(value)`.
x is.character(x) is not TRUE
Run `rlang::last_error()` to see where the error occurred.
I tried to change the value in character but still it doesnt work....
Thanks for your help!
CodePudding user response：
The 'value' column is not numeric class as there is , in the values. We can either remove the , and convert to numeric with as.numeric or use parse_number from readr to extract numeric part. Also, in the summarise, it may be better to wrap in a list (as the function returns a data.frame) which can be later unnested
library(dplyr)
library(tidyr)
TS %>%
mutate(value = readr::parse_number(value)) %>%
group_by(station) %>%
summarise(N = n(), out = list(MK(value))) %>%
unnest_wider(c(out))
-output
# A tibble: 4 x 7
station N tau sl S D varS
<int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1001 52 0.176 0.0660 234 1326. 16059.
2 1511 30 -0.361 0.00538 -157 435. 3142.
3 1512 31 0.351 0.00590 163 465. 3462.
4 1515 19 0.0352 0.861 6 170. 816
With sens.slope, the output is in a list structure. We can tidy up the output to return a tibble
library(broom)
library(trend)
TS %>%
mutate(value = readr::parse_number(value)) %>%
group_by(station) %>%
summarise(out = list(broom::tidy(sens.slope(value)))) %>%
unnest(out)
-output
# A tibble: 4 x 8
station statistic p.value parameter conf.low conf.high method alternative
<int> <dbl> <dbl> <int> <dbl> <dbl> <chr> <chr>
1 1001 1.84 0.0660 52 -0.351 15.4 Sen's slope two.sided
2 1511 -2.78 0.00538 30 -48.6 -8.5 Sen's slope two.sided
3 1512 2.75 0.00590 31 4.2 20.5 Sen's slope two.sided
4 1515 0.175 0.861 19 -7.18 13.5 Sen's slope two.sided
data
TS <- structure(list(date = c("01/05/2007", "01/06/2007", "01/07/2007",
"01/08/2007", "01/09/2007", "01/10/2007", "01/11/2007", "01/12/2007",
"01/01/2008", "01/02/2008", "01/03/2008", "01/04/2008", "01/05/2008",
"01/06/2008", "01/07/2008", "01/08/2008", "01/09/2008", "01/10/2008",
"01/11/2008", "01/12/2008", "01/01/2009", "01/02/2009", "01/03/2009",
"01/04/2009", "01/05/2009", "01/06/2009", "01/07/2009", "01/08/2009",
"01/09/2009", "01/10/2009", "01/11/2009", "01/12/2009", "01/01/2010",
"01/02/2010", "01/03/2010", "01/04/2010", "01/05/2010", "01/06/2010",
"01/07/2010", "01/08/2010", "01/09/2010", "01/10/2010", "01/11/2010",
"01/12/2010", "01/01/2011", "01/02/2011", "01/03/2011", "01/04/2011",
"01/05/2011", "01/06/2011", "01/07/2011", "01/08/2011", "01/05/2016",
"01/06/2016", "01/07/2016", "01/08/2016", "01/09/2016", "01/10/2016",
"01/11/2016", "01/12/2016", "01/01/2017", "01/02/2017", "01/03/2017",
"01/04/2017", "01/05/2017", "01/06/2017", "01/07/2017", "01/08/2017",
"01/09/2017", "01/10/2017", "01/11/2017", "01/12/2017", "01/02/2018",
"01/03/2018", "01/04/2018", "01/05/2018", "01/06/2018", "01/07/2018",
"01/08/2018", "01/09/2018", "01/10/2018", "01/11/2018", "01/10/2007",
"01/11/2007", "01/12/2007", "01/01/2008", "01/02/2008", "01/03/2008",
"01/04/2008", "01/05/2008", "01/06/2008", "01/07/2008", "01/08/2008",
"01/09/2008", "01/10/2008", "01/11/2008", "01/12/2008", "01/01/2009",
"01/02/2009", "01/03/2009", "01/04/2009", "01/05/2009", "01/06/2009",
"01/07/2009", "01/08/2009", "01/09/2009", "01/10/2009", "01/11/2009",
"01/12/2009", "01/01/2010", "01/02/2010", "01/03/2010", "01/04/2010",
"01/07/2017", "01/08/2017", "01/09/2017", "01/10/2017", "01/11/2017",
"01/12/2017", "01/02/2018", "01/03/2018", "01/04/2018", "01/05/2018",
"01/06/2018", "01/07/2018", "01/08/2018", "01/09/2018", "01/10/2018",
"01/11/2018", "01/12/2018", "01/01/2019", "01/07/2019"), station = c(1001L,
1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L,
1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L,
1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L,
1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L,
1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1001L,
1001L, 1001L, 1001L, 1001L, 1001L, 1001L, 1511L, 1511L, 1511L,
1511L, 1511L, 1511L, 1511L, 1511L, 1511L, 1511L, 1511L, 1511L,
1511L, 1511L, 1511L, 1511L, 1511L, 1511L, 1511L, 1511L, 1511L,
1511L, 1511L, 1511L, 1511L, 1511L, 1511L, 1511L, 1511L, 1511L,
1512L, 1512L, 1512L, 1512L, 1512L, 1512L, 1512L, 1512L, 1512L,
1512L, 1512L, 1512L, 1512L, 1512L, 1512L, 1512L, 1512L, 1512L,
1512L, 1512L, 1512L, 1512L, 1512L, 1512L, 1512L, 1512L, 1512L,
1512L, 1512L, 1512L, 1512L, 1515L, 1515L, 1515L, 1515L, 1515L,
1515L, 1515L, 1515L, 1515L, 1515L, 1515L, 1515L, 1515L, 1515L,
1515L, 1515L, 1515L, 1515L, 1515L), value = c("12,27", "12,33",
"6,26", "2,75", "1,15", "3,68", "4,72", "6,16", "7,80", "8,42",
"8,85", "9,55", "11,44", "12,90", "11,67", "2,56", "2,72", "3,30",
"4,82", "5,77", "6,65", "7,35", "8,53", "10,75", "12,02", "13,68",
"9,36", "9,18", "6,10", "8,36", "8,86", "9,67", "10,49", "11,36",
"12,05", "13,74", "14,75", "15,28", "13,87", "2,66", "1,84",
"2,49", "3,55", "6,27", "6,73", "8,64", "9,85", "10,91", "12,31",
"12,93", "8,89", "3,84", "21,22", "21,88", "21,15", "20,25",
"14,70", "7,63", "8,00", "8,16", "8,41", "10,27", "11,28", "13,82",
"14,77", "16,35", "17,26", "13,47", "11,06", "5,98", "6,81",
"7,82", "9,02", "9,82", "10,24", "9,23", "9,32", "9,06", "7,55",
"9,00", "8,35", "8,46", "7,08", "6,71", "7,41", "7,87", "8,27",
"8,63", "3,64", "8,91", "8,89", "9,54", "7,82", "6,76", "6,74",
"6,66", "6,06", "3,22", "7,12", "7,71", "4,25", "4,85", "8,96",
"5,48", "9,87", "9,64", "8,88", "9,12", "9,49", "10,08", "10,26",
"10,59", "11,88", "5,92", "5,64", "4,93", "4,05", "4,17", "4,64",
"4,80", "5,36", "5,72", "5,72", "6,01", "6,26", "5,33", "5,13",
"4,22", "4,19", "4,59", "4,70", "7,23")), class = "data.frame", row.names = c(NA,
-132L))