Reliability of a fitted curve?

Question

I would like to estimate the uncertainty or the reliability of a fitted curve. I intentionally don't name a precise mathematical quantity that I am looking for, since I don't know what it is.

Here $E$ (energy) is the dependent variable (response) and $V$ (volume) is the independent variable. I would like to find the Energy-Volume curve, $E(V)$, of some material. So I made some calculations with a quantum chemistry computer program to get the energy for some sample volumes (green circles in the plot).

Then I fitted these data samples with the Birch–Murnaghan function: $$ \mathbb{E}(E|V) = E_0 + \frac{9V_0B_0}{16} \left\{ \left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^3B_0^\prime + \left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^2 \left[6-4\left(\frac{V_0}{V}\right)^\frac{2}{3}\right]\right\}\;, $$ which depends on four parameters: $E_0, V_0, B_0, B_0'$. I also assume that this is the correct fitting function, so all errors just come from the noise of the samples. In what follows, the fitted function $(\hat{E})$ will be written as a function of $V$.

Here you can see the result (fitting with a least squares algorithm). The y-axis variable is $E$ and the x-axis variable is $V$. The blue line is the fit and the green circles are the sample points.

I now need some measure of the reliability (at best in dependence of the volume) of this fitted curve, $\hat{E}(V)$, because I need it to calculate further quantities like transition pressures or enthalpies.

My intutition tells me that the fitted curve is most reliable in the middle, so I guess that the uncertainty (say uncertainty range) should increase near the end of the sample data, like in this sketch:

However, what it this kind of measure that I am looking for and how can I calculate it?

To be precise, there is actually only one error source here: The calculated samples are noisy due to computational limits. So if I would calculate a dense set of data samples they would form a bumpy curve.

My idea to find the desired uncertainty estimate is to calculate the following ''error'' based on the parameters as you learn it in school (propagation of uncertainty):

$$ \Delta E(V) = \sqrt{ \left(\frac{\partial E(V)}{\partial E_0} \Delta E_0\right)^2 + \left(\frac{\partial E(V)}{\partial V_0} \Delta V_0\right)^2 + \left(\frac{\partial E(V)}{\partial B_0} \Delta B_0\right)^2 + \left(\frac{\partial E(V)}{\partial B_0'} \Delta B_0'\right)^2} $$ The $\Delta E_0, \Delta V_0, \Delta B_0$ and $\Delta B_0'$, are given by the fitting software.

Is that an acceptable approach or am I doing it wrong?

PS: I know that I could also just sum up the squares of the residuals between my data samples and the curve to get some kind of ''standard error'' but this is not volume dependent.

Answer 1

This is an ordinary least squares problem!

Defining

$$x = V^{-2/3}, \ w = V_0^{1/3},$$

the model can be rewritten

$$\mathbb{E}(E|V) = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3$$

where the coefficients $\beta=(\beta_i)^\prime$ are algebraically related to the original coefficients via

$$16 \beta = \pmatrix{16 E_0 + 54 B_0 w^3 - 9 B_0 B_0^\prime w^3\\ -144 B_0 w^5 + 27 B_0 B_0^\prime w^5\\ 126 B_0 w^7 - 27 B_0 B_0^\prime w^7\\ -36 B_0 w^9 + 9 B_0 B_0^\prime w^9}.$$

This is straightforward to solve algebraically or numerically: pick the solution for which $B_0, B_0^\prime$, and $w$ are positive. The only reason to do this is to obtain estimates of $B_0, B_0^\prime, w$, and $E_0$ and to verify they are physically meaningful. All analyses of the fit can be carried out in terms of $\beta$.

This approach is not only much simpler than nonlinear fitting, it is also more accurate: the variance-covariance matrix for $(E_0, B_0, B_0^\prime, V_0)$ returned by a nonlinear fit is only a local quadratic approximation to the sampling distribution of these parameters, whereas (for normally distributed errors in measuring $E$, anyway) the OLS results are not approximations.

Confidence intervals, prediction intervals, etc. can be obtained in the usual way without needing to find these values: compute them in terms of the estimates $\hat\beta$ and their variance-covariance matrix. (Even Excel could do this!) Here is an example, followed by the (simple) R code that produced it.

#
# The data.
#
X <- data.frame(V=c(41, 43, 46, 48, 51, 53, 55.5, 58, 60, 62.5),
                E=c(-48.05, -48.5, -48.8, -49.03, -49.2, -49.3, -49.35, 
                    -49.34, -49.31, -49.27))
#
# OLS regression.
#
fit <- lm(E ~ I(V^(-2/3)) + I(V^(-4/3)) + I(V^(-6/3)), data=X)
summary(fit)
beta <- coef(fit)
#
# Prediction, including standard errors of prediction.
#
V0 <- seq(40, 65)
y <- predict(fit, se.fit=TRUE, newdata=data.frame(V=V0))
#
# Plot the data, the fit, and a three-SEP band.
#
plot(X$V, X$E, xlab="Volume", ylab="Energy", bty="n", xlim=c(40, 60))
polygon(c(V0, rev(V0)), c(y$fit + 3*y$se.fit, rev(y$fit - 3*y$se.fit)),
        border=NA, col="#f0f0f0")
curve(outer(x^(-2/3), 0:3, `^`) %*% beta, add=TRUE, col="Red", lwd=2)
points(X$V, X$E)

---

If you are interested in the joint distribution of the original parameter estimates, then it is easy to simulate from the OLS solution: simply generate multivariate Normal realizations of $\beta$ and convert those into the parameters. Here is a scatterplot matrix of 2000 such realizations. The strong curvilinearity shows why the Delta method is likely to give poor results.

Answer 2

There is a standard approach for this called the delta method. You form the inverse of the Hessian of the log-likelihood wrt your four parameters. There is an extra parameter for the variance of the residuals, but it does not play a role in these calculations. Then you calculate predicted response for the desired values of the independent variable and calculate its gradient (the derivative wrt) these four parameters. Call the inverse of the Hessian $I$ and the gradient vector $g$. You form the vector matrix product $$-g^tIg$$ This gives you the estimated variance for that dependent variable. Take the square root to get the estimated standard deviation. then the confidence limits are the predicted value +- two standard deviations. This is standard likelihood stuff. for the special case of a nonlinear regression you can correct for the degrees of freedom. You have 10 observations and 4 parameters so you can increase the estimation of the variance in the model by multiplying by 10/6. Several software packages will do this for you. I wrote up your model in AD Model in AD Model Builder and fit it and calculated the (unmodified) variances. They will be slightly different from yours because I had to guess a bit at the values.

                    estimate   std dev
10   pred_E      -4.8495e+01 7.5100e-03
11   pred_E      -4.8810e+01 7.9983e-03
12   pred_E      -4.9028e+01 7.5675e-03
13   pred_E      -4.9224e+01 6.4801e-03
14   pred_E      -4.9303e+01 6.8034e-03
15   pred_E      -4.9328e+01 7.1726e-03
16   pred_E      -4.9329e+01 7.0249e-03
17   pred_E      -4.9297e+01 7.1977e-03
18   pred_E      -4.9252e+01 1.1615e-02

This can be done for any dependent variable in AD Model Builder. One declares a variable in the appropriate spot in the code like this

   sdreport_number dep

and writes the code the evaluate the dependent variable like this

dep=sqrt(V0-cube(Bp0)/(1+2*max(V)));

Note this is evaluated for a value of the independent variable 2 times the largest one observed in the model fitting. Fit the model and one obtains the standard deviation for this dependent variable

19   dep          7.2535e+00 1.0980e-01

I have modified the program to include code for calculating the confidence limits for the enthalpy-volume function The code (TPL) file looks like

DATA_SECTION
 init_int nobs
 init_matrix data(1,nobs,1,2)
 vector E
 vector V
 number Vmean
LOC_CALCS
 E=column(data,2);
 V=column(data,1);
 Vmean=mean(V);

PARAMETER_SECTION
 init_number E0
 init_number log_V0_coff(2)
 init_number log_B0(3)
 init_number log_Bp0(3)
 init_bounded_number a(.9,1.1)
 sdreport_number V0
 sdreport_number B0
 sdreport_number Bp0
 sdreport_vector pred_E(1,nobs)
 sdreport_vector P(1,nobs)
 sdreport_vector H(1,nobs)
 sdreport_number dep
 objective_function_value f
PROCEDURE_SECTION
  V0=exp(log_V0_coff)*Vmean;
  B0=exp(log_B0);
  Bp0=exp(log_Bp0);
  if (current_phase()<4)
  f+=square(log_V0_coff) +square(log_B0);

  dvar_vector sv=pow(V0/V,0.66666667);
  pred_E=E0 + 9*V0*B0*(cube(sv-1.0)*Bp0
    + elem_prod(square(sv-1.0),(6-4*sv)));

  dvar_vector r2=square(E-pred_E);
  dvariable vhat=sum(r2)/nobs;
  dvariable v=a*vhat;
  f=0.5*nobs*log(v)+sum(r2)/(2.0*v);

  // code to calculate the  enthalpy-volume function
  double delta=1.e-4;
  dvar_vector svp=pow(V0/(V+delta),0.66666667);
  dvar_vector svm=pow(V0/(V-delta),0.66666667);
  P = -((9*V0*B0*(cube(svp-1.0)*Bp0
      + elem_prod(square(svp-1.0),(6-4*svp))))
      -(9*V0*B0*(cube(svm-1.0)*Bp0
      + elem_prod(square(svm-1.0),(6-4*svm)))))/(2.0*delta);
  H=E+elem_prod(P,V);

dep=sqrt(V0-cube(Bp0)/(1+2*max(V)));

Then I refitted the model to get the standard devs for the estimates of H.

29   H           -3.9550e+01 5.9163e-01
30   H           -4.1554e+01 2.8707e-01
31   H           -4.3844e+01 1.2333e-01
32   H           -4.5212e+01 1.5011e-01
33   H           -4.6859e+01 1.5434e-01
34   H           -4.7813e+01 1.2679e-01
35   H           -4.8808e+01 1.1036e-01
36   H           -4.9626e+01 1.8374e-01
37   H           -5.0186e+01 2.8421e-01
38   H           -5.0806e+01 4.3179e-01

These are calculated for your observed V values, but could be easily calculated for any value of V.

It has been pointed out that this is actually a linear model for which there is simple R code to perform the parameter estimation via OLS. This is very appealing especially to naive users. However since the work of Huber over thirty years ago we know or should know that one should probably almost always replace OLS with a moderately robust alternative. The reason this is not routinely done I believe is that robust methods are inherently nonlinear. From this point of view the simple appealing OLS methods in R are more of a trap, rather than a feature. An advabtage of the AD Model Builder approach is its built in support for nonlinear modelling. To change the least squares code to a robust normal mixture only one line of the code needs to be changed. The line

    f=0.5*nobs*log(v)+sum(r2)/(2.0*v);

is changed to

f=0.5*nobs*log(v)
  -sum(log(0.95*exp(-0.5*r2/v) + 0.05/3.0*exp(-0.5*r2/(9.0*v))));

The amount of overdispersion in the models is measured by the parameter a. If a equals 1.0, the variance is the same as for the normal model. If there is inflation of the variance by outliers we expect that a will be smaller than 1.0. For these data the estimate of a is about 0.23 so that the variance is about 1/4 the variance for the normal model. The interpretation is that outliers have increased the variance estimate by a factor of about 4. The effect of this is to increase the size of the confidence bounds for parameters for the OLS model. This represents a loss of efficiency. For the normal mixture model the estimated standard deviations for the enthalpy-volume function are

 29   H           -3.9777e+01 3.3845e-01
 30   H           -4.1566e+01 1.6179e-01
 31   H           -4.3688e+01 7.6799e-02
 32   H           -4.5018e+01 9.4855e-02
 33   H           -4.6684e+01 9.5829e-02
 34   H           -4.7688e+01 7.7409e-02
 35   H           -4.8772e+01 6.2781e-02
 36   H           -4.9702e+01 1.0411e-01
 37   H           -5.0362e+01 1.6380e-01
 38   H           -5.1114e+01 2.5164e-01

One sees that there are small changes in the point estimates, while the confidence limits have been reduced to about 60% of the ones produced by OLS.

The main point I want to make is that all the modified calculations occur automatically once one changes the one line of code in the TPL file.

Answer 3

Cross-validation is a simple way to estimate reliability of your curve: https://en.wikipedia.org/wiki/Cross-validation_(statistics)

Propagation of uncertainty with partial differentials is great is you really know $\Delta E_{0},\Delta V_{0},\Delta B_{0}$, and $\Delta B'$. However, the program you are using gives only fitting errors(?). Those will be too optimistic (unrealistically small).

You can calculate 1-fold validation error by leaving one of your points away from fitting and using the fitted curve to predict value of the point that was left away. Repeat this for all points so that each is left away once. Then, calculate validation error of your final curve (curve fitted with all points) as an average of prediction errors.

This will only tell you how sensitive your model is for any new data point. For example, it will not tell you how inaccurate your energy model is. However, this will be much more realistic error estimate mere fitting error.

Also, you can plot prediction errors as a function of volume if you want.