Merge pull request #21 from kdorheim/mse_with_unc

Added ways to calculate MSE while accounting for uncertainty in temperature

results/default_stats.txt

@@ -4,11 +4,25 @@ Unnamed Hector core	NA	q10_rh	1.76	(unitless)
 Unnamed Hector core	NA	diff	2.38	cm2/s
 
 Mean of T, CO2 MSEs: 2.33
+Mean of T, CO2 NMSEs: 0.00204
+Mean of T, CO2 NMSEs (with unc): 0.00091
 
-CO2 MSE: 4.58
-T MSE:   0.0771
+Mean of smooth T, CO2 MSEs (k = 3 ): 2.32
+Mean of smooth T, CO2 NMSEs (k = 3 ): 0.00201
+Mean of smooth T, CO2 MSEs (k = 5 ): 2.32
+Mean of smooth T, CO2 NMSEs (k = 5 ): 0.00201
+Mean of smooth T, CO2 MSEs (k = 10 ): 2.32
+Mean of smooth T, CO2 NMSEs (k = 10 ): 0.00198
+
+CO2 MSE:        4.58
+T MSE:          0.0771
+T MSE with unc: 0.0343
 RMSE:    2.14
 
+CO2 NMSE: 2.74e-07
+T NMSE:   0.00409
+T NMSE with unc: 0.00182
+
 ***Key Metrics***
 TCRE: 1.51
 TCR:  1.77

results/unc_nmse.txt

@@ -0,0 +1,42 @@
+parameters	values
+beta	0.732
+q10_rh	2.64
+diff	2.4
+
+Objective Function Value: 0.000714
+Counts: 9
+Counts: 9
+
+CO2 MSE: 58.5
+T MSE:   0.0651
+RMSE:    7.65
+T MSE accounting for unc: 0.0268
+
+***Key Metrics***
+TCRE: 1.53
+TCR:  1.76
+
+***Historical Warming and ERF***
+GSAT Warming:              0.633
+Ocean Heat Content Change: 423
+Total Aerosol ERF:         -1.24
+WMGHG ERF:                 3.63
+Methane ERF:               0.538
+
+***Future Warming***
+scenario	start	end	GSAT
+ssp119	2021	2040	0.694
+ssp119	2041	2060	0.842
+ssp119	2081	2100	0.693
+ssp126	2021	2040	0.709
+ssp126	2041	2060	1.02
+ssp126	2081	2100	1.05
+ssp245	2021	2040	0.715
+ssp245	2041	2060	1.23
+ssp245	2081	2100	1.9
+ssp370	2021	2040	0.725
+ssp370	2041	2060	1.37
+ssp370	2081	2100	2.86
+ssp585	2021	2040	0.842
+ssp585	2041	2060	1.68
+ssp585	2081	2100	3.71

results/unc_nmse_big_box.txt

@@ -0,0 +1,93 @@
+parameters	values
+beta	1.08394143340344
+q10_rh	3.52
+diff	2
+
+Objective Function Value: 0.000642
+Counts: 17
+Counts: 17
+
+CO2 MSE: 233
+T MSE:   0.0592
+RMSE:    15.3
+T MSE accounting for unc: 0.0239
+
+***Key Metrics***
+TCRE: 1.51
+TCR:  1.82
+
+***Historical Warming and ERF***
+GSAT Warming:              0.538
+Ocean Heat Content Change: 347
+Total Aerosol ERF:         -1.24
+WMGHG ERF:                 3.31
+Methane ERF:               0.536
+
+***Future Warming***
+scenario	start	end	GSAT
+ssp119	2021	2040	0.651
+ssp119	2041	2060	0.754
+ssp119	2081	2100	0.601
+ssp126	2021	2040	0.666
+ssp126	2041	2060	0.933
+ssp126	2081	2100	0.916
+ssp245	2021	2040	0.672
+ssp245	2041	2060	1.14
+ssp245	2081	2100	1.73
+ssp370	2021	2040	0.682
+ssp370	2041	2060	1.28
+ssp370	2081	2100	2.7
+ssp585	2021	2040	0.805
+ssp585	2041	2060	1.61
+ssp585	2081	2100	3.58
+parameters	values
+beta	1.08393761240082
+q10_rh	3.52
+diff	2
+
+Objective Function Value: 0.000642
+Counts: 18
+Counts: 18
+
+CO2 MSE: 233
+T MSE:   0.0592
+RMSE:    15.3
+T MSE accounting for unc: 0.0239
+
+***Key Metrics***
+TCRE: 1.51
+TCR:  1.82
+
+***Historical Warming and ERF***
+GSAT Warming:              0.538
+Ocean Heat Content Change: 347
+Total Aerosol ERF:         -1.24
+WMGHG ERF:                 3.31
+Methane ERF:               0.536
+
+***Future Warming***
+scenario	start	end	GSAT
+ssp119	2021	2040	0.651
+ssp119	2041	2060	0.754
+ssp119	2081	2100	0.601
+ssp126	2021	2040	0.666
+ssp126	2041	2060	0.933
+ssp126	2081	2100	0.916
+ssp245	2021	2040	0.672
+ssp245	2041	2060	1.14
+ssp245	2081	2100	1.73
+ssp370	2021	2040	0.682
+ssp370	2041	2060	1.28
+ssp370	2081	2100	2.7
+ssp585	2021	2040	0.805
+ssp585	2041	2060	1.61
+ssp585	2081	2100	3.58
+parameters	values
+beta	1.0838541250268
+q10_rh	3.52
+diff	2
+
+Objective Function Value: 0.000642
+Counts: 20
+Counts: 20
+Convergence:  0

scripts/default_stats.R

@@ -6,6 +6,7 @@
 
 # Importing libraries
 library(hector)
+library(zoo)
 
 # Setting up file paths
 COMP_DATA_DIR <- file.path(here::here(), "comparison_data")
@@ -28,8 +29,8 @@ source(file.path(SCRIPTS_DIR, "major_functions.R"))
 
 
 ### Getting observational data ###
-co2_data <- get_co2_data(CO2_PATH)
-temp_data <- get_temp_data(TEMP_PATH)
+co2_data <- get_co2_data(CO2_PATH, include_unc = T)
+temp_data <- get_temp_data(TEMP_PATH, include_unc = T)
 obs_data <- rbind(co2_data, temp_data)
 
 
@@ -52,28 +53,85 @@ hector_data <- run_hector(ini_file = INI_FILE,
                           params = NULL, 
                           vals = best_pars, 
                           yrs = 1750:2014, 
-                          vars = c(GMST(), CONCENTRATIONS_CO2()))
+                          vars = c(GMST(), CONCENTRATIONS_CO2()),
+                          include_unc = T)
 
 
+# Outputting overall MSEs
 write_metric("Mean of T, CO2 MSEs:", 
              mean_T_CO2_mse(obs_data, hector_data), 
              OUTPUT)
+write_metric("Mean of T, CO2 NMSEs:",
+             mean_T_CO2_nmse(obs_data, hector_data),
+             OUTPUT)
+write_metric("Mean of T, CO2 NMSEs (with unc):",
+             mean_T_CO2_nmse_unc(obs_data, hector_data),
+             OUTPUT)
+write("", OUTPUT, append = TRUE)
+
+# Outputting smoothed MSEs
+for (k in c(3, 5, 10)) {
+  smoothed_temps <- rollmean(temp_data$value, 
+                             k = k, 
+                             align = "center",
+                             fill = NA)
+  smoothed_temp_data <- temp_data
+  smoothed_temp_data$variable <- "Smooth T"
+  smoothed_temp_data$value <- smoothed_temps
+
+  new_data <- rbind(co2_data, temp_data, smoothed_temp_data)
+
+  write_metric(paste("Mean of smooth T, CO2 MSEs (k =", k, "):"),
+               smooth_T_CO2_mse(new_data, hector_data),
+               OUTPUT)
+
+  write_metric(paste("Mean of smooth T, CO2 NMSEs (k =", k, "):"),
+               smooth_T_CO2_nmse(new_data, hector_data),
+               OUTPUT)
+}
 write("", OUTPUT, append = TRUE)
 
 # Outputting individual MSEs
 T_mse <- get_var_mse(obs_data = obs_data, 
                      hector_data = hector_data, 
                      var = GMST(), 
                      yrs = 1850:2014)
+T_mse_unc <- get_var_mse_unc(obs_data = obs_data, 
+                             hector_data = hector_data, 
+                             var = GMST(), 
+                             yrs = c(1850:2014),
+                             mse_fn = mse_unc)
 CO2_mse <- get_var_mse(obs_data = obs_data, 
                        hector_data = hector_data, 
                        var = CONCENTRATIONS_CO2(), 
                        yrs = c(1750, 1850:2014))
 
-write_metric("CO2 MSE:", CO2_mse, OUTPUT)
-write_metric("T MSE:  ", T_mse, OUTPUT)
+# Getting NMSEs
+T_nmse <- get_var_mse(obs_data = obs_data, 
+                     hector_data = hector_data, 
+                     var = GMST(), 
+                     yrs = 1850:2014,
+                     mse_fn = nmse)
+T_nmse_unc <- get_var_mse_unc(obs_data = obs_data, 
+                              hector_data = hector_data, 
+                              var = GMST(), 
+                              yrs = c(1850:2014),
+                              mse_fn = nmse_unc)
+CO2_nmse <- get_var_mse(obs_data = obs_data, 
+                       hector_data = hector_data, 
+                       var = CONCENTRATIONS_CO2(), 
+                       yrs = c(1750, 1850:2014),
+                       mse_fn = nmse)
+
+write_metric("CO2 MSE:       ", CO2_mse, OUTPUT)
+write_metric("T MSE:         ", T_mse, OUTPUT)
+write_metric("T MSE with unc:", T_mse_unc, OUTPUT)
 write_metric("RMSE:   ", sqrt(mean(CO2_mse, T_mse)), OUTPUT) # not 100% sure this is how we want to calculate this
 write("", OUTPUT, append = TRUE)
+write_metric("CO2 NMSE:", CO2_nmse, OUTPUT)
+write_metric("T NMSE:  ", T_nmse, OUTPUT)
+write_metric("T NMSE with unc:", T_nmse_unc, OUTPUT)
+write("", OUTPUT, append = TRUE)
 
 ### Outputting table metrics ###
 calc_table_metrics(NULL, NULL, OUTPUT)