consider adding unbiased root mean square difference to the metrics

[DanCodigaNOAA]

I did not see "unbiased root mean square difference" (URMSD), also known as "centered RMSD", in the list of metrics and if it's not there I suggest that it be included.

Using O for observation and M for model,
URMSD = stdev(M-O) = $< (M - < M > )^2 + (O - < O > )^2 >$

A reason to include it is that the following six metrics are essentially a "complete" set:

RMSD
URMSD
Mean Bias (MB)
stdev(M)
stdev(O)
CC

where the first 5 all have the same units (say, meters, for sea level) and the last one (CC) is the familiar unitless Pearson's correlation coefficient:

CC = $< (M - < M >) * (O - < O >) > / [stdev(M) * stdev(O)]$ .

These six are all part of the fundamental relationships:

$RMSD^2 = URMSD^2 + MB^2$

$URMSD^2 = stdev(M)^2 + stdev(O)^2 - 2* CC* stdev(M)*stdev(O)$

$RMSD^2 = stdev(M)^2 + stdev(O)^2 - 2* CC*stdev(M)*stdev(O) + MB^2$

[DanCodigaNOAA]

Just to continue this thought... when comparing model output to observations one could use all 6 of those parameters.

Alternatively, there is a set of 4 parameters which can get you a similarly ~complete set of information:

RMSD, MB, CC, and VD

where VD is the variance dissimilarity (unitless -1 to +1), as explained at #20 (comment)

It can be more convenient to use 4 metrics instead of 6.

[tomsail]

Hi Dan,

thanks a lot for this input.
I think the URMSD (or centered RMS) is the one actually used in the taylor diagrams (we use them for model skill assessment).

The rms metric implemented is already by default this precise metric.

And actually all metric you mentioned should already be implemented:

Issue metric	seastats key	Python implementation
RMSD	rmse	seastats/seastats/stats.py Lines 19 to 20 in 7ea1e05 def get_rmse(sim: pd.Series[float], obs: pd.Series[float]) -> float: return float(np.sqrt(get_mse(sim, obs))) where get_mse(sim, obs): seastats/seastats/stats.py Lines 15 to 16 in 7ea1e05 def get_mse(sim: pd.Series[float], obs: pd.Series[float]) -> float: return float(np.square(np.subtract(obs, sim)).mean())
URMSD (centered)	rms	seastats/seastats/stats.py Lines 41 to 43 in 7ea1e05 def get_rms(sim: pd.Series[float], obs: pd.Series[float]) -> float: crmsd = ((sim - sim.mean()) - (obs - obs.mean())) ** 2 return float(np.sqrt(crmsd.mean()))
Mean Bias (MB)	bias	seastats/seastats/stats.py Lines 11 to 12 in 7ea1e05 def get_bias(sim: pd.Series[float], obs: pd.Series[float]) -> float: return float(sim.mean() - obs.mean())
stdev(M)	sim_std	seastats/seastats/__init__.py Lines 182 to 183 in 7ea1e05 case "sim_std": stats["sim_std"] = sim.std()
stdev(O)	obs_std	seastats/seastats/__init__.py Lines 184 to 185 in 7ea1e05 case "obs_std": stats["obs_std"] = obs.std()
CC (Pearson)	cr	seastats/seastats/stats.py Lines 46 to 47 in 7ea1e05 def get_corr(sim: pd.Series[float], obs: pd.Series[float]) -> float: return float(sim.corr(obs))

So functionally, nothing new needs to be computed - all six already exist and be returned by calling get_stats(sim, obs, GENERAL_METRICS_ALL).

However I am conscious that these metrics are not following a conventional naming scheme.

Please let me know:

• if the metrics returned are correctly calculated
• if I should add those 6 metrics you mentioned to the GENERAL_METRICS:

  seastats/seastats/__init__.py

  Line 80 in 7ea1e05

  GENERAL_METRICS = ["bias", "rms", "rmse", "cr", "nse", "kge"]

• If I should rename some metrics for more clarity.

On my side:

• I will implement the VD metric you defined in consider adding "variance similarity" and "variance dissimilarity" to metrics #20

[DanCodigaNOAA]

Correct, URMSD is what Taylor diagrams show (distance of the symbol/marker from the 'perfect' star on the horizontal axis).
Good to know you did already have URMSD-- it was the naming that tripped me up.

My suggestion for naming the metrics and corresponding seasats keys would be RMSD, MB, URMSD, CC, VS, VD if possible. If you have a hard constraint against capitals in keys, then lowercase rmsd, mb, urmsd, cc, vs, vd should work. Using bias instead of mb would be fine, as would using cr instead of cc, if you prefer.

(I like D instead of E, in RMSD, because "difference" is more general than "error"-- in some cases you may be comparing observations to observations, collected by different methods, where "error" isn't really the right term, I don't think. )

Based on a quick look I think your calcs look fine but just be aware that there are two ways to do the calcs-- using population based statistics (with "delta deg of freedom" = 0, so you have n in the denom, where n is the sample size) or using sample based statistics (with ddof=1, so you have n-1 in the denominator).

Ideally the metrics calculation function can use population OR sample statistics, with one set to be the default (say, ddof=0 for population type stats) and the other implemented by passing in an input (say, ddof=1 for sample type stats). I think I would make population the default but the way to decide is based on the fact that the two methods converge as n values increase. So for low n it matters more, and sample could be the better choice.

Under GENERAL METRICS you have a mixture of error metrics and indices of agreement (nse and kge). You already have RMSD, URMSD, MB, and CC for error metrics so yes, I would suggest adding VS and VD there, as well as mean absolute error MAE. For indices of agreement you might want to add the Willmott Skill Score but only if you know people would use it (there are three versions of the Willmott Skill, original modified, and refined-- if you want to add them I can give equations and references to cite.)

Glad to know you use Taylor plots. I'm a big fan of them. I'm working on an improved diagram, actually.