模式結果分析

比較模擬結果的方法有很多種，我們可以分為兩種類類型，(1) 模擬值與量測值的相關性、(2) 模擬與量測值的誤差。

基於相關性的模擬結果評估方法

基於相關性的評估方法中，最常見的方法是進行模擬與實測值的簡單線性回歸，通常會將模擬值輸出為y軸，實測值輸出為x軸，當斜率接近1、截距接近0、R²接近1時，我們可以可以認定為良好的模式。

另一個方法就是計算 r - 相關係數(Correlation coefficient)

$r=\tfrac{\sum\limits_{i=1}^n(o_i-\bar{o})(s_i-\bar{s})} {\sqrt{\sum\limits_{i=1}^n(o_i-\bar{o})^2}\sqrt{\sum\limits_{i=1}^n(s_i-\bar{s})^2}}$

符號 $o_i$ 代表實測值、 $s_i$ 代表相對應的模擬值

基於誤差的評估方法

常見的包括RMSE, RRMSE, MAE, EF，以下分別進行說明，

所使用的符號 $o_i$ 代表實測值、 $s_i$ 代表相對應的模擬值

均方根差 (RMSE)

最常見的方法為均方根誤差(root mean squre error,RMSE)，計算公式如下:

$RMSE = \sqrt{ \tfrac{1}{n} \times \sum\limits_{i=1}^n (o_i - s_i)^2 }$

另外也可以將RMSE除以實測值的平均值，成為相對均方根誤差(relative root mean square error, RRMSE)，公式如下

$RRMSE = \tfrac{RMSE} {\bar{o_i} }$ 

通常均方根誤差都帶有單位，例如RMSE = 653.2 (kg/ha)，代表模擬與實測的產量誤差有653.2，RRMSE就可以代表相對值例如0.23，可以提供我們了解RMSE 的的佔比。

無腦指標 - EF

推薦一個無腦的指標就是EF (Nash–Sutcliffe Efficiency)，公式如下 

 $EF=1-\tfrac{ \sum\limits_{i=1}^n (o_i-s_i)^2 } { \sum\limits_{i=1}^n (o_i - \bar{o_i})^2 }$

如果是完美的模式，模擬值會與實測值相同，也就是 $s_i = o_i$，因此EF = 1。 

如果每次的模擬值都是實測的平均值，也就是 $s_i = \bar{o_i}$，計算起來EF = 0，我們可以說EF = 0 的模式並不是好模式，只要用一個平均值去猜就跟模擬值一樣了；如果EF < 0，就顯示這個模式比平均值還要差。

在這邊可做一個簡單的結論是 

EF = 1，完美模式 

EF = 0，不好 (模式預測效果和使用平均值一樣) 

EF < 0，差 (模式預測結果比使用平均值猜測更糟糕) 

 EF其實也是R²的另一種形式

我們重新整理EF值的公式 ，可以改寫成如下

$EF = 1- \tfrac{MSE}{MSE_{\bar y}}$ 

其中，

$MSE_\bar{y} = \sum\limits_{i=1}^n (y_i- \bar{y}) = var(y)(n/1)/n$

可以理解$MSE_\bar{y}$為模擬值使用實測平均值($\bar o$) 所計算出來的平均誤差，EF值是比較MSE和$MSE_\bar{y}$ 的結果，因此EF的分子項(MSE)是用來計算模式無法解釋的觀察資料變異量，而分母則衡量的是觀察值的總變異量。

這正是決定係數(R²)的其中一種定義 

$R^2 = 1- \tfrac{TSS}{RSS}$ 

兩者都是 $ 1- \tfrac{模型預測誤差平方和}{觀察資料總變異}$

程式指令

R語言指令

     
Performance <- function(sim, obs){
  S <- sim
  O <- obs
  meanO <- mean(O)
  n <- length(S)
  wl = 0;   se = 0; dif = 0; sd = 0
  
  for(i in 1:n){
    ##  difference from sim and obs (ME)
    dif <- dif + (O[i]-S[i])
    ## square value of difference between sim and obs - se (RMSE,ME,EF,d)
    se <- se + (O[i]-S[i])^2
    ## calculate deviation square for d
    dist <- abs(O[i]-meanO)+abs(S[i]-meanO)
    wl<- wl + dist^2
    ## calculate the squared distance value of observed (EF)
    sd <- sd + (O[i]-meanO)^2
  }
  ## calculate the parameter 
  ME <- dif/n
  RMSE <- sqrt(se/n)
  d <- 1 - se/wl
  EF <- 1- se/sd
  ## linear regression
  lmout <- lm(S~O)
  intercept <- summary(lmout)$coefficient[1]#intercept
  slope <- summary(lmout)$coefficient[2]#slope
  r2 <- summary(lmout)$r.squared
  p <- summary(lmout)$coefficients[8]#p value
  
  return(c(n,round(RMSE,3), round(ME,3), round(d,3), round(EF,3), 
           round(intercept,3), round(slope,3),round(r2,3),round(p,3))
         )
}

參考文獻

1. Soltani, A., T.R. Sinclair (2012) Modeling Physiology of Crop Development, Growth and Yield. CABI, USA. 
2. Wallach, D., D. Makowski, J.W. Jones, and F. Brun (2019) Working with Dynamic Crop Models: Methods, Tools, and Examples for Agriculture and Environment. Elsvier, United Kingdom