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

Question

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

数据：

我有一个 data.frame df1 每分钟测量一次室温，另一个 data.frame df2 每秒测量 32 次皮肤温度。我有 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

滚动连接或插值可能有助于在测量 skinTemp 时估算 roomTemp。以下是两者的示例。第一部分是处理多个 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)

将数据框转换为数据表。这一次，除了time，我们将ID 设置为一个键列，这样每个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"

要添加由ID 插入的roomTemp，我使用了purrr 包中的map_df。 map_df 在每个 ID 上单独运行。 approx 负责插值。在原始答案中，我首先使用approxfun 创建了一个近似函数，但在这里我只是一步直接完成了插值。 map_df 返回一个数据框，但我们只需要 y 列，它具有 roomTemp 的插值，所以我提取了 dplyr 函数链末尾的那些并将它们分配给 @987654349 @在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"

将来自df1 的原始roomTemp 测量值合并到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 列以及来自 df1 的原始 roomTemp 测量值。您现在可以使用此数据框来评估 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")

Answer 2

下面，我提供了一个关于如何实现这种关联的最小示例。

您可以在下面查看我的 cmets，但实际上我所做的是为每个室温观察时间创建箱（或“桶”）。然后，我将皮肤温度观测值（大大超过室温观测值）按相应的箱子汇总。因此，由于每 36*60*32 个皮肤 温度观察得到一个 房间 温度观察，因此前 36*60*32 皮肤温度观察被滚动到 bin '1 '。该过程从那里继续，来自 [36*60*32, 36*60*32*2] 的皮肤温度观测值被卷入 bin“2”，依此类推。

library(lubridate)
library(dplyr)

# create the times of our observations
time.room.temp <- seq.POSIXt(from = as.POSIXct('02/20/2017', format = '%m/%d/%Y'), to = as.POSIXct('02/21/2017', format = '%m/%d/%Y'), by = 36*60)
time.skin.temp <- seq.POSIXt(from = as.POSIXct('02/20/2017', format = '%m/%d/%Y'), to = as.POSIXct('02/21/2017', format = '%m/%d/%Y'), by = 1/32)

n.obs.room.temp <- length(room.temp)
n.obs.skin.temp <- length(skin.temp)

# create some "actual" temperature data
obs.room.temp <- rnorm(n.obs.room.temp, mean = 60, sd = 10)
obs.skin.temp <- rnorm(n.obs.skin.temp, mean = 95, sd = 5)

room.temp.df <- data.frame('room temp' = obs.room.temp, 'time' = time.room.temp)
skin.temp.df <- data.frame('skin temp' = obs.skin.temp, 'time' = time.skin.temp)

# Every 32 indices, seconds is incremented by one.. So our modulus calculuation should be every
# time the index evenly divides 36*60*32... there are 69120 skin-temp observations for every room-temp observation
# So we can effectively "bin" the different seconds for which we observed skin temperatures in order to create a mean temperature by bin,
# i.e. a mean skin temperature for every time at which room temp was recorded
bins <- cut(1:n.obs.skin.temp, seq(0, n.obs.skin.temp, 36*60*32), labels = 1:40)
skin.temp.df$bins <- bins

# Now, we can effectively group skin temperature observations by room temperature observations, and get the average (or median, if you like)
# temperature for each bin
shorter.skin.temp.df <- skin.temp.df %>%
  group_by(bins) %>%
  summarise(average.skin.temp = mean(skin.temp))

# Now we can get the correlation between the two types of temperatures!
cor(room.temp.df$room.temp, shorter.skin.temp.df$average.skin.temp)

编辑：对独特室温和“卷起”的皮肤温度观察数量的一点点验证：

> print(length(unique(skin.temp.df$bins)))
[1] 41
> print(length(unique(room.temp.df$time)))
[1] 41

因此，您可以高枕无忧，因为每个独特的室温观察时间都有一个相应的独特的皮肤温度观察时间箱。