Построение графика совокупной неопределенности по совокупным средним значениям данных временных рядов в 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)) %>%  # removing Date and Time columns: cannot be averaged
        group_by(hourlyDate) %>% 
        summarize(
                across(everything(),mean))

    ## make cumulative data 
    
    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" "2024-06-11 02:00:00" "2024-06-11 03: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

Полученная цифра такая:

Я могу написать код для расчета стандартного отклонения, используя summarise() в dplyr, но не могу распространить его на кумулятивные средства. Мне трудно найти пакет R или функцию tidyverse для разработки такого графика (https://doi.org/10.1038/s41561-021-00785-2):

Заранее спасибо.

РЕДАКТИРОВАТЬ.

Я пытался улучшить код, но почему-то geom_ribbon() не может отобразить неопределенность.

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