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

论文中常见的拟合散点验证图(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')
  • grid用于添加格网线,更美观
  • lm用于线性拟合
  • lines函数是在已有图绘上增加内容
  • paste和str_sub用于连接和截取字符串,来自于stringr包

结果如图:

  • 1
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

地学万事屋

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值