r - 如何执行不等频率的时间序列之间的相关性

Question

我每分钟测量室温 36 分钟，同时每秒测量皮肤温度 32 次。我有 35 次重复的实验标记为 (ID)。我需要能够查看相关性，但是样本大小不等。

数据:

我有一个每分钟测量一次室温的 data.frame df1 和另一个每秒测量 32 次皮肤温度的 data.frame df2。我有 36 分钟的数据。此外，还有另一列名为 ID 的列显示了实验编号 (1-35)，但我不知道如何在以下示例数据中表示这一点。所以从技术上讲，我正在寻找基于 ID 的每个 SkinTemp 与 RoomTemp 的相关性。

    df1 <- data.frame(
        roomTemp = rnorm(1*36),
    )

   df2 <- data.frame(
        skinTemp = rnorm(32*60*36),
        )

我试着做:

Data <- data.frame(
  Y=c(df1,df2),
  Variable =factor(rep(c("RoomTemp", "SkinTemp"), times=c(length(df1), length(df2))))
)

cor(Data$Y~Data$Variable)

但这似乎不起作用。

Answer 1

滚动连接或插值可能有助于在测量 roomTemp 时输入 skinTemp。下面是两者的例子。第一部分是处理多个ID的更新，后面是针对单个ID的情况的原始答案。

更新:处理多个 ID 的新版本

此更新解决了具有多个 ID 的数据的情况，我们希望为每个 ID 分别进行插值或滚动连接。

library(data.table)
library(reshape2)
library(dplyr)
library(purrr)
library(ggplot2)
theme_set(theme_classic(base_size=16))

首先，我们将为两个独立的 ID 创建虚假的自相关数据:

set.seed(395)
df1 <- data.frame(roomTemp = c(cumsum(rnorm(1*36)), cumsum(rnorm(1*36))),
                  ID = rep(c("A","B"), each=36))
df2 <- data.frame(skinTemp = c(cumsum(rnorm(32*60*36,0,0.01)),
                               cumsum(rnorm(32*60*36,0,0.01))),
                  ID = rep(c("A","B"), each=32*60*36))

现在我们添加一个时间列，但在这种情况下，我还在 df1 中添加了一个偏移，以便没有 df1 测量与 df2 测量同时发生，只是为了使答案更通用。

# Add time column
df1$time = rep(0:(0.5*nrow(df1)-1)*60 + 0.0438,2)
df2$time = rep(0:(0.5*nrow(df2)-1)/32, 2)

将数据框转换为数据表。这一次，除了 ID 之外，我们还使 time 成为关键列，以便滚动连接将针对每个 ID 单独发生。

# Convert data frames to data tables
setDT(df1)
setDT(df2)

# Make ID and time key columns in both data frames (for joining)
setkey(df1, ID, time)
setkey(df2, ID, time)

# Rolling join roomTemp to nearest time value of skinTemp
df2 = df1[df2, roll="nearest"]

# Rename rolling joined room temperature column
names(df2)[grep("roomTemp", names(df2))] = "roomTempRoll"

要通过 roomTemp 添加插入的 ID ，我使用了 map_df 包中的 purrr 。 map_df 对每个 ID 单独操作。 approx 负责插值。在最初的答案中，我首先使用 approxfun 创建了一个近似函数，但在这里我只是直接在一个步骤中完成了插值。 map_df 返回一个数据帧，但我们只需要 y 列，它具有 roomTemp 的内插值，所以我在 dplyr 函数链的末尾提取了这些值，并将它们分配给 roomTempInterp 中的 df2 。

# Add interpolated room temperature by ID
df2$roomTempInterp = unique(df2$ID) %>% 
  map_df(~ approx(df1$time[df1$ID==.x], df1$roomTemp[df1$ID==.x], 
                  xout=df2$time[df2$ID==.x]), .id="ID") %>% .$y

在下图中，我们通过 ID 分面，以便我们可以分别查看每个 ID 的估算温度值。

# Plot so we can see what the rolling joined room temperature and 
#  interpolated room temperature look like
ggplot(melt(df2, id.var=c("ID", "time")), aes(time, value, colour=variable)) +
  geom_line(size=0.7) +
  geom_point(data=df1, aes(time, roomTemp), colour="black") +
  facet_grid(ID ~ .)

这是通过 ID 获取相关性的一种方法:

df2 %>% group_by(ID) %>%
  summarise(r_interp = cor(skinTemp, roomTempInterp, use="pairwise.complete.obs"),
            r_roll = cor(skinTemp, roomTempRoll, use="pairwise.complete.obs"))

      ID    r_interp      r_roll
1      A -0.04853998 -0.02993207
2      B -0.53993960 -0.53092150

原答案

首先，我修改了示例数据框以添加一些自相关，因为这看起来更接近您的真实实验并且使可视化更容易。

library(data.table)
library(reshape2)
library(dplyr)
library(ggplot2)
theme_set(theme_classic(base_size=16))

# Fake data with autocorrelation
set.seed(395)
df1 <- data.frame(roomTemp = cumsum(rnorm(1*36)))
df2 <- data.frame(skinTemp = cumsum(rnorm(32*60*36,0,0.01)))

现在添加一个时间列。您可以使用实际的日期时间列，但在这里我刚刚使用了以秒为单位的数字列。

# Add time column
df1$time = 0:(nrow(df1)-1)*60
df2$time = 0:(nrow(df2)-1)/32

对于插值，我们需要一个函数来在室温测量之间测量皮肤温度时对室温进行插值。 approxfun 在点之间执行线性插值。您还可以以类似的方式使用 splinefun 来使用样条进行插值。

# Function to interpolate room temperature between measurements
roomTempInterp = approxfun(df1$time, df1$roomTemp)

将数据框转换为数据表以使用 data.table 的滚动连接功能。

# Convert data frames to data tables
setDT(df1)
setDT(df2)

# Make time a key column in both data frames (for joining)
setkey(df1, time)
setkey(df2, time)

现在对最近的时间值执行滚动连接。

# Rolling join roomTemp to nearest time value of skinTemp
df2 = df1[df2, roll="nearest"]

# Rename rolling joined room temperature column
names(df2)[grep("roomTemp", names(df2))] = "roomTempRoll"

将原始 roomTemp 测量值从 df1 合并到 df2 中。

df2 = df1[df2, ]  # Equivalent to dplyr: df2 = left_join(df2, df1)

使用我们上面创建的函数添加内插的室温。

# Add interpolated room temperature
df2$roomTempInterp = roomTempInterp(df2$time)

插值方法对我来说似乎更现实，特别是如果我们可以假设 roomTemp 在测量之间相对平滑和单调地变化。下面是 df2 的前 10 行，其中包括原始 df2 数据加上新的 roomTempRoll 和 roomTempInterp 列以及来自 roomTemp 的原始 df1 测量值。您现在可以使用此数据框来评估 roomTemp 和 skinTemp 之间的相关性和其他关系。

    roomTemp    time roomTempRoll     skinTemp roomTempInterp
 1: -1.21529 0.00000     -1.21529 -0.006511475      -1.215290
 2:       NA 0.03125     -1.21529 -0.014058076      -1.215531
 3:       NA 0.06250     -1.21529 -0.017741690      -1.215773
 4:       NA 0.09375     -1.21529 -0.030211177      -1.216014
 5:       NA 0.12500     -1.21529 -0.027105225      -1.216255
 6:       NA 0.15625     -1.21529 -0.035784295      -1.216497
 7:       NA 0.18750     -1.21529 -0.031319748      -1.216738
 8:       NA 0.21875     -1.21529 -0.033758959      -1.216979
 9:       NA 0.25000     -1.21529 -0.040667384      -1.217220
10:       NA 0.28125     -1.21529 -0.026291442      -1.217462

下面是一个图，您可以看到滚动连接和插值的样子。黑点标记原始 roomTemp 测量值。

ggplot(melt(df2 %>% select(-roomTemp), id.var="time"), aes(time, value, colour=variable)) +
  geom_line(size=1) +
  geom_point(data=df2, aes(time, roomTemp), colour="black")