JSim calculates the correlation between parameters two different ways: Normalizing the inverse of an estimation of the Hessian and a Monte Carlo approach. The results are not necessarily similar.

Description

The meaning of “correlation of model parameters”

INSTRUCTIONS:

     Run the model. Three curves will appear when you click on the Concentrations TAB
(upper right corner). The three curves are (1) Cout, the model outflow curve (black 
stars), Ctrunc (values of Cout <1e-6 set to zero (red circles), and Cnoise (blue line)
which equals Ctrunc multiplied by (1+ Gaussian noise with mean 0 and standard 
deviation 0.1),random seed = 12345.
     Switch to the optimizer GUI (TAB at lower left hand of Run Time GUI) and run the
optimizer which will vary four parameters to fit Cnoise with Cout. Click the Report TAB
(located near upper Optimizer GUI page) and either print or export the text. You will
need the correlation matrix from the optimizer to compare with the correlation matrix 
from the Monte Carlo simulation.
     Now click on the Monte Carlo TAB (bottom of Optimizer GUI) and run the Monte Carlo
run. Click on the Report TAB near the top of the Monte Carlo GUI. Find the correlation
matrix and compare the values to those returned by the optimizer.

Abstract
	A commonly used term, the “correlation of model parameters” has been used 
to describe two different mathematical procedures in JSim  which  give different results.

Introduction
	There are two  distinct meanings of the phrase “correlation of model
 parameters” in JSim. The first meaning is based on  the sensitivity functions 
for a given set of parameters. The second is the correlation of multiple 
estimates of those parameters using a Monte Carlo approach adding noise to optimized
curves. A simple model will be introduced to show the difference. The analysis was done with 
JSim[1].

        The following text is adapted from reference 3. Some parts are verbatim without
        using quotes.
     -------------------------------------------------------------------------------------
        Notation: Let the observed data be given as

        {y  = f(P,x ) + r      i=1..m }
          i        i     i

        where f(P,x ) is a mathematical model, a function of the parameter vector P
                   i
        of length n, x  is the independent variable defined at m points, and the r 
                      i                                2                          i
        are random errors with mean zero and variance s  . In vector notation,
                                                       i
        Y = F(P,X) + R.

        A weighting matrix, W, of size m x m, W is defined as 
                                             2
        W = diag(w , w , ... w ), where w = s .
                  1   2       m          i   i

        The m x 1 residual vector R is given by

        R(P) = Y-F(P,X).

        The sum of the squares of the weighted residuals, SSR(P) is denoted by
                     T
        SSR(P) = R(P)  W  R(P).


        The sensitivity matrix S(p ,x ,h  ) = (f(P+h  *I ,x ) - f(P,x ))/h
                                  j  i  ij          ij  j  i         i    ij
        for 1<=i<=m and 1<=j<=n, is an approximation of the Jacobian matrix where I is the
                                                                                   j
        j'th column of a n x n identity matrix and h   is the step size or mesh spacing
                                                    ij
        for j'th parameter and i'th observation.
                                             
        The covariance matrix of the solution parameters can be estimated by the Hessian

        matrix at the solution (i.e., the second derivative matrix).                                                  
        

	Near the solution for small residuals the co-variance matrix of P from an 
  
        optimization run is estimated by
                               T        -1              
        Cov(P) = SSR/(nh-nx)*[S * W * S]   (see references 3 and 5).

        The n x n correlation matrix is given by

        Correlation(P) = Cov(k,l)/sqrt(Cov(k,k)*Cov(l,l)). 

    ---------------------------------------------------------------------------------------------

        The definition of correlation of multiple estimates of parameters from the

        Monte Carlo simulation is

                              N
                            -----         _          _
                             \    (A(i) - A) (B(i) - B)
                              )   ---------------------------
                             /               N - 1
                            -----
                            i = 1
                   corr   = --------------------------------- .
                       AB               sigma  sigma
                                             A      B

        It is a function of the optimized parameter values and is unrelated to the sum 

        of the squares of the residuals, the number of degrees of freedom or the sensitivity 

        matrix.

Methods and Results

Consider a capillary modeled by an advection-diffusion partial differential equation
[2,4], governed by

C:t = (-F * L/V) * C:x -(G/V)  * C + D * C:x:x;

with inital and bundary conditions given by

when(t=t.min) C=0 mM;
when(x=x.min) (-F*L/V)*(C-Cin) +D*C:x =0;
when(x=x.max) {C:x=0; Cout=C; }
The model was run and the outflow concentration values were truncated to zero when 
they were  less than 1e-6 mM. Ten percent  proportional noise was added to the
truncated output., i.e.,

Cnoise = Ctrunc * ( 1 + 0.1 * randomg() );

and saved as Cnoise. The function, randomg(), generates Gaussianly distributed 
random numbers with mean equal to zero and variance equal to one. The initial 
seed was set to 12345.
 
Correlation of parameters from the estimated Hessian

JSim was set to optimize parameters F, V, G, and log10Dp, 
fitting Cout to Cnoise. The point weights for Cnoise were set to 1.0. The Sensop[3] 
optimizer was used with the maximum number of model runs set to 100. LSFEA3[4]  
was used as the partial differential equation solver. 

The correlation matrix from optimization based on the estimation of 
the Hessian yielded the following correlations between the parameters.
    
 Correlation matrix: (condition number=3.1872574E5)
                     F             V             G          log10Dp   
       F             1         .92197543    -0.07320687   -0.8004058  
       V         .92197543         1        -0.15421016   -0.82221498 
       G        -0.07320687   -0.15421016        1        -0.20963089 
    log10Dp     -0.8004058    -0.82221498   -0.20963089        1


Correlation of parameters using the Monte-Carlo method (Pearson product-moment
correlation coefficient)

The model was run in JSim's  Monte-Carlo mode for 100 optimizations, i. e. 10,000
 model runs, with 10% proportional Gaussian noise added to the truncated outflow 
concentration data. This requires changing the optimized data being fit to Ctrunc.
 
Correlation Matrix:
                  F           V           G        log10Dp  
      F           1         .8154       .0979      -0.2515  
      V         .8154         1         .2376      -0.1917  
      G         .0979       .2376         1         .0491   
   log10Dp     -0.2515     -0.1917      .0491         1      

Notice that the correlations of F with G, F with log10Dp, V with G, and G with log10Dp
change signs. 

Conclusion
If we consider the magnitudes of the values reported, clearly F with log10Dp and 
V with log10Dp are very different. The differences between the Monte-Carlo correlations
and the correlations from an individual optimization based on inverting the Hessian
indicates that they are unrelated. It is not surprising that the values are different 
because they are calculated by completely different processes. We have been misled by 
the use of the same descriptive term, “correlation of parameters”, used to describe 
very different processes.

Equations

None.

The equations for this model may be viewed by running the JSim model applet and clicking on the Source tab at the bottom left of JSim's Run Time graphical user interface. The equations are written in JSim's Mathematical Modeling Language (MML). See the Introduction to MML and the MML Reference Manual. Additional documentation for MML can be found by using the search option at the Physiome home page.