在 POSIXt 数据上拟合正弦波模型并使用 Ggplot2 绘图

Question

长期读者，第一次在这里提问:)

我在特定时间和日期收集了一些数据，有理由假设这些数据大致遵循 24 小时周期。我想在我的数据上拟合一个正弦波模型作为时间的函数，这样就可以测试未来的数据点是否落在预测的模式上。

我已阅读 this、this 和 this 的回复，但它们并没有解决我的问题，因为就我而言，我希望将 x 轴数据保留为 POSIXct 日期时间格式。这就是收集数据的方式，使用这种格式可以轻松解释图表。

这里有一些与我的真实数据相同的可重现数据：

time <- c("2022-01-01 09:20:00", "2022-01-02 11:10:00", 
          "2022-01-02 18:37:00", "2022-01-03 14:01:00", 
          "2022-01-05 06:50:00", "2022-01-06 17:03:00")

time <- as.POSIXct(time)

value <- c(3, 6, 2, 8, 4, 1)

这些在基础 R 中绘制得很好：

plot(time, value)

但是，当我尝试构建适合时间序列的正弦回归模型时，我遇到了麻烦。我也在努力完全理解 nls 函数所需的参数。根据前面的例子，我已经尝试过这种方法（用 cmets 了解它的工作原理）：

res <- nls(value ~ A * sin(omega * time + phi) + C,        # This is the basic sine-function format
           data = data.frame(time, value),                 # This defines the data used
           start = list(A = 1, omega = 1, phi = 1, C = 1)) # This gives nls the starting values?

在这里，我收到一条错误消息：“Ops.POSIXt(omega, time) 中的错误：'*' 未为“POSIXt”对象定义”我将其解释为我想使用的特定日期格式不是这种方法可以接受。我知道这一点，因为如果我只是用整数的虚拟向量替换时间变量，模型就可以正常工作，我可以将其绘制如下：

time2 <- c(1, 2, 3, 4, 5, 6)
res <- nls(value ~ A * sin(omega * time2 + phi) + C,
       data = data.frame(time, value),
       start=list(A=1, omega=1, phi=1, C=1))

coefs <- coef(res)

fit <- function(x, a, b, c, d) {a * sin(b * x + c) + d}

plot(time2, value)
curve(fit(x, a = coefs["A"], b = coefs["omega"], 
             c = coefs["phi"], d = coefs["C"]), add=TRUE, 
             lwd=2, col="red")

我知道我走在正确的轨道上，但我的主要问题是，如何在保持 POSIXct 格式的时间变量的同时完成上述过程？

如前所述，我的主要业务顺序是使用 Ggplot2 绘制数据，但在解决最初的问题之前我什至无法开始尝试。但是，非常感谢任何有关如何开始使用的指示！ :)

Answer 1

我可能只是从任意起始时间生成一个数字天数并使用它。然后，您可以修改 fit 函数，以便将日期时间转换为预测值。然后，您可以轻松地从您的模型中创建一个预测数据框并绘制它。

df <- data.frame(time = time, value = value)

origin <- as.POSIXct("2022-01-01 00:00:00")

df$days <- as.numeric(difftime(time, origin, unit = "day"))

res <- nls(value ~ A * sin(omega * days + phi) + C,  
           data = df, 
           start = list(A = 1, omega = 1, phi = 1, C = 1))

fit <- function(res, newdata) {
  
  x <- as.numeric(difftime(origin, newdata$time, units = "days"))
  C <- as.list(coef(res))
  C$A * sin(C$omega * x + C$phi) + C$C
}

new_df <- data.frame(time = origin + as.difftime(new_times, units = "days"))
new_df$value <- fit(res, new_df)

ggplot(df, aes(time, value)) +
  geom_point() +
  geom_line(data = new_df, colour = "gray") +
  theme_bw()