Condition number calculation

Diagnostics
1

Hello,

I am evaluating a PK and PK-PD model with rich clinical trials data. Overall, the condition numbers are not expected to be too high as it is relatively simple model. The model has been fitted previously with NONMEM and condition numbers were consistently < 1000. I am trying to get condition numbers from Pumas model fit with infer() and then cond(). I am getting very high condition numbers > 10^6. I have tried several variations of the model with simplification but condition number remain high.

Can you please discuss how condition number calculations differ between NONMEM COV step and what is being used by Pumas default (sandwich matrix)? How should the comparison be made? I tried using Pumas.SIR but cond() function don’t work on it. Bottomline is that I am trying to understand the differences and see how can we correct or justify large condition number from Pumas model run?

Thanks!

2

Krina - welcome to the Pumas community.

Can you please tell us if we matched the results between nomem and Pumas?

1 Like
3

The condition number we report is the ratio of the largest and smallest eigenvalues of the covariance matrix estimate. The sandwich estimator is used by default when you call infer. If you want to use the inverse of the observed Fischer information matrix instead of the sandwich estimator, you can set sandwich_estimator = false in the call to infer, e.g. infer(fpm, sandwich_estimator = false).

The condition number of a covariance matrix is a measure of how close that matrix is to being singular. There are 2 main reasons why a covariance matrix may be close to being singular:

  1. Some parameters have a large standard error compared to some other parameters
  2. Some parameters are highly correlated due to poor identifiability of the model’s parameters.

An example of the first case is:

julia> C = [1.0 0.0;
            0.0 1e-9]
2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0e-9
julia> cond(C)
9.999999999999999e8

An example of the second case is:

julia> C = [1.0   0.9999;
           0.9999    1.0]
2×2 Matrix{Float64}:
 1.0     0.9999
 0.9999  1.0
julia> cond(C)
19998.999999980093

In your example, it can be a mixture of the 2 factors. One less common third scenario that may also lead to spuriously high condition numbers is if the fitting failed to find a local optimal solution. Instead it could have converged to a saddle point or to a point where the log likelihood is not a smooth function so the Hessian is not well-defined or has no statistical meaning. This can lead to negative or high positive condition numbers but we can’t draw statistical conclusions from that.

When diagnosing a model, it is perhaps more useful to inspect the standard errors and correlation matrix directly instead of looking at the condition number which is a summary figure that can be high due to the multiple reasons mentioned above.

A high condition number (beside the statistical interpretation) also has numerical consequences. When doing various computations involving matrices, the higher the condition number of the matrix, the larger the computational relative error that will accumulate due to round-off errors. This is a purely numerical consideration but it may lead in theory to spurious failed fits. During the fitting process, the optimization algorithm uses an approximation of the Hessian to guide the selection of a search direction in every iteration. If the error is too high, bad search directions may be chosen, leading to slow or failing fits. If this happens, this may be an instance of The folk theorem of statistical computing which says: “When you have computational problems, often there’s a problem with your model.“.

2 Likes
4

Thanks, Vijay! Yes, the models have been compared between NONMEM and Pumas.

Thanks, Mohamed for the explanation! I will try without sandwich matrix.

5

Hi @vijay and @mohamed82008 ,

I have tried sandwich_estimator = false. I have carefully looked into the SEs and all SEs look very reasonable. The parameter estimates are very very close to what NONMEM had previously provided - for the same dataset and comparable model. From what it seems all other checks look good for the models. But, all the models that I have tried for this project show high condition numbers.

Is there an easy way to see correlation matrix? From infer() function, I see coeftable but can’t examine the correlation matrix.

Any other thoughts regarding condition numbers? We are working on this project with Pumas consulting team, so I have shared the latest models with Joga. I will appreciate it if you can please look into the models for condition number check.

Thanks,
Krina

6

hi @krina.mehta.bb - you can look at the correlation matrix for any inferred object using the cor2cov function. I just realized that it had not been exported so you can access it as shown below. Regarding condition numbers, I am on a different side of the fence where I would not worry too much about it, but we will expand on some more details shortly.

help?> PumasReports.cor2cov
  cor2cov(C, s)

  Compute the covariance matrix from the correlation matrix C and a vector of standard deviations s. Use StatsBase.cor2cov! for an in-place version.

julia> PumasReports.cov2cor(fit_infer.vcov)
5×5 Matrix{Float64}:
  1.0         0.0370938   0.0149482  -0.0967604  -0.120237
  0.0370938   1.0         0.105866    0.202644   -0.418523
  0.0149482   0.105866    1.0        -0.0710371  -0.14918
 -0.0967604   0.202644   -0.0710371   1.0         0.0115326
 -0.120237   -0.418523   -0.14918     0.0115326   1.0
7

Thanks Vijay! At this point, since all other checks look good, I am not worried about condition number for this project. However, in general, we use condition number as a quick check to make sure that the model is not over-parameterized. Therefore, additional checks/information regarding condition number calculations between NONMEM and Pumas will be very helpful in the future.

8

For this I would recommend looking at the condition number of the correlation matrix instead of the covariance matrix.

9

You can use the following:

cond(Pumas.cov2cor(inf.vcov.Σ))

where inf is the output from infer.