论文中常见的拟合散点验证图（R语言版）

如上图所示，是论文中常见的validation图，python也能实现相似的图绘。

今天先介绍R语言版，python改期再介绍吧

这张图需要依次实现下列功能：

data实测与data模拟的散点绘制
R方、RMSE、Bias等参量计算与绘制
拟合后的直线与y=x直线
网格线（如果需要）

使用base R的功能一步一步计算：

为了实现多个因变量和一个自变量在同一个图片里，我们要用points或者lines函数画其他因变量和自变量的值几点注意：

画图的时候，明明设置了X轴Y轴的范围，为什么坐标轴的长度还是会多出一块呢？看下图：

这里用到了xaxs和yaxs，可以使x，y轴从0开始

直接上代码

# 数据读取
library(readxl)
library(stringr)
library(data.table)
DE_Hai <- read_excel("D:/acad3_211001/data/DE-Hai.xlsx", 
    sheet = "DE-Hai")
Qsim <- DE_Hai$Gc; Qobs <- DE_Hai$Gc1

Qsim为模拟的数据，Qobs为观测的数据，自行读取

再手动写一下求RMSE等参数的函数：

{r}
NSE <- function(Qsim, Qobs){
  
  # exclude NA
  data <- data.table(Qsim, Qobs)
  data <- na.omit(data)
  Qsim <- data$Qsim
  Qobs <- data$Qobs
  
  # compute
  Qobs_mean <- mean(Qobs)
  Emod <- sum((Qsim - Qobs)^2)
  Eref <- sum((Qobs - Qobs_mean)^2)
  if (Emod == 0 & Eref == 0) {
    NSE <- 0
  } else {
    NSE <- (1 - Emod / Eref)
  }
  NSE
}
RMSE <- function(Qsim, Qobs){
  # exclude NA
  data <- data.table(Qsim, Qobs)
  data <- na.omit(data)
  Qsim <- data$Qsim
  Qobs <- data$Qobs
  
  sqrt(mean((Qsim - Qobs)^2))
}
Bias <- function(Qsim, Qobs){
  
  # exclude NA
  data <- data.table(Qsim, Qobs)
  data <- na.omit(data)
  Qsim <- data$Qsim
  Qobs <- data$Qobs
  
  # BIAS
  sumQobs <- sum(Qobs)
  sumQsim <- sum(Qsim)
  if (sumQobs == 0) {
    BIAS <- 0
  } else {
    BIAS <- (sumQsim - sumQobs) / sumQobs
  }
  BIAS
}

接下来是核心部分：

# line plot
plot(x=Qsim, y=Qobs,
     col=rgb(0.4,0.4,0.8,0.6),pch=16, cex=1.3,
     xlab="Gc by PML", ylab="Gc by LE",
     xlim=c(0,0.01) , ylim=c(0,0.01), yaxs="i",xaxs="i")
grid(nx=5,ny=5,lwd=1,lty=2,col=rgb(0.4,0.4,0.8,0.2))

line.model <- lm(Qobs~Qsim)
k <- as.numeric(line.model$coefficients[2])
b <- as.numeric(line.model$coefficients[1]) 
p <- (summary(line.model))$coef[2,4]

lines(seq(-0.001,0.011,length=100), 
      k * seq(-0.001,0.011,length=100) + b, 
      col='#c27ba0', lwd=2)

lines(seq(-0.001,0.011,length=100), 
      seq(-0.001,0.011,length=100), col="#3c78d8", lwd=2)
text(0.007, 0.002, 
     str_sub(paste('NSE=',NSE(Qsim, Qobs), sep=' '),1, 10))
text(0.009, 0.002, str_sub(paste('R^2=',
  R2 = cor(Qsim, Qobs, use = "complete")^2, sep=' '),1, 9))
text(0.007, 0.001, str_sub(paste('RMSE=',
  RMSE(Qsim, Qobs),sep=' '),1, 12))
text(0.009, 0.001, str_sub(paste('bias=',
  Bias(Qsim, Qobs), sep=' '),1, 12))
text(0.007, 0.004, str_sub(paste('k=', k, sep=' '), 1, 8))
text(0.009, 0.004, 'p < 0.05')