在R中绘制时间序列累积均值的累积不确定性图的技术问题

我在处理30分钟间隔的温室气体通量时间序列数据时，遇到了绘制带累积不确定性的累积图的难题。已通过cumsum()将原始数据聚合为小时累积数据，但无法将标准差传播到累积均值中，也没法画出类似文献中带不确定性丝带的累积图。

原始聚合与绘图代码 #

### ---- Binning Net emission to times. Averaging half-hourly to Hourly Data 
str(gp_NEE  )
# 输出：
# tibble [1,270 × 18] (S3: tbl_df/tbl/data.frame)
#  $ DateTime: POSIXct[1:1270], format: "2024-06-11 00:00:00" "2024-06-11 00:30:00" ...
#  $ daytime : num [1:1270] 0 0 0 0 0 0 0 0 0 0 ...
#  $ NEE     : num [1:1270] 1.831 0.702 2.017 8.33 3.872 ...

gp_NEE_hour <- gp_NEE %>% 
        mutate(hourlyDate = floor_date(DateTime,unit='hour')) %>%
        select(-c(DateTime)) %>%  # 移除无法平均的日期时间列
        group_by(hourlyDate) %>% 
        summarize(
                across(everything(),mean))

## 生成累积数据

gp_NEE_hour<- gp_NEE_hour %>%
        mutate(
            cum_NEE = cumsum(replace_na(NEE, 0)))

str(gp_NEE_hour)
# 输出：
# tibble [635 × 19] (S3: tbl_df/tbl/data.frame)
#  $ hourlyDate: POSIXct[1:635], format: "2024-06-11 00:00:00" "2024-06-11 01:00:00" ...
#  $ daytime   : num [1:635] 0 0 0 0 0 0 1 1 1 1 ...
#  $ NEE       : num [1:635] 1.27 5.17 2.58 NA 1.62 ...
#  $ cum_NEE   : num [1:635] 1.27 6.44 9.02 9.02 10.64 ...


cumplot <- gp_NEE_hour %>%
            ggplot( aes(x=hourlyDate, y=cum_NEE)) +
                    geom_line( color="grey") +
                    geom_point(shape=21, color="black", fill="#69b3a2", size=1) +
                    theme_bw() +theme(legend.position="none", 
        text = element_text(size=15))+
    labs(y=expression(CO[2]~Flux~Netto~Kumulatif~~"("~g~m^{-2}~~jam^{-1}~")"))+xlab ( "Waktu")
cumplot

这段代码生成了基础的累积通量曲线，但没有不确定性区间。

我能用dplyr::summarise()计算每小时数据的标准差，但不知道如何将不确定性传播到累积均值中，也找不到合适的R包或tidyverse函数绘制类似文献中带阴影不确定性区间的累积曲线。

尝试添加不确定性的代码 #

gp_NEE_hour_no_NA <-gp_NEE %>% 
                    mutate(NEE_Mgha = NEE * 0.24) %>% 
                    select(DateTime, NEE_Mgha) %>% 
                    mutate(hourlyDate = floor_date(DateTime,unit='hour')) %>%
                    group_by(hourlyDate) %>% 
                    summarise(
                            meanH_NEE = mean(NEE_Mgha), 
                            sd_NEE = sd(NEE_Mgha), 
                            n = n(),
                            se_NEE = sd_NEE / sqrt(n),.groups = "drop") %>%
                    mutate(cum_NEE = cumsum(replace_na(meanH_NEE, 0)))%>% 
                    mutate(cum_se_NEE = sqrt(cumsum(se_NEE^2)))
                    
str(gp_NEE_hour_no_NA)
# 输出：
# tibble [635 × 7] (S3: tbl_df/tbl/data.frame)
#  $ hourlyDate: POSIXct[1:635], format:  ...
#  $ meanH_NEE : num [1:635] 0.304 1.242 0.619 NA 0.39 ...
#  $ sd_NEE    : num [1:635] 0.1915 1.0714 0.4392 NA 0.0646 ...
#  $ n         : int [1:635] 2 2 2 2 2 2 2 2 2 2 ...
#  $ se_NEE    : num [1:635] 0.1354 0.7576 0.3106 NA 0.0457 ...
#  $ cum_NEE   : num [1:635] 0.304 1.546 2.164 2.164 2.554 ...
#  $ cum_se_NEE: num [1:635] 0.135 0.77 0.83 NA NA ...

plot_NEE2 <-  ggplot(data = gp_NEE_hour_no_NA, aes(x = hourlyDate, y = cum_NEE)) +                  
                geom_line() + 
                geom_ribbon(aes(ymin = cum_NEE - 1.96 * cum_se_NEE, ymax = cum_NEE + 1.96 * cum_se_NEE), alpha = .5,
                            fill = "darkseagreen3", color = "transparent") +                    
                #geom_point(aes(y = mean),color="aquamarine4")+                
                ylim(-1.5, 1.5) + theme(text = element_text(size=15))+
                labs(y=expression(CO[2]~Flux~Netto~Kumulatif~~"("~Mg~ha^{-1}~~hari^{-1}~")"))+xlab("Waktu") + 
                theme_bw()
plot_NEE2

这段代码中geom_ribbon无法显示不确定性区间，特此寻求解决方法。

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