Model ensemble mean, bias, and RMSE

This example shows how to calculate: 1) the model ensemble mean 2) the model simple bias 3) the model root-mean-squared error

Let’s start by creating 3 example datasets. All 3 represent annual mean temperatures from 1979 to 2014

• dat1 = observation
• dat2 = model A
• dat3 = model B

dat1 = c(23.27, 23.09, 22.93, 23.50, 23.65, 23.98, 24.09,
         23.51, 24.43, 24.14, 23.40, 23.71, 23.60, 23.67,
         24.76, 23.66, 22.99, 23.42, 22.92, 24.55, 24.32,
         23.49, 24.33, 24.10, 24.00, 23.72, 23.65, 23.76,
         23.46, 23.17, 24.02, 23.35, 22.96, 24.80, 24.41,
         24.91)

dat2 = c(24.11, 22.43, 22.82, 23.65, 24.63, 23.22, 23.08,
         23.37, 23.22, 23.02, 23.81, 22.80, 23.20, 22.77,
         23.98, 23.15, 22.75, 23.09, 23.88, 24.01, 23.15,
         24.00, 23.71, 23.50, 23.53, 23.51, 23.41, 24.10,
         23.35, 24.44, 24.53, 23.44, 23.91, 24.16, 24.28,
         23.66)

dat3 = c(23.77, 22.69, 22.86, 23.59, 24.24, 23.52, 23.49,
         23.43, 23.70, 23.46, 23.65, 23.17, 23.36, 23.13,
         24.29, 23.36, 22.85, 23.22, 23.49, 24.22, 23.62,
         23.80, 23.96, 23.74, 23.72, 23.59, 23.51, 23.97,
         23.39, 23.93, 24.32, 23.40, 23.53, 24.41, 24.33,
         24.16)

Let’s convert the 3 datasets into time series

#convert monthly data into time series 
library(xts)
nt = length(dat1)
dates = seq(as.Date("1979-01-01"), length=nt, by="years")
obs = xts(dat1,order.by=dates,frequency=12)
mdl.A = xts(dat2,order.by=dates,frequency=12)
mdl.B = xts(dat3,order.by=dates,frequency=12)

The model ensemble mean is simply the mean of all model results. In this example, we only have 2 models: A and B

#ensemble mean
ens_mean = (mdl.A+mdl.B)/2

The model simple bias is the mean of the differences between model and observation

bias.A = mean(mdl.A-obs)
bias.B = mean(mdl.B-obs)

print("Bias:")

## [1] "Bias:"

print(bias.A)

## [1] -0.2236111

print(bias.B)

## [1] -0.1347222

And the model RSME is the square-root of the mean of the squared differences between model and observation

rmse.A = sqrt(mean((mdl.A-obs)^2))
rmse.B = sqrt(mean((mdl.B-obs)^2))

print("RMSE:")

## [1] "RMSE:"

print(rmse.A)

## [1] 0.7126379

print(rmse.B)

## [1] 0.4276713

You can also calculate the bias and RMSE for the ensemble mean.