Ваш код можно ускорить, используя join s для сопоставления id-пар в dat2 с измерениями в dat1 вместо for -циклов. Кстати: по крайней мере, на мой взгляд, использование объединений делает код также более лаконичным и понятным. И более надежный.
Примечание далее : Я нашел ошибку в вашем коде. Настройка матрицы корреляций с использованием correlations[x,1] приводит к присвоению корреляций неправильным id-парам.
Тест
Для сравнения с вашим подходом я настроил две функции: cor_join, которая реализует простой подход к этой идее, и cor_loop, который является оболочкой для вашего кода.
Микробенчмаркинг обеих функций показывает, что использование объединений ускоряет вычисления в ~ 2,5 раза. Я не знаком с data.table, но могу предположить, что использование data.table или dtplyr (бэкэнд таблицы данных для dplyr) может еще больше повысить производительность, особенно для вашего реального набора данных.
<!-- language-all: lang-r -->
library(data.table)
library(dplyr)
library(purrr)
library(ggplot2)
library(microbenchmark)
dat1 <- fread("https://www.dropbox.com/s/d2s61du255vzu7g/dat1.csv?dl=1") # ~2 MB
dat2 <- fread("https://www.dropbox.com/s/7n0z0gbeoifss4j/dat2.csv?dl=1") # ~5 KB
# fix column classes
dat1$date <- as.Date(dat1$date)
dat1$ID <- as.character(dat1$ID)
dat2[, (c("common_date_begin","common_date_end")) := lapply(.SD, as.Date), .SDcols = c("common_date_begin","common_date_end")]
dat2[, (c("ID1","ID2")) := lapply(.SD, as.character), .SDcols = c("ID1","ID2")]
cor_join <- function(dat1, dat2) {
# We want to get a dataframe with
#
# 1. pairs of sites,
# 2. dates where we have measurements for both
# 3. the measurements at each site
#
# This could be achieved via left_joins
dat3 <- dat2 %>%
# Join dates and measurements for ID1
left_join(dat1, by = c("ID1" = "ID")) %>%
rename(value1 = value) %>%
# Join dates and measurements for ID2 on the same date
left_join(dat1, by = c("ID2" = "ID", "date" = "date")) %>%
rename(value2 = value, ID = ID1, ID_nearest = ID2)
dat3
# Compute correlations
dat3 %>%
# Drop missings, i.e. observations with no common dates
filter(date >= common_date_begin & date <= common_date_end) %>%
group_by(ID, ID_nearest, dist, common_date_begin, common_date_end, diff_days) %>%
summarise(correl = cor(value1, value2)) %>%
ungroup()
}
cor_loop <- function(dat1, dat2) {
# get list of unique stations
ids <- unique(dat2$ID1)
# initialize matrix to hold correlations
correlations <- matrix(NA, nrow = nrow(dat2), ncol=1)
# initialize data frame to hold results
results <- as.data.frame(dat2[, -c(4:5)])
# initialize loop counters
x <- 1
# loop over the main ID's
for (i in ids) {
tmp <- dat2[ID1==i]
#loop over the ID's of the neighbour stations
for (id in 1:nrow(tmp)){
# get ID of the neighbours
near_id <- as.numeric(tmp[id, 2])
# get common dates
beg_date <- tmp[id, 4]
end_date <- tmp[id, 5]
# calculate correlations
correlations[x,1] <- cor(dat1[ID==i & date %between% c(beg_date, end_date)]$value,
dat1[ID==near_id & date %between% c(beg_date, end_date)]$value)
# increment loop counter
x <- x + 1
}
}
# assemble final data frame
results <- data.table(ID=results[, 1],
ID_nearest=results[, 2],
distance=results[, 3],
overlapping_days=results[, 4],
correl=as.vector(correlations))
results
}
# microbenchmark
microbenchmark::microbenchmark(cor_join(dat1, dat2), cor_loop(dat1, dat2), times = 10)
#> Unit: milliseconds
#> expr min lq mean median uq max
#> cor_join(dat1, dat2) 247.4106 286.1556 301.6367 296.6921 302.2751 400.8654
#> cor_loop(dat1, dat2) 773.5274 784.9197 807.3767 798.4800 842.3080 854.1716
#> neval
#> 10
#> 10
Результаты проверки
Чтобы проверить, что обе функции дают одинаковые результаты, я создал диаграмму рассеяния
# Check result
results <- list(join = cor_join(dat1, dat2), loop = cor_loop(dat1, dat2))
# Plot
check <- results %>%
purrr::reduce(left_join, by = c("ID", "ID_nearest"), suffix = c("_join", "_loop"))
check %>%
ggplot(aes(correl_join, correl_loop, color = ID)) +
geom_point()
OOPS: диаграмма рассеяния показывает разные результаты? Чтобы проверить, что я использовал простой набор данных, где я отфильтровал наборы данных для сайтов 1183, 1550 и 1551:
dat1a <- dat1 %>% filter(ID %in% c(1183, 1550, 1551)) %>% as.data.table()
dat2a <- dat2 %>% filter(ID1 %in% c(1183, 1550, 1551), ID2 %in% c(1183, 1550, 1551)) %>% as.data.table()
# For the simple dataset I get the same correlations
cor_join(dat1a, dat2a)
#> # A tibble: 2 x 7
#> ID ID_nearest dist common_date_begin common_date_end diff_days correl
#> <chr> <chr> <dbl> <date> <date> <int> <dbl>
#> 1 1183 1550 1576360. 2010-02-23 2017-06-16 2670 0.0456
#> 2 1183 1551 1513356. 2010-02-23 2017-06-16 2670 -0.0251
cor_loop(dat1a, dat2a)
#> ID ID_nearest distance overlapping_days correl
#> 1: 1183 1550 1576360 2670 0.04564506
#> 2: 1183 1551 1513356 2670 -0.02513991
После проверки вашего кода я догадался, что различия возникают из-за присвоения корреляций неправильным id-парам из-за correlations[x,1]. Для проверки я настроил cor_loop. Помимо возврата df results он также возвращает второй df correlations2, который установлен в l oop и содержит не только корреляцию, но и соответствующее значение id и near_id:
cor_loop_check <- function(dat1, dat2) {
# get list of unique stations
ids <- unique(dat2$ID1)
# initialize matrix to hold correlations
correlations <- matrix(NA, nrow = nrow(dat2), ncol=1)
correlations2 <- data.frame(id1 = rep(NA, nrow(dat2)),
id2 = rep(NA, nrow(dat2)),
correl = rep(NA, nrow(dat2)))
# initialize data frame to hold results
results <- as.data.frame(dat2[, -c(4:5)])
# initialize loop counters
x <- 1
# loop over the main ID's
for (i in ids) {
tmp <- dat2[ID1==i]
#loop over the ID's of the neighbour stations
for (id in 1:nrow(tmp)){
# get ID of the neighbours
near_id <- as.numeric(tmp[id, 2])
# get common dates
beg_date <- tmp[id, 4]
end_date <- tmp[id, 5]
# calculate correlations
correlations[x,1] <- cor(dat1[ID==i & date %between% c(beg_date, end_date)]$value,
dat1[ID==near_id & date %between% c(beg_date, end_date)]$value)
# Put correlation in df together with current id and near id
correlations2[x, "id1"] <- i
correlations2[x, "id2"] <- near_id
correlations2[x, "correl"] <- correlations[x,1]
# increment loop counter
x <- x + 1
}
}
# assemble final data frame
results <- data.table(ID=results[, 1],
ID_nearest=results[, 2],
distance=results[, 3],
overlapping_days=results[, 4],
correl=as.vector(correlations))
list(results, correlations2)
}
results_check <- cor_loop_check(dat1, dat2)
# Check results for e.g. row 20: Same value for correlation but differing id-pair ):
results_check[[1]][20,]
#> ID ID_nearest distance overlapping_days correl
#> 1: 1315 1551 1193032 2670 -0.06323207
results_check[[2]][20,]
#> id1 id2 correl
#> 20 1315 1559 -0.06323207
Создано в 2020-03-14 пакетом представ (v0.3.0)
Как видите. В строке 20 оба df содержат одинаковую корреляцию, но разные пары идентификаторов.