Как вы выполняете тест Манна-Кендалла на нескольких станциях, используя петли? - PullRequest
Новая
       23

Как вы выполняете тест Манна-Кендалла на нескольких станциях, используя петли?

1 голос
/ 03 апреля 2019

Я только начал использовать R и хотел бы использовать модифицированный пакет mk для тестирования месячных данных об уровне подземных вод. Мой фрейм данных (GL) выглядит примерно так 

GL

well    year    month   value
684     1994    Jan     8.53
684     1995    Jan     8.74
684     1996    Jan     8.88
684     1997    Jan     8.24
1001    2000    Jan     9.1
1001    2001    Jan     9.2
1001    2002    Jan     9.54
1001    2003    Jan     9.68
2003    1981    Jan     55.2
2003    1982    Jan     55.8
2003    1983    Jan     56.4
2003    1984    Jan     53.2

Сначала я создал список лунок и файл результатов для печати результатов 

well_list <- unique(GL$well)
results_list <- vector("list", length(well_list))

Тогда я создал то, что я считаю циклом 

for(i in well_list){
  results_list[[grep(i, well_list)]] <- MannKendall(GL[,4])
}
names(results_list) <- well_list

но я продолжал получать эту ошибку 

Error in Kendall(1:length(x), x) : length(x)<3

Я могу заставить этот код работать (который я взял из другого поста здесь), но я не хочу выполнять метод предварительного отбеливания 

for(i in well_list){
tempDF <- GL1[GL$well == i, ]
c<-acf(tempDF$value,lag.max=1)
t <- dim(c$acf)
tempDF$prewhit1<-c$acf[[t[1], t[2], t[3]]]*tempDF$value
prewhitseries<-data.frame(with(tempDF,(tempDF$value[-1] - prewhit1[ 
length(prewhit1)])))
autocordata<-cbind(tempDF[-1,],prewhitseries)
results_list[[grep(i, well_list)]] <- MannKendall(autocordata[,5])}
names(results_list) <- well_list

Ответы [ 2 ]

1 голос
/ 03 апреля 2019

Нечто подобное должно делать.split() разбивается по столбцу well, создавая список с вектором для каждой лунки.Сохраняются только векторы длиной 3 или более.MannKendall() затем запускается на каждом из оставшихся векторов, используя lapply()

library(Kendall)

tt <- read.table(text="
well    year    month   value
684     1994    Jan     8.53
684     1995    Jan     8.74
684     1996    Jan     8.88
684     1997    Jan     8.24
1001    2000    Jan     9.1
1001    2001    Jan     9.2
1001    2002    Jan     9.54
1001    2003    Jan     9.68
2003    1981    Jan     55.2
2003    1982    Jan     55.8
2003    1983    Jan     56.4
2003    1984    Jan     53.2
2004    1984    Jan     53.2", header=TRUE)

tt.wells <- split(tt$value, tt$well)
tt.wells <- tt.wells[lengths(tt.wells) >= 3]

lapply(tt.wells, MannKendall)

# $`684`
# tau = 0, 2-sided pvalue =1

# $`1001`
# tau = 1, 2-sided pvalue =0.089429

# $`2003`
# tau = 0, 2-sided pvalue =1
0 голосов
/ 03 апреля 2019

Может быть, с split/lapply это проще и удобочитаемее.

Сначала запускаются тесты. 

sp <- split(GL, GL$well)
results_list <- lapply(sp, function(DF){
  tryCatch(MannKendall(DF[, 4]),
           error = function(e) e)
})

Теперь найдите хорошие, без ошибок. 

bad <- sapply(results_list, inherits, "error")

И осмотреть их. 

results_list[bad]
#named list()

results_list[!bad]
#$`684`
#tau = 0, 2-sided pvalue =1
#
#$`1001`
#tau = 1, 2-sided pvalue =0.089429
#
#$`2003`
#tau = 0, 2-sided pvalue =1
...