按小时积分theMass结果异常：线性增长问题排查及实现方案问询

Hey George, let's work through this problem step by step—that weird linear growth you're seeing is definitely a clue that something's off with either your grouping logic, averaging step, or integration calculation. Let's break this down and fix it:

The most common cause of linear growth here is incorrect time grouping. If you're only grouping by the hour of day (0-23) instead of date + hour, you're aggregating data across multiple days into a single hourly bucket. That would make each "hourly" result the sum of every same-hour slot across all days, leading to a steady upward trend.

Fix this by using lubridate to create a proper hourly time window that includes the date:

library(lubridate)
library(dplyr)

# Make sure your date_time is a proper POSIXct time format first
df$date_time <- as.POSIXct(df$date_time, format = "%Y-%m-%d %H:%M:%S") # Adjust format to match your data

# Create a group for each full hour (date + hour)
df$hour_window <- floor_date(df$date_time, unit = "hour")

You mentioned vs and t_in need hourly averages to compute theMass accurately. Make sure you're calculating these averages per hourly window, not globally or with rolling averages. Use dplyr grouping to keep this clean:

# Get hourly averages for vs and t_in
hourly_metrics <- df %>%
  group_by(hour_window) %>%
  summarise(
    avg_vs = mean(vs, na.rm = TRUE),
    avg_t_in = mean(t_in, na.rm = TRUE)
  )

# Merge these averages back into your original minute-level data
# This lets you use the hourly average for every minute in that window
df_with_hourly_avg <- df %>%
  left_join(hourly_metrics, by = "hour_window")

Now recalculate theMass using the hourly averages, then integrate over each hour. Since your data is minute-level, each interval is 1/60 of an hour—you can use simple summation (for discrete values) or the trapezoidal rule (for smoother integration):

Option 1: Simple Summation (Discrete Minute Values) #

hourly_integral <- df_with_hourly_avg %>%
  # Replace this with your actual theMass formula using avg_vs and avg_t_in
  mutate(calculated_mass = avg_vs * avg_t_in) %>%
  group_by(hour_window) %>%
  summarise(
    total_mass = sum(calculated_mass * (1/60), na.rm = TRUE)
    # Multiply by 1/60 to convert minute-level values to hourly total
  )

Option 2: Trapezoidal Rule (Smoother Integration) #

If you want a more precise integral, use the trapezoidal method with the pracma package:

library(pracma)

hourly_integral <- df_with_hourly_avg %>%
  mutate(calculated_mass = avg_vs * avg_t_in) %>% # Your theMass formula
  group_by(hour_window) %>%
  summarise(
    total_mass = trapz(date_time, calculated_mass) / 3600
    # Divide by 3600 to convert seconds (from POSIXct) to hours
  )

If you still see linear trends after fixing grouping, check these common issues:

• Accidental cumulative sum: Make sure you're using sum() per group, not cumsum() (which would accumulate totals across hours).
• Time format errors: If date_time was stored as a string instead of POSIXct, grouping might sort incorrectly and mess up calculations.
• Global averages instead of hourly: Double-check that you're using avg_vs/avg_t_in from the hourly group, not a global mean of the entire dataset.

Here's a complete, reproducible example to test with:

library(lubridate)
library(dplyr)

# Simulate minute-level test data
set.seed(123)
df <- tibble(
  date_time = seq(as.POSIXct("2024-01-01 00:00:00"), as.POSIXct("2024-01-02 23:59:00"), by = "1 min"),
  vs = rnorm(nrow(df), mean = 10, sd = 2),
  t_in = rnorm(nrow(df), mean = 5, sd = 1)
)

# Create hourly windows
df$hour_window <- floor_date(df$date_time, "hour")

# Get hourly averages
hourly_metrics <- df %>%
  group_by(hour_window) %>%
  summarise(avg_vs = mean(vs), avg_t_in = mean(t_in))

# Merge and calculate theMass
df_with_avg <- df %>%
  left_join(hourly_metrics, by = "hour_window") %>%
  mutate(calculated_mass = avg_vs * avg_t_in) # Replace with your formula

# Compute hourly integral
hourly_integral <- df_with_avg %>%
  group_by(hour_window) %>%
  summarise(total_mass = sum(calculated_mass * (1/60)))

# View results (no linear growth here!)
print(hourly_integral)

内容的提问来源于stack exchange，提问作者George