如何使用facet_wrap将NSE和PBIAS结果添加到ggplot中？

Question

我是一个初学者，我创建了一个函数（如下所示）来计算模拟数据与观察数据的百分比偏差（PBIAS）和纳什-舒特克里夫效率（NSE）。但是，我只能使用整个数据集来计算这些测试。

model.assess <- function(Sim, Obs) { 
  rmse = sqrt( mean( (Sim - Obs)^2, na.rm = TRUE) ) #Formula to calculate RMSE
  RSR <- rmse / sd(Obs) #object producing RSR test from the RMSE formula
  PBIAS <- 100 *(sum((Sim - Obs)/sum(Obs), na.rm =TRUE)) #object producing PBIAS test
  NSE <- 1 - sum((Obs - Sim)^2)/sum((Obs - mean(Obs))^2, na.rm =TRUE) #object producing NSE test
  stats <- print(paste0("RSR = ", sprintf("%.3f", round(RSR, digits=3)), "    PBIAS = ", sprintf("%.3f",round(PBIAS, digits=3)),"    NSE = ", sprintf("%.3f",round(NSE, digits=3))))  
  return(stats) #returns the results of the tests with 3 decimals and spacing in between

这是我的数据集，四个不同电台（SNS，MRC，TLG，SJF）的每月流量：

StationID Date      Obs_flow    Sim_flow      Month     Year
SNS    1950-10-01   0.010170    0.030687967 October 1950-01-01      
SNS    1950-11-01   0.366260    0.416466741 November 1950-01-01     
SNS    1950-12-01   0.412210    0.496136731 December 1950-01-01     
SNS    1951-01-01   0.119520    0.182072570 January 1951-01-01      
SNS    1951-02-01   0.113480    0.142611192 February 1951-01-01     
SNS    1951-03-01   0.127090    0.176350274 March   1951-01-01  
SNS    1951-04-01   0.175120    0.193221389 April   1951-01-01      
SNS    1951-05-01   0.208940    0.275980903 May     1951-01-01  
SNS    1951-06-01   0.114420    0.144675317 June    1951-01-01      
SNS    1951-07-01   0.032280    0.018057796 July    1951-01-01  

用等式和R平方绘制Obs vs Sim的散点图：

dataset %>%
  filter(StationID == "SNS") %>%
ggplot(aes(x = Obs_flow, y = Sim_flow)) + 
   geom_point(aes(Obs_flow, Sim_flow), alpha = 0.3)+
     stat_smooth(aes(x = Obs_flow, y = Sim_flow), 
                method = "lm", se = TRUE, colour="#FC4E07", fullrange = TRUE) + 
   stat_poly_eq(formula = "y~x", 
             aes(label = paste0(..eq.label..)),  #adding the equation on the top
             parse = TRUE, label.x.npc = "center", label.y.npc = 0.97, size = 3.45, family= "Times New Roman")+
     stat_poly_eq(formula = "y~x", 
             aes(label = paste0(..rr.label..)), #adding the Rsquared at the bottom
             parse = TRUE, label.x.npc = 0.95, label.y.npc = 0.05, size = 3.45, family= "Times New Roman")+

  annotate("text", x = 0, y = 1.3,, label = paste0(model.assess(dataset$Sim_flow, dataset$Obs_flow)),  collapse = "\n", hjust = 0, size=2.4, family= "Times New Roman") +

   facet_wrap(~ Month, ncol=4,  labeller = labeller(StationID = c("MRC" = "Merced River", "SJF"= "Upper San Joaquin River", "SNS" = "Stanislaus River", "TLG" = "Tuolumne River")), scales = "fixed")```

"stat_poly_eq" added an equation and Rsquared for each facet, but the annotate adds the same number for all facets. Is there a way to add NSE and PBIAS for each facet separately? I tried the package HydroGOF, but I got the same result. Excuse the aesthetics.

  [1]: https://i.stack.imgur.com/rDZ4z.png

Answer 1

样本数据将有助于其他人为您提供帮助。请查看此link以供将来查询。

您有一些问题。您的model.assess()函数提供一条记录，而您需要每个方面的值。因此，我使用代码

创建了一个虚拟对象

ll <- data.frame(Month=c(),label=c())
nM <- length(Month)
lapply(1:nM, function(i){
  a <- Sim_flow*i*i*0.5
  b <- Obs_flow*i
  m <- model.assess(a,b)
  ll <<- rbind(ll,data.frame(Month=Month[i],label=m))
})
labels <- ll

接下来，您需要使用geom_label代替here进行注释。下面的代码

ggplot(data=dataset, aes(x = Obs_flow, y = Sim_flow)) + 
  geom_point(aes(Obs_flow, Sim_flow), alpha = 0.3)+
  stat_smooth(aes(x = Obs_flow, y = Sim_flow), 
              method = "lm", se = TRUE, colour="#FC4E07", fullrange = TRUE) + 
  stat_poly_eq(formula = "y~x", 
               aes(label = paste0(..eq.label..)),  #adding the equation on the top
               parse = TRUE, label.x.npc = "center", label.y.npc = 0.97, size = 3.45, family= "Times New Roman") +
  stat_poly_eq(formula = "y~x", 
               aes(label = paste0(..rr.label..)), #adding the Rsquared at the bottom
               parse = TRUE, label.x.npc = 0.95, label.y.npc = 0.05, size = 3.45, family= "Times New Roman") +
  facet_wrap(~Month, ncol=4,  labeller = labeller(StationID = c("MRC" = "Merced River", "SJF"= "Upper San Joaquin River", "SNS" = "Stanislaus River", "TLG" = "Tuolumne River")), scales = "fixed") +
  geom_label(data = labels, aes(label=label, x = Inf, y = -Inf), 
             hjust=1, vjust=0, size=1.8,
             inherit.aes = FALSE)

给出以下