Блог

Главная - Вопросы по программированию - Как вы выполняете тест Манна-Кендалла на нескольких станциях с использованием циклов?

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

#r #loops #time-series #lapply

#r #циклы #временные ряды #lapply

Вопрос:

Я только начал использовать R и хотел бы использовать пакет modifiedmk для выполнения тестов на ежемесячных данных об уровне грунтовых вод. Мой фрейм данных (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

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

 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

Ответ №1:

Что-то вроде этого должно сработать. 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

Комментарии:

1. Привет, ребята, это было потрясающе! первое решение сработало, но, к сожалению, у меня все еще были проблемы с ошибками. Второе решение сработало, но без ошибок. Огромное спасибо — я целую вечность пытался разобраться в этом, и посмотрите, как это было просто!

Ответ №2:

Возможно, с 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