Reliable estimation of prediction uncertainty for physico-chemical property models

TL;DR

This paper introduces a bootstrapping-based calibration method for physico-chemical property models, providing reliable uncertainty estimates and systematic error diagnostics, demonstrated on Moessbauer isomer shift predictions with density functional theory.

Contribution

It presents the first statistically rigorous calibration analysis for Moessbauer spectroscopy, enhancing uncertainty estimation and model validation for physico-chemical property predictions.

Findings

Prediction uncertainties vary with the choice of density functional.

Systematic error diagnostics improve model reliability.

Reference data composition significantly affects model parameters.

Abstract

The predictions of parameteric property models and their uncertainties are sensitive to systematic errors such as inconsistent reference data, parametric model assumptions, or inadequate computational methods. Here, we discuss the calibration of property models in the light of bootstrapping, a sampling method akin to Bayesian inference that can be employed for identifying systematic errors and for reliable estimation of the prediction uncertainty. We apply bootstrapping to assess a linear property model linking the 57Fe Moessbauer isomer shift to the contact electron density at the iron nucleus for a diverse set of 44 molecular iron compounds. The contact electron density is calculated with twelve density functionals across Jacob's ladder (PWLDA, BP86, BLYP, PW91, PBE, M06-L, TPSS, B3LYP, B3PW91, PBE0, M06, TPSSh). We provide systematic-error diagnostics and reliable, locally resolved uncertainties for isomer-shift predictions. Pure and hybrid density functionals yield average prediction uncertainties of 0.06-0.08 mm/s and 0.04-0.05 mm/s, respectively, the latter being close to the average experimental uncertainty of 0.02 mm/s. Furthermore, we show that both model parameters and prediction uncertainty depend significantly on the composition and number of reference data points. Accordingly, we suggest that rankings of density functionals based on performance measures (e.g., the coefficient of correlation, r2, or the root-mean-square error, RMSE) should not be inferred from a single data set. This study presents the first statistically rigorous calibration analysis for theoretical Moessbauer spectroscopy, which is of general applicability for physico-chemical property models and not restricted to isomer-shift predictions. We provide the statistically meaningful reference data set MIS39 and a new calibration of the isomer shift based on the PBE0 functional.

Figures

Tables

Table 1. Table 1: Reference set of molecular iron compounds compiled for this study (to be continued).

no.	compound111Abbreviations in alphabetical order: ac, CH3CO-; ArNC, 2,6-Me2-C6H3NC; ArS, 2,6-Me2-C6H3S; $\beta$-DKI${}^{\text{Me}}$, bulky $\beta$-diketiminate; cyclam, 1,4,8,11-tetraazacyclotetradecane; dae, N–CH2–CH2–N; Dipp, 2,6-diisopropylphenyl; DTSQ, bis(1,2-dithiosquarato-S,S’); HBpz3, hydrotris-1-(pyrazolyl)borate; MAC∗, tetraanion of 1,4,8,11-tetraaza-13,13-diethyl-2,2,5,5,7,7,10,10-octamethyl-3,6,9,12,14-pentaoxocyclotetradecane; MBTHx, bis(N-methylbenzothiohydroxamato) anion; Me3TACN, 1,4,7-trimethyl-1,4,7-triazacyclononane; OEC, dianion of trans-7,8-dihydro-octaethylporphyrin; OEP, dianion of octaethylporphyrin; opda, o-phenylenediamine; PhTt${}^{\text{t}\text{Bu}}$, phenyltris((tert-butylthio)-methyl)borate; por, porphyrin; PyO, 6-methyl-2-pyridonate; PyS, 6-methyl-2-thiolate; pyS4, dianion of 2,6-bis(2-mercaptophenylthiomethyl)pyridine; “S2”, 1,2-benzenedithiolato-S,S’ dianion; ${}^{\text{t}\text{Bu}}$py, 4-tert-butylpyridine; tmc, 1,4,8,11-tetramethyl-1,4,8,11-tetraazacyclotetradecane; trop, 5H-dibenzo[a,d]cyclo-hepten-5-yl.	$2S+1$	ox. state	CN	${\delta }_{\text{exp}}$ (mm s-1)222All experimental isomer shifts reported in this study refer to bulk $\alpha$-iron at room temperature.	$T_{\text{exp}}$ (K)	ref. ${\delta }_{\text{exp}}$	ref. coords.	code333The six-letter codes refer to the Cambridge Structural Database identifiers. The N09-X code refers to the structure #X in ascending order as found in the Supporting Information of the study by Römelt, Ye, and Neese [49]. The R13-Y code refers to the structure #Y as found in Table 1 of our previous isomer-shift calibration study [54]. For two compounds, no code is available (–).
1444This compound refers to an inconsistent input–target pair according to the jackknife-after-bootstrapping method.	Fe(NCS)${}_{}{}^{2-}{}_4{}^{}$	5	+2	4	0.97($-$)	4.2	[77]	[78]	KEFFEG
2\footreffn:critical	Fe(opda)2Cl2	5	+2	6	0.88($-$)	80	[79]	[80]	FUJQOQ
3555The molecular structure was truncated prior to structure optimization.	Fe(${}^{\text{t}\text{Bu}}$py)${}_2{}^{}\beta$-DKI${}^{\text{Me}}$	4	+1	4	0.79($-$)	80	[68]	[68]	YEZSAZ
4\footreffn:trunc	Fe(PhTt${}^{\text{t}\text{Bu}}$)(PEt3)	4	+1	4	0.76(3)	4.5	[65]	[65]	VIYHER
5\footreffn:trunc	Fe(${\eta }^6$-benzene)$\beta$-DKI${}^{\text{Me}}$	2	+1	8	0.68($-$)	80	[69]	[69]	$-$
6	Fe(DTSQ)${}_{}{}^{2-}{}_2{}^{}$	5	+2	4	0.67($-$)	4.2	[81]	[81]	PTSQFE10
7\footreffn:critical	Fe(por)(O2)-	6	+3	5	0.67($-$)	4.2	[82]	[49]	N09-14
8	Fe(SPh)${}_{}{}^{2-}{}_4{}^{}$	5	+2	4	0.66($-$)	4.2	[81]	[81]	PTHPFE10
9	Fe(OEC)	3	+2	4	0.63(1)	4.2	[83]	[83]	BUYKUB10
10\footreffn:trunc	Fe(PhTt${}^{\text{t}\text{Bu}}$)(PhCCPh)	4	+1	5	0.62(3)	4.5	[65]	[65]	VIYHOB
11	Fe(OEP)	3	+2	4	0.62(1)	4.2	[83]	[83]	DEDWUE
12\footreffn:trunc	Fe(PMe)(PhBP${}_3{}^{}{}_{}{}^{\text{i}\text{Pr}}$)	4	+1	4	0.57(2)	80	[64]	[84]	QAJNUL
13\footreffn:critical	Fe2O(HBpz3)2(OAc)2	1	+3	6	0.52(3)	4.2	[85]	[85]	CACZIP10
14\footreffn:trunc	Fe(HCCPh)(HC(C${}^{\text{t}\text{Bu}}$N–Dipp)2)-	4	+1	4	0.50(1)	4.2	[63]	[63]	MATVOT
15	Fe2O(OAc)2(Me3TACN)${}_{}{}^{2+}{}_2{}^{}$	1	+3	6	0.47(3)	4.2	[86]	[86]	DIBXAN10
16	Fe(SMe3)(pyS4)	1	+2	6	0.44(2)	80	[87]	[49]	N09-24
17	FeCl(MBTHx)2	6	+3	5	0.43(1)	4.2	[88]	[88]	CELVEU
18\footreffn:trunc	Fe(SiP${}_3{}^{}{}_{}{}^{\text{i}\text{Pr}}$)N2	2	+1	5	0.38($-$)	77	[66]	[89]	NIHQIF
19	Fe(PH3)(pyS4)	1	+2	6	0.34(2)	80	[87]	[49]	N09-23
20	Fe(NO)(pyS4)	2	+2	6	0.33(2)	80	[87]	[49]	N09-22
21	trans-Fe(cyclam)(N3)${}_{}{}^{+}{}_2{}^{}$	2	+3	6	0.28(6)	80	[90]	[49]	N09-6
22	LiFe(trop2dae)	2	+1	6	0.28(1)	77	[70]	[70]	RUDHAB

Table 2. Continuation of Table 1 .

no.	compound666Abbreviations in alphabetical order: ac, CH3CO-; ArNC, 2,6-Me2-C6H3NC; ArS, 2,6-Me2-C6H3S; $\beta$-DKI${}^{\text{Me}}$, bulky $\beta$-diketiminate; cyclam, 1,4,8,11-tetraazacyclotetradecane; dae, N–CH2–CH2–N; Dipp, 2,6-diisopropylphenyl; DTSQ, bis(1,2-dithiosquarato-S,S’); HBpz3, hydrotris-1-(pyrazolyl)borate; MAC∗, tetraanion of 1,4,8,11-tetraaza-13,13-diethyl-2,2,5,5,7,7,10,10-octamethyl-3,6,9,12,14-pentaoxocyclotetradecane; MBTHx, bis(N-methylbenzothiohydroxamato) anion; Me3TACN, 1,4,7-trimethyl-1,4,7-triazacyclononane; OEC, dianion of trans-7,8-dihydro-octaethylporphyrin; OEP, dianion of octaethylporphyrin; opda, o-phenylenediamine; PhTt${}^{\text{t}\text{Bu}}$, phenyltris((tert-butylthio)-methyl)borate; por, porphyrin; PyO, 6-methyl-2-pyridonate; PyS, 6-methyl-2-thiolate; pyS4, dianion of 2,6-bis(2-mercaptophenylthiomethyl)pyridine; “S2”, 1,2-benzenedithiolato-S,S’ dianion; ${}^{\text{t}\text{Bu}}$py, 4-tert-butylpyridine; tmc, 1,4,8,11-tetramethyl-1,4,8,11-tetraazacyclotetradecane; trop, 5H-dibenzo[a,d]cyclo-hepten-5-yl.	$2S+1$	ox. state	CN	${\delta }_{\text{exp}}$ (mm s-1)777All experimental isomer shifts reported in this study refer to bulk $\alpha$-iron at room temperature.	$T_{\text{exp}}$ (K)	ref. ${\delta }_{\text{exp}}$	ref. coords.	code888The six-letter codes refer to the Cambridge Structural Database identifiers. The N09-X code refers to the structure #X in ascending order as found in the Supporting Information of the study by Römelt, Ye, and Neese [49]. The R13-Y code refers to the structure #Y as found in Table 1 of our previous isomer-shift calibration study [54]. For two compounds, no code is available (–).
23\footreffn:trunc	Fe(OEP)CO	1	+2	5	0.27($-$)	4.2	[91]	[91]	YEQPOA
24	Fe(SEt)${}_{}{}^{-}{}_4{}^{}$	6	+3	4	0.25($-$)	4.2	[92]	[92]	CANDAW10
25\footreffn:trunc	Fe(SiP${}_3{}^{}{}_{}{}^{\text{i}\text{Pr}}$)CO	2	+1	5	0.21($-$)	77	[67]	[67]	UTOXOR
26	NaFe(trop2dae)	2	+1	6	0.20(1)	77	[70]	[70]	RUDGUU
27\footreffn:trunc	Fe(OEC)Cl	3	+4	5	0.19($-$)	77	[93]	[93]	SUMWUS
28\footreffn:critical	Fe(NO)2(S(p-Me)Ph)${}_{}{}^{-}{}_2{}^{}$	6	+3	4	0.18(2)	4.2	[94]	[94]	SONMUE
29	Fe(PPh3)2(“S2”)2	3	+4	6	0.17(1)	4.2	[95]	[95]	SOCVUB
30	Fe$=$O(tmc)(NCCH3)2+	3	+4	6	0.17(1)	4.2	[96]	[49]	N09-18
31	Fe(PPh3)(“S2”)2	3	+4	5	0.12(1)	4.2	[95]	[95]	SOCWAI
32	Fe(${\eta }^4$-butadiene)(CO)3	1	0	7	0.12($-$)	4.2	[97]	[98]	$-$
33	Fe(PyO)I(CO)2PPh3	1	+2	6	0.10(2)	55	[99]	[99]	PUCFID
34	Fe(PyS)I(CO)2PPh3	1	+2	6	0.10(2)	77	[100]	[54]	R13-L
35	Fe(PyO)(ArS)(CO)2PPh3	1	+2	6	0.06(2)	100	[100]	[54]	R13-K
36	Fe2(PyS)2(ac)2(CO)4	1	+2	6	0.06(1)	77	[100]	[101]	FUMZAP
37	Fe(NO)(pyS4)+	1	+2	6	0.04(2)	80	[87]	[49]	N09-21
38	Fe(SC5H4N-CO)I(CO)2(ArNC)	1	+2	6	0.01(2)	77	[100]	[100]	YAKZAN
39	Fe(PyS)(ac)(CO)2(ArNC)	1	+2	6	0.00(2)	100	[100]	[54]	R13-O
40	Fe(PyS)(ac)(CO)2PPh3	1	+2	6	0.00(2)	77	[100]	[101]	FUMZET
41	Fe(PyS)(ac)(CN)(CO)2	1	+2	6	$-$0.02(2)	77	[100]	[54]	R13-Q
42\footreffn:trunc	FeCl(${\eta }^4$-MAC∗)-	5	+4	5	$-$0.04(2)	4.2	[102]	[102]	JESGUJ
43	Fe(CO)5	1	0	5	$-$0.09(1)	78	[103]	[49]	N09-20
44	Fe(OEC)C6H5	3	+4	5	$-$0.11($-$)	77	[93]	[93]	SUMXED

Table 3. Table 2: Performance measures (RMSE, R632, RMPV, RLOO) of density functionals regarding isomer-shift prediction. The reference data set is the original one ( N = 44 𝑁 44 N=44 ). RMSE, R632, and RLOO were calculated on the basis of ordinary (unweighted) least-squares regression. The RMPV was obtained from Bayesian linear regression based on the evidence approximation. For the calculation of the .632 estimator, R632, B = 10 5 𝐵 superscript 10 5 B=10^{5} bootstrap samples were generated. All performance measures are reported in mm s -1 . The corresponding ACED values were obtained with the def2-TZVP basis set for ligand atoms.

functional	RMSE	R632	RMPV	RLOO
PWLDA	0.09	0.10	0.09	0.10
BP86	0.09	0.09	0.09	0.09
BLYP	0.10	0.10	0.10	0.10
PW91	0.09	0.09	0.09	0.09
PBE	0.09	0.09	0.09	0.09
M06-L	0.08	0.08	0.08	0.08
TPSS	0.08	0.09	0.09	0.09
B3LYP	0.07	0.07	0.07	0.07
B3PW91	0.07	0.07	0.07	0.07
PBE0	0.07	0.07	0.07	0.07
M06	0.07	0.07	0.07	0.07
TPSSh	0.07	0.08	0.08	0.08

Table 4. Table 3: Parameters of the linear isomer-shift model, δ ​ ( ρ absorber , 𝐰 ) 𝛿 subscript 𝜌 absorber 𝐰 \delta(\rho_{\text{absorber}},\mathbf{w}) , optimized with different calibration procedures based on linear least-squares regression. Only those reference isomer shifts with reported experimental uncertainty ( N = 30 𝑁 30 N=30 ) were employed. The corresponding ACED values were obtained from B3LYP/def2-TZVP calculations. Standard deviations in parentheses were obtained from bootstrapping ( B = 10 4 𝐵 superscript 10 4 B=10^{4} ).

calibration procedure	$w_0$ (mm s-1)	$w_1$ (mm s-1 bohr3)
ordinary (unweighted)	4454.80(15303)	$-$0.38(1)
weighted	4397.39(14260)	$-$0.37(1)
iteratively reweighted	4451.15(14963)	$-$0.38(1)

Table 5. Table 4: Estimated MPU (R632) for four property models of increasing polynomial degree, M 𝑀 M , on the basis of B = 10 4 𝐵 superscript 10 4 B=10^{4} bootstrap samples ( N = 44 𝑁 44 N=44 ), respectively. The corresponding ACED results were obtained from B3LYP/def2-TZVP calculations. The linear and quadratic models reveal the lowest MPU.

property model	$M$	R632 (mm s-1)
linear	1	0.07
quadratic	2	0.07
cubic	3	0.08
quartic	4	0.10

Table 6. Table 5: Statistics of model parameters corresponding to the distributions shown in Fig.  2 . w ¯ m subscript ¯ 𝑤 𝑚 \bar{w}_{m} and σ w m subscript 𝜎 subscript 𝑤 𝑚 \sigma_{w_{m}} ( m = { 0 , 1 } 𝑚 0 1 m=\{0,1\} ) correspond to mean and standard deviation of the inferred parameter distributions, respectively. corr w 0 , w 1 subscript corr subscript 𝑤 0 subscript 𝑤 1 \text{corr}_{w_{0},w_{1}} refers to the correlation between w 0 subscript 𝑤 0 w_{0} and w 1 subscript 𝑤 1 w_{1} as measured by the square root of r 2 superscript 𝑟 2 r^{2} . For better comparability, the results have not been rounded to the experimental resolution (0.01 mm s -1 ).

inference method	$N$	${\overline{w}}_0$	${\overline{w}}_1$	${\sigma }_{w_0}$	${\sigma }_{w_1}$	${\text{corr}}_{w_0,w_1}$
bootstrapping	44	0.3295	$-$0.3674	0.0101	0.0157	$-$0.4720
	39	0.2888	$-$0.3805	0.0056	0.0082	$-$0.0776
Bayesian	44	0.3288	$-$0.3658	0.0102	0.0137	0.0000
	39	0.2886	$-$0.3801	0.0057	0.0085	0.0000

Table 7. Table 6: Performance measures (RMSE, R632, RMPV, RLOO) of density functionals regarding isomer-shift prediction. The reference data set is the original one with the inconsistent data points (#1, #2, #7, #13, #28, cf. Table 1 and Fig.  1 , left) removed ( N = 39 𝑁 39 N=39 ). RMSE, R632, and RLOO were calculated on the basis of ordinary (unweighted) least-squares regression. The RMPV was obtained from Bayesian linear regression based on the evidence approximation. For the calculation of the .632 estimator, R632, B = 10 5 𝐵 superscript 10 5 B=10^{5} bootstrap samples were generated. All performance measures are reported in mm s -1 . The corresponding ACED values were obtained with the def2-TZVP basis set for ligand atoms.

functional	RMSE	R632	RMPV	RLOO
PWLDA	0.07	0.07	0.07	0.07
BP86	0.07	0.07	0.07	0.07
BLYP	0.08	0.08	0.08	0.08
PW91	0.07	0.07	0.07	0.07
PBE	0.07	0.07	0.07	0.07
M06-L	0.06	0.06	0.06	0.06
TPSS	0.06	0.06	0.06	0.06
B3LYP	0.03	0.04	0.04	0.04
B3PW91	0.03	0.04	0.04	0.04
PBE0	0.03	0.04	0.04	0.04
M06	0.05	0.05	0.05	0.05
TPSSh	0.04	0.04	0.04	0.04

Table 8. Table 7: Parameters of the linear isomer-shift model, δ ​ ( ρ absorber , 𝐰 ) 𝛿 subscript 𝜌 absorber 𝐰 \delta(\rho_{\text{absorber}},\mathbf{w}) , optimized with different calibration procedures based on linear least-squares regression. Only those reference isomer shifts with reported experimental uncertainty and consistent status ( N = 28 𝑁 28 N=28 ) were employed. The corresponding ACED values were obtained from B3LYP/def2-TZVP calculations. Standard deviations in parentheses were obtained from bootstrapping ( B = 10 4 𝐵 superscript 10 4 B=10^{4} ).

calibration procedure	$w_0$ (mm s-1)	$w_1$ (mm s-1 bohr3)
ordinary (unweighted)	4510.05(11754)	$-$0.38(1)
weighted	4434.25(12639)	$-$0.38(1)
iteratively reweighted	4495.55(10932)	$-$0.38(1)

Equations

$$\text{MSE}\equiv\text{MSE}_{\mathcal{D},\mathbf{w}}=\frac{1}{N}\sum_{n=1}^{N}\big{(}y_{n}-f(\mathbf{x}_{n},\mathbf{w})\big{)}^{2}\,$$

$$R^2\equiv R_{\mathcal{D}}^2=1-N\frac{{\min }_{\mathbf{w}}({\text{MSE}}_{\mathcal{D},\mathbf{w}})}{{\sum }_{n=1}^N(y_n-\bar{y}{)}^2},$$

$$\text{MSE}_{\mathcal{D},\mathbf{w}_{b}^{\ast}}=\frac{1}{N}\sum_{n=1}^{N}\big{(}y_{n}-f(\mathbf{x}_{n},\mathbf{w}_{b}^{\ast})\big{)}^{2}\geq\text{MSE}_{\mathcal{D}}\,.$$

$$\mathbb{E}[\text{MSE}]=\frac{1}{B}\sum _{b=1}^B{\text{MSE}}_{\mathcal{D},{\mathbf{w}}_b^{\ast }}\ge {\text{MSE}}_{\mathcal{D}},$$

$${\bar{\mathbf{w}}}^{\ast }=\frac{1}{B}\sum _{b=1}^B{\mathbf{w}}_b^{\ast }$$

$${\sigma }_{{\mathbf{w}}^{\ast }}^2=\frac{1}{B-1}\sum _{b=1}^B({\mathbf{w}}_b^{\ast }-{\bar{\mathbf{w}}}^{\ast })({\mathbf{w}}_b^{\ast }-{\bar{\mathbf{w}}}^{\ast }{)}^{\top },$$

$$\mathbb{E}[\text{MSE}{]}_{.632}=0.368{\text{MSE}}_{\mathcal{D}}+0.632\frac{1}{N}\sum _{n=1}^N\mathbb{E}[\text{MSE}{]}_{-n},$$

$$\mathbb{E}[\text{MSE}]_{-n}=\frac{1}{|\mathcal{B}_{-n}|}\sum_{b\in\mathcal{B}_{-n}}\big{(}y_{n}-f(\mathbf{x}_{n},\mathbf{w}_{b}^{\ast})\big{)}^{2}\,.$$

$$\mathbb{E}[\text{MSE}]_{\text{LOO}}=\frac{1}{N}\sum_{n=1}^{N}\big{(}y_{n}-f(\mathbf{x}_{n},\mathbf{w}_{\mathcal{D}_{-n}})\big{)}^{2}\,,$$

$$f({\mathbf{x}}_n,\mathbf{w})=\sum _{m=0}^M{w}_m(x_m{)}_n={\mathbf{x}}_n^{\top }\mathbf{w},$$

$$\mathbf{w}_{\mathcal{D}}=\Big{(}\mathbf{X}^{\top}\mathbf{X}\Big{)}^{-1}\mathbf{X}^{\top}\mathbf{y}\,,$$

$$\text{WMSE}\equiv\text{WMSE}_{\mathcal{D},\mathbf{w},\mathbf{U}}=\frac{1}{N}\sum_{n=1}^{N}u_{n}^{-2}\big{(}y_{n}-f(\mathbf{x}_{n},\mathbf{w_{\mathbf{U}}})\big{)}^{2}\,,$$

$$\mathbf{U}=\left(\begin{array}{ccc}{u}_1^2& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & u_N^2\end{array}\right).$$

$$\mathbf{w}_{\mathcal{D},\mathbf{U}}=\Big{(}\mathbf{X}^{\top}\mathbf{U}^{-1}\mathbf{X}\Big{)}^{-1}\mathbf{X}^{\top}\mathbf{U}^{-1}\mathbf{y}\,.$$

$$⟨u⟩\equiv \sqrt{\frac{1}{N}\sum _{n=1}^N{u}_n^2}.$$

$$\mathbf{w}_{b}^{\ast}=\Big{(}(\mathbf{X}_{b}^{\ast})^{\top}(\mathbf{U}_{b}^{\ast})^{-1}\mathbf{X}_{b}^{\ast}\Big{)}^{-1}(\mathbf{X}_{b}^{\ast})^{\top}(\mathbf{U}_{b}^{\ast})^{-1}\mathbf{y}_{b}^{\ast}\,,$$

$$u({\mathbf{x}}_0)=\sqrt{\frac{N}{N-M-1}{\text{MSE}}_{\mathcal{D},{\bar{\mathbf{w}}}^{\ast }}+{\mathbf{x}}_0^{\top }{\sigma }_{{\mathbf{w}}^{\ast }}^2{\mathbf{x}}_0}.$$

$$\mathbb{E}[\text{MSE}{]}_{.632}\approx \frac{1}{N}\sum _{n=1}^N{u}^2({\mathbf{x}}_n),$$

$$\text{MPV}=\frac{1}{N}\sum _{n=1}^N{s}^2({\mathbf{x}}_n),$$

$$\delta_{\text{exp}}=g(r)\Bigg{(}\frac{\Delta r}{r}\Bigg{)}[\rho_{\text{absorber}}-\rho_{\text{source}}]\,.$$

$$\delta ({\rho }_{\text{absorber}},\mathbf{w})=w_0+w_1{\rho }_{\text{absorber}}\approx {\delta }_{\text{exp}}.$$

$${\delta }_M({\rho }_{\text{absorber}},\mathbf{w})=\sum _{m=0}^M{w}_m{\rho }_{\text{absorber}}^m,$$

$$\Delta w_{0,n}\equiv \frac{|{\bar{w}}_{0,-n}^{\ast }-{\bar{w}}_0^{\ast }|}{{\sum }_{n=1}^N|{\bar{w}}_{0,-n}^{\ast }-{\bar{w}}_0^{\ast }|},$$

$$\Delta w_{1,n}\equiv \frac{|{\bar{w}}_{1,-n}^{\ast }-{\bar{w}}_1^{\ast }|}{{\sum }_{n=1}^N|{\bar{w}}_{1,-n}^{\ast }-{\bar{w}}_1^{\ast }|}.$$

$$\Delta\delta_{n}\equiv\frac{\sqrt{\sum_{i=1}^{N}\Big{(}\delta(\rho_{i,\text{absorber}},\bar{\mathbf{w}}^{\ast}_{-n})-\delta(\rho_{i,\text{absorber}},\bar{\mathbf{w}}^{\ast})\Big{)}^{2}}}{\sum_{n=1}^{N}\sqrt{\sum_{i=1}^{N}\Big{(}\delta(\rho_{i,\text{absorber}},\bar{\mathbf{w}}^{\ast}_{-n})-\delta(\rho_{i,\text{absorber}},\bar{\mathbf{w}}^{\ast})\Big{)}^{2}}}\,.$$

$$\delta ({\rho }_0,{\bar{\mathbf{w}}}_{\text{PBE0}}^{\ast })$$

$$\bar{\rho }$$

$${\sigma }_{{\mathbf{w}}_{\text{PBE0}}^{\ast }}^2=\left(\begin{array}{cc}2.959(4)\times 10^{-5}& 5.994(48)\times 10^{-6}{\text{ bohr}}^3\\ 5.994(48)\times 10^{-6}{\text{ bohr}}^3& 6.641(9)\times 10^{-5}{\text{ bohr}}^6\end{array}\right){\text{mm}}^2{\text{s}}^{-2},$$

$$\frac{N}{N-M-1}{\text{MSE}}_{\mathcal{D},{\bar{\mathbf{w}}}_{\text{PBE0}}^{\ast }}=1.185(1)\times 10^{-3}{\text{ mm}}^2{\text{ s}}^{-2}.$$

$$u({\mathbf{\rho }}_0)=\sqrt{\frac{N}{N-M-1}{\text{MSE}}_{\mathcal{D},{\bar{\mathbf{w}}}_{\text{PBE0}}^{\ast }}+{\mathbf{\rho }}_0^{\top }{\sigma }_{{\mathbf{w}}_{\text{PBE0}}^{\ast }}^2{\mathbf{\rho }}_0}.$$

$${\mathbf{w}}_{\text{MAP}}=\beta \mathbf{S}{\mathbf{X}}^{\top }\mathbf{y},$$

Full text

ETH Zürich, Laboratorium für Physikalische Chemie,

Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland.

1 Introduction

Predicting observables from a combination of scientific (knowledge-based) and statistical (data-based) information is at the heart of any parameteric property model applied in computational science.[1] In chemical physics, property models are applied whenever it is unfeasible or too demanding to calculate a target observable for a desired range of chemical systems with benchmark methods. To resolve this issue, the target observable is represented by a property model, which is a parametric representation of the former. The statistical variables of a property model (parameters) represent its unknown part, and their optimization (calibration of the model) requires both a reference data set and an objective, e.g., minimization of the mean squared error. The reference data set comprises pairs of values for (a) the target observable (obtained from measurements or benchmark results) and (b) the corresponding input variable. The input variable can represent the target observable itself or another physically motivated variable. It is also possible that the input variable is a vectorial quantity representing, e.g., nuclear coordinates.

The purpose of calibration is the estimation of parameter values that maximize the transferability of a property model to measurements or benchmark results of its target observable not included in the reference data set. In the field-specific literature, there are numerous applications of property models, e.g., for the calibration of force fields[2, 3], exchange–correlation density functionals[4, 5], dispersion-corrected potentials for density functional theory[6], semi-empirical electronic structure methods[7], vibrational frequencies[8, 9, 10], kinetic models[11, 12, 13], ionization potentials[14], thermochemical properties [15], properties of semi-conductors and insulators[16], linear free energy relationships[17], or melting point models[18], to name only a few.

Concomitant with the maximization of transferability is the assessment of model prediction uncertainty (MPU), i.e., the expected random deviation of a prediction from a measurement or benchmark result. MPU can be estimated analytically (or at least iteratively) if certain parametric assumptions are made on the population distribution underlying the reference data set. For a continuous variable, the most likely of all other parametric population distributions is the normal distribution, which is parameterized by mean and variance.[19] In that case, Bayesian inference is an efficient way to estimate MPU,[19] which we will also address in this paper. Generally, the more input–target pairs are included in the reference data set, the more reliable it does represent the underlying population distribution, and the less ambiguous is the selection of an adequate property model. Consequently, MPU estimation becomes increasingly reliable for an increasing number of data points given a reasonable property model and a specific input domain (the interval of the input variable studied). However, in the usual case of a limited number of data points, MPU estimation is error-prone and requires a thorough analysis of parametric population assumptions. Moreover, the property model under consideration may be inadequate such that systematic deviations of predictions from measurements or benchmark results are observed.[20, 16, 21, 22] Another source of systematic errors are inconsistent data (e.g., outliers).[23, 24, 21]

Here, we study nonparametric bootstrapping introduced by Efron[25], a statistical inference method that meets the challenge of unknown population distributions by sampling from available data. This approach produces an empirical population distribution, where the reference data set itself represents the underlying population. By drawing samples from the empirical population distribution, uncertainty in the parameters can be inferred, which is an essential part in MPU estimation and, hence, for the assessment of transferability.

We apply nonparametric bootstrapping to estimate the uncertainty of predictions of the 57Fe Mössbauer isomer shift (target observable) inferred from a linear property model. The corresponding input variable is the contact electron density at the iron nucleus of a molecular iron compound, which is obtained from calculations based on density functional theory. For 44 molecular iron compounds, we generate input data with 12 exchange–correlation density functionals across Jacob’s ladder (PWLDA, BP86, BLYP, PW91, PBE, M06-L, TPSS, B3LYP, B3PW91, PBE0, M06, TPSSh). We examine (a) systematic errors in both the reference data and the property model, (b) the effect of experimental uncertainty on the model parameters, (c) the reliability of different performance measures (e.g., the squared coefficient of correlation, $r^{2}$, or the root-mean-square error, RMSE), and (d) the dependency of density functional rankings on the number and composition of reference data points.

This paper is organized as follows: In Section 2, we introduce the idea of sampling methods to estimate MPU, explain the essential concepts of nonparametric bootstrapping, discuss its advantages for a specific class of property models (linear least-squares regression), and demonstrate its ability to unravel different sources of systematic errors. In Section 3, we discuss challenges for the prediction of the 57Fe Mössbauer isomer shift on the basis of density functional theory, and reevaluate these challenges in the light of a statistically rigorous calibration analysis.

2 Prediction uncertainty from sampling methods

Since statistically rigorous calibration of physico-chemical property models is the objective of this paper, we first provide a brief and concise review of the relevant concepts and notation needed for this purpose.

To discuss the general concepts of sampling methods, we focus on the common case of a single target observable, $y$, linked to a scalar input variable, $x$, or to a vectorial input variable, $\mathbf{x}$ (in the following, we choose the more general notation $\mathbf{x}=(x_{0},x_{1},...)^{\top}$).

Given a reference data set, $\mathcal{D}$, comprising $N$ data points ($\mathcal{D}\equiv\{(\mathbf{x}_{n},y_{n})\}$ with $n=1,...,N$) we would like to learn predictions of the target observable by calibration of the underlying property model, $f(\mathbf{x},\mathbf{w})\approx y$, where $\mathbf{w}=(w_{0},w_{1},...)^{\top}$ is the vector of parameters.

The mean squared error, MSE, and the coefficient of determination, $R^{2}$, are determined with respect to a reference data set, $\mathcal{D}$, and represent common performance measures of a property model,

[TABLE]

and

[TABLE]

where $\bar{y}$ is the arithmetic mean of the target values, $\{y_{n}\}$, and $\mathbf{w}$ is the vector of parameters. Minimizing the MSE with respect to the parameters, $\partial\text{MSE}_{\mathcal{D},\mathbf{w}}/\partial w_{m}=0\ \forall\ m$, is equivalent to the method of least squares and yields $\min_{\mathbf{w}}(\text{MSE}_{\mathcal{D},\mathbf{w}})$. In the following, we refer to the corresponding parameter vector as $\mathbf{w}_{\mathcal{D}}$, and the shorthand notation $\text{MSE}_{\mathcal{D}}\equiv\min_{\mathbf{w}}(\text{MSE}_{\mathcal{D},\mathbf{w}})$ implies the least-squares objective.

It is important to distinguish the coefficient of determination, $R^{2}$, from the squared coefficient of correlation, $r^{2}$, the latter being independent of the parametric model. Only in some cases (such as linear least-squares regression with a single input variable), both quantities are equivalent.

MSED and $R^{2}$ are established measures of the model performance conditioned on $\mathcal{D}$, respectively. However, in terms of transferability we would like to know the model performance independent of a specific data set. In the following thought experiment, we assume that we can generate an asymptotically large number, $B$, of new (training) samples, $\mathcal{D}_{b}^{\ast}\equiv\{(\mathbf{x}_{n_{b}}^{\ast},y_{n_{b}}^{\ast})\}$ with $n=1,...,N$ and $b=1,...,B$, where the asterisk means “drawn from the population distribution underlying $\mathcal{D}$”. If the data points are randomly drawn from a smooth population distribution, we can safely assume that all samples have no data in common (the importance of independent samples for MPU estimation will be discussed in Section 2.1). For every training sample, $\mathcal{D}_{b}^{\ast}$, we learn the least-squares parameters, $\mathbf{w}_{b}^{\ast}$, of the corresponding property model, $f(\mathbf{x},\mathbf{w}_{b}^{\ast})$, and evaluate its deviation from the target values of the reference sample, $\mathcal{D}$,

[TABLE]

In Eq. (3), the target values $y_{n}$ are elements of $\mathcal{D}$, whereas the predictions $f(\mathbf{x}_{n},\mathbf{w}_{b}^{\ast})$ have been learned from $\mathcal{D}_{b}^{\ast}$.

To estimate the model performance independent of $\mathcal{D}$, $\mathbb{E}[\text{MSE}]$, we average over all training samples, $\mathcal{D}_{b}^{\ast}$,

[TABLE]

where the equality only holds in the artificially ideal case of $f(\mathbf{x}_{n},\mathbf{w}_{b}^{\ast})=f(\mathbf{x}_{n},\mathbf{w}_{\mathcal{D}})\ \forall\ b$. Compared to the $\text{MSE}_{\mathcal{D}}$, the $\mathbb{E}[\text{MSE}]$ additionally incorporates uncertainty in the parameters.

Note that under the assumption of normally distributed data with respect to a parametric model, the $\mathbb{E}[\text{MSE}]$ can be calculated analytically, i.e., sampling is not required in that case. Here, however, we examine the implications of such a critical assumption, which is why we explicitly refrain from considering a parameteric population distribution underlying $\mathcal{D}$.

2.1 Nonparametric bootstrapping — Sampling from available data

The bootstrapping class of sampling methods by Efron[25] has been continuously developed[26, 27]. In nonparametric bootstrapping, the reference sample, $\mathcal{D}$, itself acts as population and, hence, new samples are drawn from $\mathcal{D}$. Consequently, given $N$ input–target pairs in $\mathcal{D}$, these pairs are drawn with equal probability, $p(\mathbf{x}_{n},y_{n})=N^{-1}$. This procedure of sampling from available data is referred to as resampling. The term nonparametric refers to the exclusion of a parameteric population distribution. We drop this term in the following and bootstrapping, if not otherwise mentioned, implies its (original) nonparametric variant.

To generate a bootstrap sample, $N$ elements are drawn from $\mathcal{D}$ with replacement. This procedure is repeated $B$ times, say $B=1000$. There exist two variants of bootstrapping if the reference data set is composed of input–target pairs.[26, 27] The first variant is independent of the property model under consideration and referred to as pair resampling. In that case, the input–target pairs themselves are drawn with replacement, i.e., the bootstrap samples are constructed as $\mathcal{D}_{b}^{\ast}\equiv\{(\mathbf{x}_{n},y_{n})_{b}^{\ast}\}$. The second variant requires a pre-calibrated property model conditioned on the reference data set, $\mathcal{D}$. The resulting residuals, $r_{n}=y_{n}-f(\mathbf{x}_{n},\mathbf{w}_{\mathcal{D}})$, are then subject to resampling, which is referred to as residual resampling. Since independent and identically distributed residuals are assumed in that case, they have the same probability to occur anywhere along the input domain for which data is available. Consequently, the input values, $x_{n}$, are fixed in residual resampling, whereas the residuals, $r_{n}$, are randomly drawn (with replacement) and added to the pre-calibrated property model, $f(\mathbf{x}_{n},\mathbf{w}_{\mathcal{D}})$. Hence, the corresponding bootstrap samples are constructed as $\mathcal{D}_{b}^{\ast}\equiv\{(\mathbf{x}_{n},f(\mathbf{x}_{n},\mathbf{w}_{\mathcal{D}})+r_{n_{b}}^{\ast})\}$. Here, we will exclusively apply pair resampling since calibration prior to bootstrapping imposes further critical assumptions.

Bootstrapping allows us to sample mean and variance of the parameters contained in a calibration model (and of arbitrary other statistics) without relying on parametric population assumptions. Mean and covariance of $\mathbf{w}$ read

[TABLE]

and

[TABLE]

respectively.

The model performance can be estimated from bootstrapping according to Eq. (4), where the square root of $\mathbb{E}[\text{MSE}]$ represents an estimate of the MPU. Note that this bootstrapped variant of the $\mathbb{E}[\text{MSE}]$ is different from that derived in the thought experiment introduced in the beginning of this section. Here, the reference data set, $\mathcal{D}$, and the training data sets, $\mathcal{D}_{b}^{\ast}$, have many data points in common. This overlap promotes underestimation of MPU, and there exist straightforward corrections to the bootstrapped $\mathbb{E}[\text{MSE}]$, for instance, the .632 estimator, [28, 29]

[TABLE]

where $\mathbb{E}[\text{MSE}]_{-n}$ refers to the $\mathbb{E}[\text{MSE}]$ with respect to the $n$-th data point, which is not involved in the calibration of the corresponding property models $f(\mathbf{x},\mathbf{w}_{b}^{\ast})$,

[TABLE]

Here, $\mathcal{B}_{-n}$ and $|\mathcal{B}_{-n}|$ represent the set and number of bootstrap samples not comprising the $n$-th input–target pair, respectively. The constant 0.632 $\approx 1-\text{e}^{-1}$ relates to the probability of a data point to be included in a bootstrap sample.[29] In overfitting situations (when the model is too complex such that data noise is fitted), the .632 estimator can be biased. In such a case, the constant 0.632 requires correction, resulting in the improved .632+ estimator.[30]

Another resampling method that is most popular with respect to MPU estimation is referred to as cross-validation.[29]. In $k$-fold cross-validation, one splits the reference data set into $k$ subsets of most similar size. The model is trained on $k-1$ subsets and validated with respect to the remaining one. This procedure can be performed in $k$ distinct ways. Leave-one-out cross-validation is a frequently applied variant of $k$-fold cross-validation where $k=N$. The corresponding measure of model performance, $\mathbb{E}[\text{MSE}]_{\text{LOO}}$, reads

[TABLE]

where $\mathbf{w}_{\mathcal{D}_{-n}}$ refers to the least-squares estimate of $\mathbf{w}$ with the $n$-th data point removed from the reference data set. Since cross-validation usually deals with a limited number of $k$ training samples, it is, in general, computationally more efficient than bootstrapping. One of the advantages of bootstrapping over cross-validation is, however, the direct assessment of variability for estimated parameters. Furthermore, cross-validation may suffer from a poor bias–variance tradeoff.[29] This issue relates to the observations that (a) for large values of $k$ (in particular for $k=N$), the variance of model parameters is overestimated, while (b) for small values of $k$ (in particular for $k=2$), the expected value of model parameters is biased. Both effects may strongly misestimate the $\mathbb{E}[\text{MSE}]$. In bootstrapping, this issue is resolved by the .632 and .632+ estimators.

We would like to highlight that the sampled parameter distributions resemble empirical variants of posterior distributions employed in Bayesian inference.[31] The latter approach allows for a direct estimation of MPU, which bypasses the need for corrections. For the inference of arbitrary posterior distributions, a Bayesian approach is generally much more involved than bootstrapping. However, if Gaussian posterior distributions are enforced (by choosing Gaussian prior distributions and a Gaussian likelihood function), the MPU can be estimated efficiently through Bayesian inference,[19] which we will discuss in the next subsection. In this study, we will compare the .632 bootstrap estimate of MPU to that obtained from Bayesian inference, assuming normality of the population distribution in the latter case.

Bootstrapping is an appealing alternative to statistical methods implying parametric assumptions on population distributions. On the one hand, modern general-purpose computers can generate and analyze thousands of bootstrap samples in a few seconds given a light fitting problem, i.e., one where calibration is not the limiting step such as in linear least-squares regression. On the other hand, it is an objective approach as it allows for inferring statistical quantities solely from available information (the reference data set). Clearly, if the available data is biased in a way that it badly represents the true underlying population distribution, application of bootstrapping or any other statistical method is not sensible. Major sources of data bias are small sample sizes[26] and gross outliers. For the latter exist established detection methods[29]. Strictly speaking, even bootstrapping builds upon a population assumption, i.e., the reference data set itself being the population. To express it in Chernick’s words[26], bootstrapping does not mean “getting nothing for something”, but “getting the most from the little that is available”.

2.2 Prediction uncertainty of linear regression models

Linear regression refers to a class of calibration procedures, where a target observable, $y$, is estimated to be a linear combination, $f(\mathbf{x},\mathbf{w})$, of $M$ input variables, $x_{m}$ ($m=1,...,M$), and $M+1$ parameters, $w_{m}$ ($m=0,1,...,M$, where the zero-index refers to the intercept, which is quasi-multiplied by $x_{0}=1$),

[TABLE]

where $n=1,...,N$ enumerates the input–target pairs of the reference data set, $\mathcal{D}$.

Given a least-squares objective, the parameter vector $\mathbf{w}_{\mathcal{D}}$ reads

[TABLE]

where the so-called design matrix[19] $\mathbf{X}=(\mathbf{x}_{1},...,\mathbf{x}_{N})^{\top}$ contains all instances of the input vector $\mathbf{x}$ contained in $\mathcal{D}$, and $\mathbf{y}=(y_{1},...,y_{N})^{\top}$ is the vector of target values.

Even though linear least-squares regression is a well-established approach with eligibility for many applications, it implies certain critical assumptions[26]; (a) independent and identically distributed residuals, $r_{n}=y_{n}-f(\mathbf{x}_{n},\mathbf{w}_{\mathcal{D}})$, (b) finite variance of residuals, and (c) variance-free input variables. Only if these assumptions are valid, the least-squares approach yields the best linear unbiased estimate of the regression parameters. For instance, violation of assumption (c) is a ubiquitous phenomenon in the calibration of property models. It is a central topic of this paper and will be discussed in detail in Section 3. Furthermore, if assumption (a) is violated because the variance of the residuals is correlated with one or more input variables, one should instead minimize the weighted MSE, WMSE,[32] for calibration,

[TABLE]

where $\mathbf{U}$ is a diagonal weight matrix (assumption of independent residuals),

[TABLE]

For property models, the elements $u_{n}^{2}$ usually represent experimental variances. The parameter vector obtained from weighted least-squares regression, $\mathbf{w}_{\mathcal{D},\mathbf{U}}$, reads

[TABLE]

This expression is a generalization of the special case covered by Eq. (11) where all elements $u_{n}$ are equal, which is why we drop the subscript $\mathbf{U}$ in the following (it will be evident from the context to which definition of $\mathbf{U}$ we are referring). For the analysis of MPU, it is sufficient to consider the average (root-mean-square) experimental uncertainty, $\langle u\rangle$, in the special case of equality,

[TABLE]

In bootstrapped linear least-squares regression, each bootstrap sample yields a parameter vector

[TABLE]

where the asterisk indicates that $\mathbf{X}$, $\mathbf{y}$ and $\mathbf{U}$ are affected by the resampling procedure. Properly speaking, when explicitly considering uncertainty on the target observable, pair resampling is replaced by triple resampling where bootstrap samples are constructed as $\mathcal{D}_{b}^{\ast}\equiv\{(\mathbf{x}_{n},y_{n},u_{n})_{b}^{\ast}\}$. Locally resolved MPU for a prediction at the input value $\mathbf{x}_{0}$, $u(\mathbf{x}_{0})$, can be estimated from the reduced MSE—by a factor of $N/(N-M-1)$—and the covariance matrix of the model parameters, $\sigma^{2}_{\mathbf{w}^{\ast}}$, introduced in Eq. (6),[16]

[TABLE]

Note that $u^{2}(\mathbf{x}_{0})$ estimates the prediction variance of the target observable at the input value $\mathbf{x}_{0}$ with respect to the average experimental variance, $\langle u\rangle^{2}$, which would converge to the reduced MSE if all systematic errors have been removed. For a new series of measurements at input value $\mathbf{x}_{0}$ with experimental variance $u_{0}^{2}$, it is necessary to add the difference $u_{0}^{2}\,-\,\langle u\rangle^{2}$ to the reduced MSE in Eq. (17). The relation between the $\mathbb{E}[\text{MSE}]_{.632}$ and the locally resolved MPU is given by

[TABLE]

where the approximation sign arises, i.a., from the assumption of normality of the parameter distributions as indicated by the second term on the right-hand side of Eq. (17). The entire right-hand side of Eq. (17) can be replaced by its sampled analog (bootstrapping of prediction intervals[27]), but for the sake of user-friendliness, we decided to employ a limited number of distinct characteristic values (in this study, two for a prediction, $\bar{\mathbf{w}}^{\ast}$, and three for the corresponding uncertainty, $\sigma^{2}_{\mathbf{w}^{\ast}}$) to facilitate comparisons of different computational methods (e.g., density functionals) with respect to their performance.

Bayesian linear regression is an alternative to bootstrapped linear least-squares regression. If there is good reason to assume that the parameter distributions are Gaussian, Bayesian linear regression can be much more efficient than sampling-based regression. Bayesian linear regression provides analytical posterior distributions of parameters,[19] the maxima of which represent the best-fit parameter vector, $\mathbf{w}_{\text{MAP}}$ (MAP, maximum posterior). The procedure is outlined in the Appendix and yields a measure of model performance referred to as mean predictive variance (MPV),

[TABLE]

where $s^{2}(\mathbf{x}_{n})$ is the analogue of $u^{2}(\mathbf{x}_{n})$ obtained from bootstrapping. Alternatively, the MPV can be obtained by summing over a dense grid of input values (with a number of grid points $\gg N$),[16] which yields a smoother result in the sense that it is less dependent on the particular choice of reference input values. In that case, it is important to specify the bounds of the input interval and the spacing between or distribution of grid points.

In the following, we will report the square roots (RMSE, R632, RMPV, RLOO) of the performance measures introduced ($\text{MSE}_{\mathcal{D}}$, $\mathbb{E}[\text{MSE}]_{.632}$, MPV, $\mathbb{E}[\text{MSE}]_{\text{LOO}}$) for the sake of better comparability with the experimental uncertainty.

2.3 Jackknife-after-bootstrapping — Data diagnostics

So far, the discussion of statistical inference was built on the implicit assumption of data sets representing their underlying population distributions well. Verifying the validity of this assumption is a tedious task, but a diagnostic referred to as jackknife-after-bootstrapping provides a good approximation to the problem[26, 33]. In the first step of that diagnostic (given a data set with $N$ input–target pairs), the sampled ensemble of parameters, $\{\mathbf{w}_{b}^{\ast}\}$, is decomposed into $N$ different sets, $\{\mathbf{w}_{b}^{\ast}\}_{-n}$ ($n=1,...,N$), the $n$-th of them containing only those parameters $\mathbf{w}_{b}^{\ast}$ learned from bootstrap samples in which the $n$-th input–target pair is not included. For instance, one obtains a bootstrap estimate of the mean of a parameter with the $n$-th data point removed, $\bar{\mathbf{w}}_{-n}^{\ast}$. If this jackknifed mean deviates significantly from that inferred from the complete reference data set, $\bar{\mathbf{w}}^{\ast}$, we have an indication that the $n$-th data point biases the calibration.

An appealing feature of the jackknife-after-bootstrapping method is its efficiency. Instead of running $N$ extra bootstrap simulations, one only performs a single bootstrap simulation on the complete reference data set. The reason is that some bootstrap samples do not contain certain data points. In this way, the decomposition of parameters into $N$ subsets can be performed simultaneously to the bootstrap simulation on the complete reference data set. Since the probability of a data point to be excluded from a bootstrap sample is roughly[29] 0.368 $\approx\text{e}^{-1}$, we recommend to increase the default number of bootstrap samples by a factor of approximately 3 to preserve the intended calibration accuracy.

3 Case study: Isomer-shift calibration in theoretical ${}^{\text{57}}$Fe Mössbauer spectroscopy

Calibration has been frequently applied to predict the isomer shift observed in 57Fe Mössbauer spectroscopy.[34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60] The corresponding theory[61] postulates a linear relationship between the measurable isomer shift, $\delta_{\text{exp}}$, and the difference in the contact electron density (CED), $\rho_{\text{absorber}}-\rho_{\text{source}}$, of the 57Fe isotope embedded in two different chemical environments (referred to as absorber and source),

[TABLE]

$g(r)$ is a function of the average charge radius, $r$, of an iron nucleus, and $\Delta r$ is the difference between the charge radii of the excited state and the ground state of an iron nucleus. In the corresponding property model, $\delta(\rho_{\text{absorber}},\mathbf{w})$ with $\mathbf{w}=(w_{0}\ w_{1})^{\top}$, all quantities of the right-hand side of Eq. (20) except for $\rho_{\text{absorber}}$ are hidden in the regression parameters, $w_{0}$ (intercept) and $w_{1}$ (slope),

[TABLE]

The absorber CED (ACED), $\rho_{\text{absorber}}$, is determined from an electronic structure method, typically based on density functional theory. Since every specification employed in an electronic structure calculation (e.g., density functionals, basis sets, integration grids, convergence criteria) may affect the value of $\rho_{\text{absorber}}$, the property model needs to be calibrated every time a specification is changed. As noted in Section 2.2, one of the key assumptions in applying linear least-squares regression are variance-free input values. While electronic structure calculations yield virtually variance-free results (neglecting numerical errors and convergence threshold effects), the expectation value of an observable remains unpredictable due to method-inherent systematic errors, which are collectively referred to as method inadequacy.[20, 16] Method inadequacy causes the average uncertainty of the reference isomer shifts to be non-reproducible, no matter which values we choose for $w_{0}$ and $w_{1}$. Consequently, one should have serious doubt on the validity of the free-of-variance assumption.

Another frequent source of systematic errors are inconsistent data. On the one hand, the reference isomer shifts employed in this study have been recorded at different temperatures (4–100 K), which can lead to signed deviations of about $-$0.02 mm s*-1* (second-order Doppler shift).[61] Bochevarov, Friesner, and Lippard[50] proposed to consider only those isomer shifts recorded at liquid helium temperature (4.2 K). This situation is clearly desirable, but it would have limited our reference data set in terms of chemical diversity. On the other hand, it is difficult to ensure that the molecular-structure representations employed are sufficiently accurate for reliable ACED calculations. Even though crystal structures guide the search for the correct minimum on the Born–Oppenheimer surface, there is no guarantee that structure optimization yields reliable results. Furthermore, not only the iron-containing compound itself may be important for isomer-shift calibration, but also the closer environment of the solid sample such as adjacent iron complexes, counter-ions, or solvent molecules. The effect of the molecular-structure representation on the ACED remains an issue to be studied, which is beyond the scope of this work.

In essentially all previous calibration studies of the isomer shift,[34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60] the squared coefficient of correlation, $r^{2}$, served as a measure to assess the performance of an electronic structure method (note that in linear least-squares regression with a single input variable, as applied here, $r^{2}$ is identical to the coefficient of determination, $R^{2}$). However, according to Eqs. (2) and (4), there is no guarantee that $R^{2}$ ($=r^{2}$) allows for a reliable comparison of two electronic structure methods with respect to their transferability, because this performance measure does not take into account uncertainty in the model parameters.

The incompleteness of $r^{2}$ (or the RMSE) as model performance measure, the possible existence of inconsistent reference data, and the unpredictable variability in the input variable has motivated us to reexamine isomer-shift calibration in the light of a statistically rigorous analysis. We assess the reliability of different performance measures and study the transferability of 12 density functionals with respect to isomer-shift predictions for 44 iron compounds of considerable chemical diversity (cf. Table 1).

Details on the computational protocol employed for both statistical calibration analysis and quantum chemical calculations are provided in the Appendix. Our statistical calibration program reBoot developed for the analysis presented in this paper will be made available on our webpage[62]. In combination with the data provided in the Supporting Information, reBoot allows one to reproduce all results of our paper. Furthermore, reBoot can be harnessed to apply the statistical calibration methods presented in our paper to arbitrary polynomial property models that are linear with respect to their parameters. Note that the statistical methods implemented in reBoot are not limited to this family of models. For instance, implementation of non-polynomial models or models being nonlinear in their parameters would be straightforward.

3.1 Reference set of molecular iron compounds

In previous calibration studies of the isomer shift[34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], a variety of iron complexes was considered. Here, we chose a diverse set of molecular iron compounds (Table 1 and Figs. S6–S16) representing wide ranges of formal oxidation states, spin multiplicities, total charges, and ligand environments (type, number, and spatial orientation). We blended parts of previous reference sets by Friesner and Lippard[50], Neese[49], Noodleman[41], and our group[54] with iron(I) complexes[63, 64, 65, 66, 67, 68, 69, 70]. We explicitly excluded linear and T-shaped iron compounds with a formal oxidation state of +1, since these species are known to reveal pronounced spin–orbit coupling[71, 72, 73, 74, 75] not considered in our computational approach. While the effect of strong spin–orbit coupling on iron CEDs remains a subject to be studied in more detail, scalar-relativistic effects have been found to induce only a constant shift of nonrelativistic CEDs for iron-containing molecules, which is why isomer-shift calibration is frequently based on nonrelativistic calculations.[61] In our previous study on 57Fe isomer-shift calibration [54], we found that scalar-relativistic calculations lead to a slightly higher correlation between the target observable and the input variable compared to nonrelativistic calculations. However, in preparation of the present work, we detected a specification error for iodine in the EMSL basis-set database[76] for the def2-TZVP basis set and the Molcas computer program that turned out to be the source of the slightly better performance. Correction of this error reveals equivalent input–target correlation for scalar-relativistic and nonrelativistic results.

[FIGURE:]

[FIGURE:]

3.2 Effect of experimental uncertainty on model parameters

When applying linear least-squares regression to all reference isomer shifts ($N=44$, see Fig. 1, left), we find an RMSE of 0.07 mm s*-1* for all hybrid density functionals (B3LYP, B3PW91, PBE0, M06, TPSSh). For the pure density functionals (PWLDA, BP86, BLYP, PW91, PBE, M06-L, TPSS), the RMSE ranges from 0.08 mm s*-1* (M06-L, TPSS) to 0.10 mm s*-1* (BLYP); see also Table 2.

Even though the RMSE is a lower bound to the MPU inferred from this specific set of reference isomer shifts, it is already significantly larger than the average experimental uncertainty of $\langle u\rangle=0.02$ mm s*-1* found for molecular iron compounds.[87] However, it is possible that $\langle u\rangle$ is larger than 0.02 mm s*-1* for the compounds studied here, since several isomer-shift measurements were reported without uncertainty (cf. Table 1). Selecting only those measurements for which uncertainty has been reported ($N=30$), we also find an average experimental uncertainty of $\langle u\rangle=0.02$ mm s*-1*. The RMSE still ranges from 0.06 mm s*-1* (B3LYP, B3PW91, PBE0, M06) to 0.10 mm s*-1* (BLYP). This discrepancy between the average experimental uncertainty and the RMSE indicates that explicit consideration of experimental uncertainty (through weighted least-squares regression) may have a minor effect on the model parameters. This hypothesis can be examined by iteratively reweighted linear least-squares regression [16, 32] (see the Appendix for more details). Table 3 summarizes the results for B3LYP, which was repeatedly found to perform superior to other density functionals in isomer-shift calibration.[35, 36, 42, 49, 50, 55]

While the slope, $w_{1}$, is almost invariant to the approaches applied, the intercept, $w_{0}$, changes by 57.41 mm s*-1* when switching from ordinary (unweighted) least-squares regression to weighted least-squares regression, whereas least-squares regression and iteratively reweighted least-squares regression yield almost identical results (deviation of 3.65 mm s*-1*). In the first instance, these results suggest that experimental uncertainty does not need to be taken into account since method inadequacy and/or data inconsistency appear to play a dominant role. However, the bootstrapped standard deviation of the intercept, which is larger than 100 mm s*-1* in each case, indicates that the differences between all calibration procedures are insignificant. Qualitatively analogous results were obtained for the remaining density functionals. As a consequence, we assume that calibration of our property model, $\delta(\rho_{\text{absorber}},\mathbf{w})$, is not perturbed by the inclusion of reference isomer shifts for which experimental uncertainty was not reported. If not mentioned otherwise, the application of bootstrapping will imply ordinary (unweighted) least-squares regression, which is equivalent to assigning the average experimental uncertainty of $\langle u\rangle=0.02$ mm s*-1* to all reference isomer shifts.

In the next two subsections, we will explore the possibilities to unravel the effects of method inadequacy and data inconsistency on the model residuals.

3.3 Selection of the property model with Occam’s razor

First of all, we examine whether systematic method inadequacy is present. Systematic method inadequacy would result in residuals which show a trend with respect to the underlying property model instead of random scatter. While Gaussian process regression[104] is a reliable approach to infer the model complexity with the highest transferability, we will apply a simple alternative based on Occam’s razor. In addition to the linear model, which is based on a physical theory[61], we choose quadratic, cubic, and quartic models,

[TABLE]

with $M=2,3,4$. In bootstrapping procedures, the model parameters are calibrated with respect to a bootstrap sample and then validated at a reference sample. Therefore, if the model is too rigid or too flexible (i.e., when it features too few or too many parameters, respectively), we will observe decreased transferability compared to a balanced model through an increase of the MPU. Table 4 summarizes the results for B3LYP.

The linear and quadratic models reveal the lowest MPU as measured by the .632 estimator. According to Occam’s razor, we choose the simpler of both models, which is also the only one built on physical grounds[61]. By inspection of Fig. 1, the residuals of the linear model appear randomly distributed. Qualitatively identical results were obtained for the remaining density functionals. Consequently, we may assume that the discrepancy between RMSE and $\langle u\rangle$ is rooted in inconsistent data and/or random method inadequacy. The latter effect would be a consequence of non-constant systematic errors introduced by the electronic structure method under consideration, which cause an increase in data noise.[16]

3.4 Assessment of data inconsistency based on jackknife-after-bootstrapping

In Fig. 1 (left), uncertainty in the regression parameters is represented by the 95% confidence band (turquoise area, see the Supporting Information for details) obtained from bootstrapping ($B=10^{5}$). The corresponding ACED values were obtained from B3LYP calculations. The black line represents the mean of regression parameters over all bootstrap samples, $\bar{\mathbf{w}}^{\ast}$, as defined in Eq. (5). R632, RMPV, and RLOO equal the RMSE (0.07 mm s*-1*); see also Table 2. Note that the experimental resolution is limited to 0.01 mm s*-1*. If instead it would be artificially increased to 0.001 mm s*-1*, which can be approximated by adding a trailing zero to the reference isomer shifts, we obtain RMSE $=$ 0.066 mm s*-1*, R632 $=$ 0.070 mm s*-1*, RMPV $=$ 0.069 mm s*-1*, and RMPV $=$ 0.070 mm s*-1*. For details on a statistically valid increase of the experimental resolution, see the Appendix. Hence, the low experimental resolution masks the effect of parameter uncertainty on the MPU as measured by R632, RMPV, and RLOO, at least for this specific composition and number of isomer shifts. Moreover, the RLOO appears to be an efficient alternative to the R632 for isomer-shift calibration. Qualitatively similar results were also obtained for the remaining density functionals (Table 2).

In Fig. 2 (top left and top right), histograms of the bootstrapped parameter distributions ($B=10^{5}$) are shown. The solid curves are Gaussian posterior distributions obtained from Bayesian linear regression. We find systematic deviation of the posterior slope distribution from its bootstrapped variant, in particular with respect to higher moments (skewness), which may explain the slight difference between R632 and RMPV when artificially increasing the experimental resolution. The skewness of the bootstrapped slope distributions may arise from data points with an above-average effect on the objective (here, sum of least squares) such as those at the boundaries of the input domain or apparent outliers. By contrast, the agreement between the posterior intercept distribution and its bootstrapped variant appears to be sufficiently high. Q–Q plots of bootstrapped and Gaussian parameter distributions (Fig. S1) support these findings.

To examine the reliability of the sampled parameter distributions and the MPU associated with them, we apply the jackknife-after-bootstrapping method. An overview in Fig. 3 (left) of the RMPV versus the normalized unsigned deviation in the intercept, $\Delta w_{0,n}$,

[TABLE]

optimized with respect to both $\mathcal{D}$ and $\mathcal{D}_{-n}$, reveals no clear correlation between the quantities (here, $\mathcal{D}_{-n}$ refers to the reference data set with the $n$-th data point removed, and the RMPV is given with respect to $\mathcal{D}_{-n}$). The same holds true for the correlation between the RMPV and the normalized unsigned deviation in the slope (Fig. 3, center),

[TABLE]

As the parameters are correlated in the property model, we also calculated the normalized root-mean-square deviation of the property model, $\Delta\delta_{n}$, calibrated with respect to both $\mathcal{D}$ and $\mathcal{D}_{-n}$,

[TABLE]

A plot (Fig. 3, right) of the RMPV versus $\Delta\delta_{n}$ reveals a more distinct correlation. At smaller deviations, the RMPV is almost constant, whereas it decreases overlinearly at larger deviations. Hence, if the removal of a data point significantly changes the functional form of the property model, it also has a significant effect on the MPU.

The plots of Fig. 3 reveal four data points (#1, #2, #13, #28, cf. Table 1 and Fig. 1, left) that lead to a distinct decrease of the MPU and high values of $\Delta\delta_{n}$, respectively. The question arises whether these data points are inconsistent in the sense that if they are present, the underlying population distribution may not be well-represented by the data set under consideration. Since we are studying several density functionals, we can examine whether we find the same potentially inconsistent data points in the remaining cases. We find that three of the four potentially inconsistent data points identified with B3LYP are also identified as most likely inconsistent with all other density functionals (we considered the four highest values of $\Delta\delta_{n}$ for each density functional). This finding indicates that these data points (#1, #2, #28) are not affected by inconsistent ACED calculations as they have been identified to be inconsistent irrespective of the density functional employed. Rather, systematic measurement errors or deficiencies in the molecular-structure representation may be responsible for their inconsistent status. Moreover, the fourth potentially inconsistent data point identified with B3LYP (#13) has also been identified with three other hybrid density functionals (B3PW91, PBE0, M06) and is ranked fifth with respect to $\Delta\delta_{n}$ for the remaining hybrid density functional (TPSSh). In this case, we interpret data point #13 as affected by inconsistent ACED calculations (method inadequacy) as it was only identified inconsistent by a particular category of density functionals. Likewise, data point #7 (cf. Table 1 and Fig. 1, left) is one of the first four data points identified potentially inconsistent by all pure density functionals (PWLDA, BP86, BLYP, PW91, PBE, M06-L, TPSS) and one hybrid density functional (TPSSh), and is ranked fifth with respect to $\Delta\delta_{n}$ for all remaining hybrid density functionals but M06. We do not find obvious similarities in the corresponding compounds regarding the different categories in Table 1, but in all five cases, iron is coordinated to nitrogen (Fig. 4). We also compared the calculated versus ideal expectation value of the $\langle S^{2}\rangle$ operator for all open-shell complexes (Table S5), but find no anomalies for the inconsistent data points. Likewise, we cannot confirm that the root-mean-square deviation of atomic positions is particularly high for those iron complexes corresponding to inconsistent input–target pairs (Table S1).

To obtain reliable estimates of the MPU, we decided to remove all inconsistent data points (#1, #2, #7, #13, #28, cf. Table 1 and Fig. 1, left). Repeating bootstrapped ($B=10^{5}$) linear least-squares regression (Fig. 1, right), the updated property model (solid line) clearly deviates from that one calibrated with respect to the reference data set including the inconsistent input–target pairs (dashed line). Fig. 2 (bottom left and bottom right) reveals that mean and standard deviation of the bootstrapped parameter distributions (histograms) change significantly, the latter decreases by a factor of 2 (cf. Table 5). The corresponding Gaussian posterior distributions obtained from Bayesian linear regression (solid curves) are now both quite similar to their sampled counterparts. This finding is also supported by Q–Q plots of bootstrapped versus Gaussian parameter distributions (Fig. S1). Even though the 95% confidence band (Fig. 1, right) is narrower than before, we now find a difference between the RMSE (0.03 mm s*-1*) and R632 as well as RMPV and RLOO (0.04 mm s*-1*, respectively). This finding is an artifact resulting from the low experimental resolution of 0.01 mm s*-1*. Increasing the experimental resolution artificially from 0.01 mm s*-1* to 0.001 mm s*-1*, we obtain RMSE $=$ 0.035 mm s*-1*, R632 $=$ 0.036 mm s*-1*, RMPV $=$ 0.036 mm s*-1*, and RLOO $=$ 0.036 mm s*-1*. In this case, we see that the difference between RMSE and R632 decreases from 0.004 mm s*-1* to 0.001 mm s*-1*, which we would expect for a narrowing of the confidence band. Furthermore, due to the higher similarity of the bootstrapped and Gaussian posterior distributions after removal of the inconsistent data points, we also find that the difference between R632 and RMPV vanishes (before: 0.001 mm s*-1*) for the increased experimental resolution. This result indicates that the normal-population assumption appears to be reasonable for isomer-shift calibration after removal of inconsistent data points. Again, the RLOO appears to be an efficient alternative to the R632 for isomer-shift calibration. Note that the updated RMSE (0.03 mm s*-1* instead of 0.07 mm s*-1*) is now much closer to the average experimental uncertainty ($\langle u\rangle=0.02$ mm s*-1*). For the remaining density functionals, the RMSE ranges from 0.03 mm s*-1* (B3PW91, PBE0) to 0.08 mm s*-1* (BLYP), while R632, RMPV, and RLOO range from 0.04 mm s*-1* (B3PW91, PBE0, TPSSh) to 0.08 mm s*-1* (BLYP); see also Table 6.

Noteworthy, while inclusion of the inconsistent data points in the reference data set ($N=44$) leads to overestimation of the MPU as measured by the R632 (86% of the data points lie within the 68% prediction band obtained from $u(\mathbf{x})$ defined in Eq. (17)), their exclusion ($N=39$) results in a significant improvement of MPU estimation (72% of the data points lie within the 68% prediction band); see also Fig. 5. Note that the prediction bands shown in Fig. 5 (bottom) can be directly related to $u(\mathbf{x})$ defined in Eq. (17).

Because of the decrease of the RMSE, we examined once again the importance of explicitly considering experimental uncertainty. The results are summarized in Table 7 and lead us again to the conclusion that calibration of $\delta(\rho_{\text{absorber}},\mathbf{w})$ is not perturbed by the inclusion of reference isomer shifts for which experimental uncertainty has not been reported. Furthermore, we applied Occam’s razor to our consistent reference data set ($N=39$). In all cases ($M=1,2,3,4$), we obtained an R632 of 0.04 mm s*-1*, which confirms once again the validity of the linear (physical) model.

Combining all these findings, we are confident that the updated, consistent reference data set represents the underlying population distribution sufficiently well if hybrid density functionals generate the input data. Hence, our reference set of $N=39$, already pruned by statistically critical outliers, provides a well-defined starting ground for further parametrization studies or for systematic extensions of the reference data. Hence, we may call this special set of data the MIS39 data set (for 39 reference data points of Mössbauer isomer shifts).

3.5 How reliable are density functional rankings based on a specific data set?

In the previous subsection, we found for the particular composition and number of data points that the effect of parameter uncertainty on the MPU is very small ($<$0.01 mm s*-1* for all density functionals). Hence, given this specific reference data set, we are confident that the RMSE (and with it the squared coefficient of correlation, $r^{2}$) is a good approximation to performance measures such as R632, RMPV, and RLOO, and therefore, suited to set up a ranking of density functionals. In the following, we only distinguish between first and other places in a ranking as only those density functionals will be considered for actual applications that reveal highest transferability. Regarding our reference set of isomer shifts ($N=39$), B3LYP, B3PW91, and PBE0 are placed first in all rankings studied (RMSE, R632, RMPV, RLOO), with TPSSh being placed first in the rankings based on R632, RMPV, and RLOO (cf. Table 6). This finding is consistent with those of other calibration studies of the isomer shift.[35, 36, 42, 49, 50, 55]

However, in practice, it is relevant to know which density functional yields the most accurate predictions independent of the reference data employed for calibration of the property model. So far, we considered rankings of density functionals conditioned on a specific reference data set, but what we aim at is an unconditional ranking of density functionals. Otherwise, we cannot assess the transferability of a property model trained on a specific density functional to data not involved in its calibration. Clearly, the dependency of any statistical measure on a specific choice of reference data cannot entirely be removed, but bootstrapping can yield a good approximation to the problem given the data under consideration is representative of the underlying population distribution, which we have assessed in detail by the jackknife-after-bootstrapping method.

To study the reliability of density functional rankings with respect to the composition and number of data points, we have drawn samples of different size ($N=5,10,20,39$) with replacement ($B=2.5\times 10^{3}$, respectively) from the empirical population distribution of the consistent reference data set ($N=39$). For every bootstrap sample, the isomer shifts have been perturbed randomly according to their experimental uncertainty (for details, see the Appendix) to allow for statistically justifiable variation between the synthetic data sets. In Fig. 6, we show the percentage of first places that a density functional has reached for 2500 different data sets and four different data-set sizes. The rankings are based on the RMSE.

For all different data sets considered, B3LYP, B3PW91, and PBE0 are still most frequently placed first, but the respective percentage varies significantly. For $N=5$ data points, PBE0 reveals the highest number of first places, but only in about 60% of the cases. Even all pure density functionals are placed first in at least 10% of the cases, although the different performance measures are clearly in favor of hybrid density functionals regarding the reference data set (cf. Table 6). Hence, calibration studies based on 5 data points are highly susceptible to random conclusions about the transferability of density functionals. With an increasing number of data points, the percentage of first places continuously decreases for all pure density functionals and for M06, while TPSSh reveals relatively constant results for all data points considered (between 30% and 40% first places). By contrast, the percentage of first places for B3LYP, B3PW91, and PBE0 increases continuously, with PBE0 being ranked first in $>95\%$ of the cases. Qualitatively identical and quantitatively similar results were obtained for density functional rankings based on the RMPV (cf. Fig. S2). This finding indicates that for the given experimental resolution, parameter uncertainty plays no significant role, even for 5 data points where the standard deviation in the parameters is about 4 times larger compared to 39 data points, the number of reference data points considered in this study (Fig. S5). Consequently, both RMSE and $r^{2}$ can be considered stable performance measures for density functional rankings applied in isomer-shift predictions of molecular iron compounds. Note that we did not employ the R632 as performance measure as it is very costly to sample an error for every bootstrap sample, which is equivalent to bootstrapping bootstrap samples (double bootstrapping[27]). However, as shown above, the RMPV can be expected a good approximation to the R632 (as well as RLOO).

The fact that more than one density functional is placed first on average is clearly an indicator of the low experimental resolution. When increasing the experimental resolution from 0.01 mm s*-1* to 0.001 mm s*-1* in a statistically sound way (for details, see the Appendix), the results based on the RMSE change significantly (Fig. 7). The dispersion of the rankings is much more pronounced for 5 data points (all density functionals are placed first in less than 40% of the cases), and a clear preference for a density functional (more than 95% first places) can be expected only for a number of data points significantly larger than that employed in this study ($N=39$). Similar results were obtained for density functional rankings based on the RMPV (cf. Fig. S3).

Reducing the data-set dependency of the MPU by the application of bootstrapping and increasing the experimental resolution in a statistically sound way, we find a clear preference for PBE0 as the density functional with the highest transferability. This preference even remains (a) when resubstituting the 5 inconsistent input–target pairs (#1, #2, #7, #13, #28, cf. Table 1 and Fig. 1, left) into our consistent reference data set ($N=44$), and (b) when selecting only those reference isomer shifts for which experimental uncertainties were reported (without resubstitution of inconsistent data points, $N=28$); cf. Fig. S4. Therefore, we suggest that PBE0 is the density functional of choice for applications (with the CP(PPP) basis for Fe and def2-TZVP for all other elements; see Appendix), which leads to the following linear isomer-shift model obtained from bootstrapped ($B=10^{6}$) linear least-squares regression,

[TABLE]

the corresponding covariance matrix,

[TABLE]

and the corresponding reduced MSE,

[TABLE]

We omitted the subscript “absorber” for the ACED, $\rho$. For reasons of reproducibility, we specify four significant figures for the characteristic values contained in $\bar{\mathbf{w}}_{\text{PBE0}}^{\ast}$ and $\sigma_{{\mathbf{w}}_{\text{PBE0}}^{\ast}}$, and for the reduced MSE. For the mean ACED, $\bar{\rho}$, we employed the raw-data precision provided in Tables S2 and S3. The standard deviations in parentheses were obtained from bootstrapping ($B=10^{3}$) the $10^{6}$ parameter estimates.

Given a single new ACED, $\rho_{0}$, the corresponding prediction uncertainty, $u(\boldsymbol{\rho}_{0})$ with $\boldsymbol{\rho}_{0}=(1,\rho_{0}-\bar{\rho})^{\top}$, would, according to Eq. (17), read

[TABLE]

3.6 Effect of exact exchange on the MPU

In this last section, we briefly discuss our observation why hybrid density functionals yielded significantly lower MPU estimates compared to pure density functionals. A systematic (even linear) behavior of relative energies of states of different spin multiplicity on the admixture of exact exchange (i.e., Hartree–Fock-type exchange) in an energy density functional has already been observed more than 15 years ago[105, 106, 107]. Since then, the exact exchange admixture, measured by the linear parameter $c_{3}$ following the notation in Ref.  105, has been well recognized as one of the most crucial parameters determining the accuracy of (hybrid) density functionals (see Refs.  108, 109, 110, 111 for recent systematic investigations). The linear dependence of such relative energies in many (but not all) cases is remarkable and results from a linear dependence of the absolute electronic energies on $c_{3}$, in which the effect of a self-consistently optimized (and hence, changing) electron density plays a negligible role.

Here, we calculated the ACED of all 44 reference compounds (cf. Table 1) for varying amounts of exact exchange from $c_{3}=0.00$ to $c_{3}=0.50$ in steps of 0.05 in the B3LYP density functional (for which $c_{3}=0.20$ was originally set). Interestingly, also the contact density features a linear dependence on $c_{3}$ ($r^{2}>0.9995$ in all compounds), although with different slopes. In other words, the first derivative of the ACED with respect to $c_{3}$ is constant for a given compound (Fig. 8), but different for all complexes. This complex-specificity is also confirmed by bootstrapping ($B=10^{3}$) the uncertainty of the first derivative of the ACED (see the error bars in Fig. 8, which represent three standard deviations).

Consequently, a change in exact exchange admixture does not lead to a unique ACED shift for all compounds, which is why the MPU is a nonlinear function of $c_{3}$. Plotting the R632 against the exact exchange admixture reveals that this performance measure is minimized for $c_{3}=0.20$–$0.25$ (cf. Fig. 9 where $N=39$), which indicates why B3LYP, B3PW91 (both featuring $c_{3}=0.20$), and PBE0 ($c_{3}=0.25$) are the most transferable density functionals with respect to isomer-shift calibration.

The observation of linear dependence of the ACED is not straightforward to explain. It is clear that only basis functions with angular momentum quantum number $l=0$ ($s$-type functions) can contribute to the nonrelativistic contact density in an atomic system (note the short-range behavior $\propto r^{l}$ of the radial function in a one-electron atom[112]). This short-range behavior is (also in the spherically averaged case of an atom in a molecule[113]), determined by the lowest-order contributions of the potential to a Taylor series expansion in terms of the radial distance in the short-range quantum mechanical differential equations for a single electron. The lowest-order constributions are the centrifugal and the nuclear point-charge potentials, but not the electron–electron interaction potentials, which provide a constant contribution at the nuclear position[114]. Hence, for a spherically symmetric atom, this is also true for the Hartree–Fock exchange interaction, which contributes a constant term to the potential at the nuclear origin[114] and would therefore not, in this case, affect the radial function and hence not the electron density at the nucleus.

For a (nonspherical) molecular system, additional contributions to the contact density need to be considered. Following the partitioning scheme of the ACED by Neese,[35] we find for compound #43 (iron pentacarbonyl) that the iron $s$-functions in the CP(PPP) basis set contribute dominantly to the ACED ($>99.9\%$). Hence, their contribution to the ACED alone already reveals the linear trend of the ACED with respect to exact exchange admixture. However, the three $s$-functions in the CP(PPP) basis set with the largest exponent are constant with respect to a change of $c_{3}$, whereas the smallest 11 $s$-functions change linearly with $c_{3}$, even though their slope is sometimes positive (for larger exponents) and sometimes negative (for smaller exponents), which makes a detailed analysis difficult. Interestingly, considering not only the charge at one specific point, i.e., the ACED, but the integral of the electron density over some space, i.e., a partial charge of the iron atom in iron pentacarbonyl, also reveals a linear behavior with $c_{3}$ although with negative slope (see also Ref.  108). Hence, Mulliken and Löwdin partial charges decrease with increasing $c_{3}$, while the ACED increases.

4 Conclusions and outlook

To reliably estimate the prediction uncertainty of a property model, it is important to identify and correct for systematic errors due to, e.g., inconsistent measurements, parametric population assumptions, or inadequate computational methods. Here, we studied this issue at the example of 57Fe Mössbauer isomer-shift predictions based on a linear model. Twelve density functionals across Jacob’s ladder were considered for the calculation of the ACED for 44 chemically diverse molecular iron compounds (formal oxidation states: 0, +1, +2, +3, +4), whereas the corresponding target data refer to measured isomer shifts. We explicitly considered uncertainty in the model parameters, which may be a crucial ingredient for the estimation of MPU. For this purpose, we employed both bootstrapping[25, 27, 26, 29] and Bayesian linear regression based on the evidence approximation.[19].

First of all, we found that the RMSE, which measures the standard deviation of the model residuals, is significantly larger than the average experimental uncertainty (0.07–0.10 mm s*-1* versus 0.02 mm s*-1*). This discrepancy cannot be explained by the simplicity of the linear model as more complex property models (quadratic, cubic, quartic) were found to yield either equal or larger MPU estimates when applying bootstrapping. However, with the jackknife-after-bootstrapping approach, which probes the sensitivity in the MPU with respect to small changes in the reference data set, we could identify 5 inconsistent data points. The reliability of the jacknife-after-bootstrapping approach can be assessed by determining the fraction of data points lying in a certain prediction interval. When including the inconsistent data points, 86% of the reference isomer shifts lie in the 68% prediction band (representing $u(\mathbf{x})$ as defined in Eq. (17)), whereas excluding them leads to only 72% of the reference isomer shifts contained in the 68% prediction band, which reveals a significant improvement of the statistical reliability. Furthermore, the high agreement of bootstrapped and Gaussian posterior parameter distributions (the latter were obtained from Bayesian linear regression) suggests that the normal-population assumption is reasonable for isomer-shift calibration given the reference data set is carefully selected.

The new (consistent) reference data set still leads to overestimation of the average experimental uncertainty (RMSE of 0.03–0.08 mm s*-1*), which can be assigned to non-constant systematic errors of the density functionals leading to random shifts of the (unknown) true input values. This random method inadequacy suggests that weighted least-squares regression, where experimental uncertainties are explicitly introduced as a weight, may bias the calibration of the property model. However, on the basis of bootstrapped regression we did not find a significant difference between the ordinary (unweighted) and the weighted setup due to pronounced parameter uncertainty, whereas Grandjean and Long[59] recommend the application of weighted least-squares regression for the calibration of isomer-shift models.

Comparison of the RMSE with more reliable performance measures, which consider uncertainty in the parameters (such as R632, RMPV, or RLOO), reveals nearly equal results in most of the cases. This finding can be explained by the low experimental resolution of 0.01 mm s*-1* compared to the width of the input domain (here, 1.08 mm s*-1*), which almost completely masks parameter uncertainty. Therefore, simple performance measures such as the RMSE and the squared coefficient of correlation, $r^{2}$, are expected to be reliable for the construction of density functional rankings, where relative rather than absolute MPU estimates are required. The same holds for the root-mean-square deviation (RMSD) and the mean unsigned error (MUE) as defined and discussed in a calibration study on physical properties of crystals by Pernot et al.[16]

Next, we examined the sensitivity of density functional rankings with respect to the composition and number of data points. For this purpose, we generated 10000 synthetic data sets (based on bootstrapping) with 5, 10, 20, and 39 (the size of our reference data set) data points. For 5 data points, false conclusions about the performance of a density functional are likely. For instance, the “best” density functional, PBE0, is placed first in only about 60% of the cases. In the case of 39 data points, B3LYP and B3PW91 are placed first in 85–90% of the cases, whereas PBE0 is placed first in more than 95% of the cases. Still, there remains a serious chance of favoring a subprime density functional for usual data-set sizes ($N<40$)[34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 56, 57, 58, 59] employed in calibration studies of the isomer shift. RMSE and RMPV produced equivalent density functional rankings. However, due to the low experimental resolution, important effects stemming from parameter uncertainty may be hidden. We examined this hypothesis by artificially increasing the experimental resolution from 0.01 mm s*-1* to 0.001 mm s*-1*. While RMSE and RMPV still yielded similar rankings, their dispersion increased significantly. For instance, the “best” density functional, PBE0, is placed first in less than 40% of the cases for 5 data points. Even for 39 data points, the identification of PBE0 as “best” density functional is only successful in 70–75% of the cases.

Finally, we discussed our observation that hybrid density functionals yielded significantly lower MPU estimates compared to pure density functionals. When varying the exact exchange admixture of the B3LYP density functional (for which the linear exact exchange parameter $c_{3}$[105] was originally set to $c_{3}=0.20$), we found a linear dependence of the ACED on $c_{3}$ for all 44 iron compounds studied ($r^{2}>0.9995$). However, the first derivative of the ACED with respect to $c_{3}$ is complex-specific, which is why the MPU is a nonlinear function of $c_{3}$. It is minimized for $c_{3}=0.20$–$0.25$, matching with our observation that B3LYP, B3PW91 ($c_{3}=0.20$), and PBE0 ($c_{3}=0.25$) yielded isomer-shift models with the lowest prediction uncertainty.

When selecting a density functional for actual applications, one is interested in a reliable uncertainty estimation of an isomer-shift prediction for a given ACED. For this purpose, the prediction bands employed in this study represent locally resolved MPU based on the .632 estimator. The validity of the MPU is highest at the mean ACED for a given density functional and decreases continuously from there, which is why extrapolations outside the input domain studied are not recommended. This limitation motivated us to cover a wide range of possible ACED values such that chemically diverse iron complexes can be investigated on the basis of our calibration analysis.

Another possibility to approach accurate isomer-shift predictions is to calculate the isomer shift directly from first principles. Filatov derived the corresponding computational scheme[115] and applied it to examine the isomer shift of molecular iron compounds [116]. Post-calibration of the resulting isomer-shift pairs (measured versus calculated) would allow for an assessment of the more adequate input variable (isomer shift versus ACED).

We discussed several issues in this paper ranging from error diagnostics, MPU estimation, and the role of experimental resolution to the reliability of conclusions drawn on the basis of a specific data set. Bootstrapping clearly shows that performance assessment of density functionals is error-prone for a small number of data points, and for univariate linear regression models, one should consider at least about 40 data points. While the main message of previous calibration studies of the isomer shift remains (hybrid density functionals perform superior to pure density functionals[35, 36, 42, 49, 50, 55]), we now provide a solid statistical framework to examine the certainty of such findings. This framework allowed us to suggest a new, statistically well-justified property model for the isomer shift, given in Eqs. (3.5) and (28), based on PBE0 calculations. Moreover, we introduced a new reference data set for such parametrizations, MIS39, which may be employed and extended for future work. In particular, this calibration study presents the first statistically rigorous analysis for theoretical Mössbauer spectroscopy providing the practitioner with reliable, locally resolved uncertainties for isomer-shift predictions. Moreover, our calibration analysis is of general applicability and not restricted to property models applied in Mössbauer spectroscopy. We will provide the scripts which were written for this work on our webpage so that the methodology can be applied in other parametrization studies of physical and chemical property models.

Acknowledgments

Financial support from the Swiss National Science Foundation (project no. 200020_169120) is gratefully acknowledged.

Appendix

Statistical calibration analysis

All statistical results presented in this work were produced with our suite of scripts reBoot[62] developed in the GNU Octave[117] programming language that is mostly compatible with Matlab. The basic functionalities of the calibration methods employed are already described in Section 2. In this subsection, we discuss some further details.

Bayesian linear regression

Bayesian linear regression[19] is an efficient calibration procedure for normally distributed parameters. The best-fit parameter vector, $\mathbf{w}_{\text{MAP}}$, is obtained from

[TABLE]

where $\mathbf{S}$ is the covariance matrix of the parameters,

[TABLE]

$\mathbf{I}$ is the $(M+1)\times(M+1)$ unit matrix ($M$ $+$ $1$ is the number of parameters contained in the property model), and $\alpha$ and $\beta$ are so-called hyperparameters,

[TABLE]

[TABLE]

which need to be iteratively refined, a procedure referred to as evidence approximation or generalized maximum likelihood.[19] For this purpose, the effective number of parameters, $\gamma$, is calculated according to

[TABLE]

where the $\lambda_{m}$ are eigenvalues of the $(\mathbf{X}^{\top}\mathbf{X})$ matrix of dimension $(M+1)\times(M+1)$,

[TABLE]

and the $\mathbf{v}_{m}$ are the corresponding eigenvectors. For the initialization of $\alpha$ and $\beta$, one can simply choose $\gamma=M+1$ and $\mathbf{w}_{\text{MAP}}=\mathbf{w}_{\mathcal{D}}$. Finally, the locally resolved MPU at input value $\mathbf{x}_{0}$, $s(\mathbf{x}_{0})$, is estimated as

[TABLE]

The hyperparameter $\beta$ represents the inverse noise variance, whereas $\alpha$ represents the belief on the parameter distributions prior to considering the actual data. Choosing $\alpha=0$, one assumes prior parameter distributions with infinite variance, which is equivalent to assigning each value on the real line the same probability density. In that case, Eq. (31) reduces to Eq. (11) resulting in $\mathbf{w}_{\text{MAP}}=\mathbf{w}_{\mathcal{D}}$, i.e., one simply performs linear least-squares regression and additionally obtains the covariance matrix of the model parameters. Here, we observe that $\alpha$ optimized with the evidence approximation and $\alpha$ constrained to zero lead to values for the RMPV deviating by less than 0.05% (i.e., $\mathbf{w}_{\mathcal{D},\text{MAP}}\approx\mathbf{w}_{\mathcal{D}}$). Hence, we can directly compare the RMPV with the R632 obtained from bootstrapped linear least-squares regression where $\alpha=0$ (the same holds for comparisons between RMPV and RLOO).

The evidence approximation is appealing as it provides, for a given initial model complexity (here, the polynomial degree, $M$), a maximum-transferable set of model parameters, $\mathbf{w}_{\text{MAP}}$, by optimizing the hyperparameters $\alpha$ and $\beta$. Hence, the evidence approximation is an approach suited for model selection.[19] With this approach, one also obtains the model parameters that minimize the regularized MSE,[19]

[TABLE]

where $\varepsilon=\alpha/\beta$ is a penalty factor. Given the normal-population assumption basically holds, the evidence approximation might yield a good initial value for the penalty factor to be learned in bootstrapped regularized regression.

In all cases examined in this study, the RMPV (in Eq. (19), the RMPV is defined as its squared variant, MPV) did not change by more than 0.002% after the first iteration of the evidence approximation. This deviation is completely masked by the low experimental resolution of 0.01 mm s*-1* for measurements of the 57Fe Mössbauer isomer shift. Therefore, we conclude that a single iteration is sufficient to obtain a converged RMPV.

Iteratively reweighted linear least-squares regression

In the variant of iteratively reweighted linear least-squares regression[32] by Pernot et al. [16], one starts with a weighted linear least-squares regression with respect to a reference data set, $\mathcal{D}$. The elements of the weight matrix $\mathbf{U}$ introduced in Eq. (13) are updated,

[TABLE]

where the term $d^{2}$ accounts for the discrepancy between the reduced MSE and the average experimental variance, $\langle u\rangle^{2}$,

[TABLE]

Subsequently, one performs another weighted linear least-squares regression with the updated weight matrix. This procedure is repeated until the change in $d^{2}$ becomes negligible between two iterations. Here, we chose a threshold of 0.001%.

Note that the approach outlined here is only of limited applicability as the constant discrepancy factor comes at the expense of two critical assumptions, i.e., a normal-population distribution and homogeneous residual variance (homoscedasticity). For our case study however, both assumptions can be well justified after removal of the inconsistent data points (#1, #2, #7, #13, #28; cf. Table 1 and Fig. 1). On the one hand, the bootstrapped parameter distributions (Fig. 2, bottom) resemble normal distribution, which is indicative of a normal-population distribution. On the other hand, we do not observe a trend of residual variance with respect to the input variable. In general, however, one should critically assess these assumptions by bootstrapping prediction intervals.[27]

Consideration of experimental uncertainty in bootstrapping

Bootstrap samples were drawn by pair resampling (cf. Section 2.1). Each target value, $\delta_{\text{exp},n}$, in a bootstrap sample was randomly perturbed according to its experimental uncertainty (assumed to be 0.02 mm s*-1* for compounds without reported experimental uncertainty). For instance, if a data point considered refers to the $n$-th input–target pair in the reference data set, a random value was drawn from a normal distribution with zero-mean and a standard deviation of $u_{n}$. This random value was rounded to the second decimal place to account for the experimental resolution of 0.01 mm s*-1*, and subsequently added to $\delta_{\text{exp},n}$.

Statistically valid increase of experimental resolution

Bootstrap samples were drawn by pair resampling (cf. Section 2.1). To increase the experimental resolution from its actual value of 0.01 mm s*-1* to 0.001 mm s*-1*, one must respect that any value of the isomer shift with three decimal places is equiprobable as long as the result rounded to two decimal places equals the reported isomer shift. For this purpose, random values were drawn from a uniform distribution with boundaries of $-$0.005 mm s*-1* and $+$0.004 mm s*-1*, and subsequently added to the actual isomer shifts, $\delta_{\text{exp},n}$, assuming the third decimal place equals zero prior to addition. If experimental uncertainty, $u_{n}$, was explicitly considered in bootstrapping (see above), the same procedure was repeated for $u_{n}$.

Quantum chemical calculations

References to the original crystal structures are listed in Table 1. First of all, solute molecules and counterions were removed from all molecular structures. The only exceptions are structures 22 and 26 (Table 1), where we kept the very small counterions due to their proximity to the iron nucleus (lithium and sodium, respectively). All molecular structures were fully optimized with the TPSS density functional [118], Ahlrichs’ def2-TZVP basis set on all atoms [119, 120], Grimme’s DFT-D3 dispersion correction [6], and the conductor-like screening model (COSMO)[121] for electrostatic screening with a dielectric constant of $\varepsilon=78$ (water). For iodine, the only group-5 element in the reference set of iron compounds, the effective core potential def2-ecp was employed.[122] Structure optimizations were performed with Turbomole 6.4.0[123] applying the resolution-of-the-identity approximation, which invokes an auxiliary basis set. The convergence thresholds were set to $10^{-7}$ hartree for the electronic energy difference and $10^{-4}$ hartree/bohr for the length of the gradient of the electronic energy with respect to the nuclear coordinates. Some molecular structures were further truncated prior to structure optimization (see Table 1). For instance, larger alkyl groups at aromatic rings not directly attached to iron were replaced by methyl groups.

All subsequent ACED calculations were performed with Orca 3.0.3[124] (see Table S2). For every molecular iron compound contained in our reference set, we determined the ACED with 12 different density functionals; one LDA functional (PWLDA[125]), four GGA functionals (BP86[126, 127], BLYP[127, 128, 129], PW91[130], PBE[131]), two meta-GGA functionals (M06-L[132], TPSS), three hybrid-GGA functionals (B3LYP[128, 129, 133], B3PW91[133, 131], PBE0[134]), and two meta-hybrid-GGA functionals (M06[135], TPSSh[136]). Neese’s CP(PPP) basis set [35] was employed for iron, whereas Ahlrichs’ def2-TZVP basis set was employed for all other elements. All calculations were based on electrostatic screening (COSMO) with a dielectric constant of $\varepsilon=80$. The convergence threshold for the electronic energy difference was set to $10^{-6}$ hartree (default). In case of error code 16384 (the only type of error that occurred in our ACED calculations), we added either the keyword slowconv to the first line (starting with !) or %scf maxiter <no.> end in a separate line (<no.> $=$ 500 always worked in the cases studied).

We repeated all ACED calculations with def2-TZVP replaced by def2-SVP[119] (see Table S3), resulting in slightly higher MPU for the latter (cf. Tables S6 and S7). We did not perform further statistical analyses based on the def2-SVP results, but exploited this smaller basis set to probe the sensitivity of the ACED with respect to the convergence threshold for the electronic energy difference (see Table S4). We decreased the convergence threshold for the electronic energy difference to $10^{-8}$ hartree (B3LYP/def2-SVP). For both thresholds, mean and maximum of the unsigned difference in the ACED are only 0.1% and 3.7% as large as the standard deviation of the ACED itself, respectively. We may conclude that default self-consistent field convergence criteria lead to stable results, especially taking into account the low experimental resolution of 0.01 mm s*-1*.

Supporting Information

All Cartesian coordinates, contact density data, and further statistical analysis can be found in the supporting information. This information is available free of charge via the Internet at http://pubs.acs.org/.