Identifiability and Estimability Analysis#
This notebook demonstrates the three identifiability diagnostics added in response to issue #45:
diagnose_identifiability— Belsley/Kuh/Welsch regression diagnostics plus the Gutenkunst sloppy-model eigenvalue spectrum [Belsley et al., 1980, Gutenkunst et al., 2007].estimability_rank— Yao-style orthogonalization ranking with the Brun collinearity index on top of discopt’s existingcompute_fimJacobian [Brun et al., 2001, Yao et al., 2003].profile_likelihood— Raue/Kreutz profile-likelihood confidence intervals and shape classification [Kreutz et al., 2013, Raue et al., 2009].
The worked example is a four-parameter series reaction kinetics model where two rate constants enter only as their product, producing a classic practical non-identifiability. We use this to demonstrate each diagnostic and show that they agree on which direction is non-identifiable [Miao et al., 2011, Villaverde, 2019].
Reparameterization note. The Yao ranking and the Belsley condition indices are not invariant under a change of variable such as \(\theta \to \log \theta\). Profile likelihood is invariant. If a diagnostic flags non-identifiability, try a log-reparameterization before concluding the parameter itself is at fault.
import os
os.environ.setdefault("JAX_PLATFORMS", "cpu")
os.environ.setdefault("JAX_ENABLE_X64", "1")
import discopt.modeling as dm
import matplotlib.pyplot as plt
import numpy as np
from discopt.doe import (
diagnose_identifiability,
estimability_rank,
profile_all,
profile_likelihood,
)
from discopt.estimate import Experiment, ExperimentModel, estimate_parameters
np.random.seed(42)
The model#
Series first-order reactions \(A \xrightarrow{k_1} B \xrightarrow{k_2} C\) observed at a set of time points, with an additional constant baseline and a quadratic nuisance term.
We build the response as
so that:
The product \(k_1 k_2\) (and \(\alpha\)) appears linearly through the single combined parameter \(\alpha k_1 k_2\); only two of \(\{\alpha, k_1, k_2\}\) are algebraically identifiable.
\(c_0\) and \(c_2\) are fully identifiable on their own.
This is the canonical “practical non-identifiability by overparameterization” scenario [Yao et al., 2003].
class KineticsExperiment(Experiment):
def __init__(self, t_data):
self.t_data = t_data
def create_model(self, **kwargs):
m = dm.Model("kinetics")
alpha = m.continuous("alpha", lb=0.1, ub=5.0)
k1 = m.continuous("k1", lb=0.01, ub=5.0)
k2 = m.continuous("k2", lb=0.01, ub=5.0)
c0 = m.continuous("c0", lb=-5.0, ub=5.0)
c2 = m.continuous("c2", lb=-1.0, ub=1.0)
responses = {}
errors = {}
for i, t in enumerate(self.t_data):
responses[f"y_{i}"] = alpha * k1 * k2 * t + c0 + c2 * t * t
errors[f"y_{i}"] = 0.1
return ExperimentModel(
m,
unknown_parameters={"alpha": alpha, "k1": k1, "k2": k2, "c0": c0, "c2": c2},
design_inputs={},
responses=responses,
measurement_error=errors,
)
# True parameter values (note: alpha*k1*k2 = 1.5)
alpha_true = 1.0
k1_true = 1.5
k2_true = 1.0
c0_true = 0.2
c2_true = 0.05
t_data = np.linspace(0.5, 5.0, 10)
exp = KineticsExperiment(t_data)
# Generate synthetic data with Gaussian noise.
data = {}
for i, t in enumerate(t_data):
mu = alpha_true * k1_true * k2_true * t + c0_true + c2_true * t * t
data[f"y_{i}"] = mu + np.random.normal(0, 0.1)
fig, ax = plt.subplots(figsize=(6, 3.5))
ax.plot(t_data, [data[f"y_{i}"] for i in range(len(t_data))], "o", label="observed")
ax.set_xlabel("t")
ax.set_ylabel("y")
ax.set_title("Synthetic kinetics data")
ax.legend()
<matplotlib.legend.Legend at 0x108b34440>
Step 1 — Fit the full model#
Use estimate_parameters to obtain maximum-likelihood estimates
\(\hat{\theta}\) and the deviance \(D(\hat{\theta})\) at the optimum.
est = estimate_parameters(
exp,
data,
initial_guess={
"alpha": alpha_true,
"k1": k1_true,
"k2": k2_true,
"c0": c0_true,
"c2": c2_true,
},
)
print(est.summary())
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
Parameter Estimation Results
==================================================
alpha = 1.34885 ± 5.203e+05 [-1.337e+06, 1.337e+06]
k1 = 1.06522 ± 4.337e+05 [-1.115e+06, 1.115e+06]
k2 = 1.06522 ± 3.707e+05 [-9.53e+05, 9.53e+05]
c0 = 0.216697 ± 0.1176 [-0.08564, 0.519]
c2 = 0.0441992 ± 0.01741 [-0.0005489, 0.08895]
Objective = 4.5897
N obs = 10
FIM det = 1.571e-13
FIM cond = 2.518e+17
Note that the point estimates of \(\alpha\), \(k_1\), \(k_2\) will not agree with the true values — because only their product is identifiable, the optimizer picks some triple lying on the ridge \(\alpha k_1 k_2 = 1.5\).
Step 2 — Diagnose identifiability#
diagnose_identifiability returns the full Belsley/Gutenkunst
diagnostic bundle. The condition indices, variance-inflation factors,
variance-decomposition proportions, and the log-eigenvalue spectrum all
point at the same offending direction
[Belsley et al., 1980, Gutenkunst et al., 2007].
diag = diagnose_identifiability(exp, est.parameters)
print(f"FIM rank: {diag.fim_rank}/{diag.n_parameters}")
print(f"Condition number: {diag.condition_number:.3g}")
print()
print("Condition indices:", np.round(diag.condition_indices, 3))
print()
print("VIFs:")
for name, v in diag.vif.items():
print(f" {name:>6s}: {v:.3g}")
print()
print("Warnings:")
for w in diag.warnings:
print(" ", w)
FIM rank: 3/5
Condition number: 6.55e+17
Condition indices: [1.00000000e+00 4.34900000e+00 2.04860000e+01 2.22238242e+16
5.38109799e+16]
VIFs:
alpha: inf
k1: inf
k2: inf
c0: 13.8
c2: 48
Warnings:
mild collinearity: condition index eta_3 = 20.5 > 10
serious collinearity: condition index eta_4 = 2.22e+16 > 30
serious collinearity: condition index eta_5 = 5.38e+16 > 30
VIF[alpha] is infinite (parameter lies in a null direction)
VIF[k1] is infinite (parameter lies in a null direction)
VIF[k2] is infinite (parameter lies in a null direction)
VIF[c0] = 13.8 > 10
VIF[c2] = 48 > 10
|rho[alpha,k2]| = 1.83e+07 > 0.95
|rho[k2,c0]| = 1.76 > 0.95
|rho[k2,c2]| = 1.8 > 0.95
# Sloppy-model spectrum (Gutenkunst)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3.5))
ax1.plot(np.arange(1, len(diag.normalized_log_spectrum) + 1), diag.normalized_log_spectrum, "o-")
ax1.set_xlabel("eigenvalue index (descending)")
ax1.set_ylabel(r"$\log_{10}(\lambda / \lambda_\max)$")
ax1.set_title("FIM eigenvalue spectrum")
ax1.axhline(-6, color="red", linestyle="--", alpha=0.3, label="sloppy (-6 decades)")
ax1.legend()
# VIF bar chart
names = list(diag.vif.keys())
vifs = [diag.vif[n] for n in names]
vifs_plot = [min(v, 1e3) for v in vifs]
ax2.bar(names, vifs_plot, color=["red" if np.isinf(v) or v > 10 else "steelblue" for v in vifs])
ax2.axhline(10, color="black", linestyle="--", alpha=0.4, label="VIF=10")
ax2.set_ylabel("VIF (clipped at 1e3)")
ax2.set_title("Variance Inflation Factors")
ax2.set_yscale("log")
ax2.legend()
plt.tight_layout()
The parameter \(k_2\) and \(k_1\) show infinite VIFs (or large finite ones depending on the fit), and the log-eigenvalue spectrum spans multiple decades — the classic sloppy signature.
Inspect the null-space direction directly:
if diag.null_space:
print("Null direction(s) in (normalized) parameter space:")
for d in diag.null_space:
for name, coef in d.items():
print(f" {name:>6s}: {coef:+.4f}")
else:
print("FIM is full rank (no strict null direction at this tolerance).")
Null direction(s) in (normalized) parameter space:
alpha: +0.7871
k1: -0.5815
k2: -0.2056
c0: -0.0000
c2: -0.0000
alpha: -0.2170
k1: -0.5732
k2: +0.7902
c0: +0.0000
c2: +0.0000
Step 3 — Estimability ranking (Yao / Brun)#
estimability_rank ranks parameters from most to least estimable using
rank-revealing QR on the Brun-scaled sensitivity matrix
[Yao et al., 2003]. The projected 2-norm at step \(k\) equals
\(|R_{kk}|\) from the QR and matches Yao’s paper numerically.
rank = estimability_rank(exp, est.parameters)
print("Ranking (most -> least estimable):")
for name, pn in zip(rank.ranking, rank.projected_norms):
bar = "=" * max(1, int(50 * pn / max(rank.projected_norms[0], 1e-30)))
print(f" {name:>6s} |R_kk| = {pn:10.4g} {bar}")
print()
print(f"Recommended subset at cutoff 0.04: {rank.recommended_subset}")
print(f"Collinearity index of the recommended subset: {rank.collinearity_index:.3g}")
Ranking (most -> least estimable):
alpha |R_kk| = 150.2 ==================================================
c2 |R_kk| = 4.372 =
c0 |R_kk| = 1.842 =
k1 |R_kk| = 3.007e-14 =
k2 |R_kk| = 5.206e-15 =
Recommended subset at cutoff 0.04: ['alpha']
Collinearity index of the recommended subset: 1
Step 4 — Profile likelihood#
profile_likelihood steps \(\theta_i\) outward from the estimate in both
directions, re-optimizing all other parameters at each step, and records
the deviance \(D(\theta)\). The 95% confidence region is
\(\{\theta_i : D(\theta_i) - D(\hat\theta) \le \chi^2_{1, 0.95}\}\)
(deviance convention — no factor of 1/2, see module docstring)
[Kreutz et al., 2013, Raue et al., 2009].
profiles = profile_all(exp, data, confidence_level=0.95, max_steps=15)
for name, p in profiles.items():
ci_lo = f"{p.ci_lower:.3g}" if p.ci_lower is not None else "-inf"
ci_hi = f"{p.ci_upper:.3g}" if p.ci_upper is not None else "+inf"
print(f" {name:>6s} shape={p.shape:17s} 95% CI = [{ci_lo}, {ci_hi}]")
for w in p.warnings:
print(f" warn: {w}")
alpha shape=flat 95% CI = [-inf, +inf]
warn: profile is flat: |D(theta) - D(theta_hat)| stays below 0.384; parameter is not practically identifiable
k1 shape=one_sided_lower 95% CI = [0.798, +inf]
warn: upper arm hit parameter bound ub=5 without crossing the threshold
k2 shape=one_sided_lower 95% CI = [0.798, +inf]
warn: upper arm hit parameter bound ub=5 without crossing the threshold
c0 shape=flat 95% CI = [-inf, +inf]
warn: neither arm crossed the threshold; parameter may be practically non-identifiable
c2 shape=flat 95% CI = [-inf, +inf]
warn: neither arm crossed the threshold; parameter may be practically non-identifiable
# Plot profiles for each parameter.
fig, axes = plt.subplots(1, len(profiles), figsize=(3 * len(profiles), 3.5), sharey=True)
for ax, (name, p) in zip(axes, profiles.items()):
ax.plot(p.theta_values, p.neg_log_lik - p.objective_hat, "o-", markersize=3)
ax.axhline(
p.threshold - p.objective_hat, color="red", linestyle="--", label=r"$\chi^2_{1,0.95}$"
)
ax.axvline(p.theta_hat, color="black", linestyle=":", alpha=0.5)
ax.set_title(f"{name} ({p.shape})")
ax.set_xlabel(r"$\theta$")
axes[0].set_ylabel(r"$D(\theta) - D(\hat\theta)$")
axes[0].legend()
plt.tight_layout()
Round-trip consistency#
All three diagnostics should identify the same non-identifiable direction:
diagnose_identifiabilityreports FIM rank < \(n_\text{params}\).estimability_rankputs the redundant parameter last with a projected norm \(\approx 0\).profile_likelihoodshows a flat or one-sided profile for that parameter.
This pipeline-level agreement is what makes the tools trustworthy; if one flags a problem, the others should too [Villaverde, 2019].
non_identifiable = {p.parameter for p in profiles.values() if p.shape != "bounded"}
print("Parameters with non-bounded profile:", non_identifiable)
excluded = set(est.parameter_names) - set(rank.recommended_subset)
print("Parameters excluded by Yao ranking:", excluded)
null_param = set()
for d in diag.null_space:
null_param |= {n for n, c in d.items() if abs(c) > 0.1}
print("Parameters with weight in FIM null space:", null_param)
Parameters with non-bounded profile: {'k2', 'alpha', 'c0', 'k1', 'c2'}
Parameters excluded by Yao ranking: {'k2', 'c0', 'k1', 'c2'}
Parameters with weight in FIM null space: {'k2', 'alpha', 'k1'}
Worked example 2 — Langmuir adsorption isotherm#
The Langmuir isotherm
has two parameters: the saturation capacity \(q_m\) and the affinity \(K\). Their identifiability depends critically on the concentration grid:
In the linear regime (\(Kc \ll 1\)) the model collapses to \(q \approx q_m K c\), so only the product \(q_m K\) is identifiable and the FIM is rank-1.
In the saturated regime (\(Kc \gtrsim 1\)) the curvature exposes \(q_m\) and \(K\) separately and the FIM is full rank.
This is a textbook practical-identifiability case [Brun et al., 2001, Raue et al., 2009].
class LangmuirExperiment(Experiment):
def __init__(self, c_data, sigma=0.05):
self.c_data = list(c_data)
self.sigma = float(sigma)
def create_model(self, **kwargs):
m = dm.Model("langmuir")
qm = m.continuous("qm", lb=0.01, ub=20.0)
K = m.continuous("K", lb=1e-4, ub=1e4)
responses = {f"q_{i}": qm * K * c / (1.0 + K * c) for i, c in enumerate(self.c_data)}
errors = {k: self.sigma for k in responses}
return ExperimentModel(m, {"qm": qm, "K": K}, {}, responses, errors)
qm_true, K_true = 5.0, 1.0
c_sat = np.array([0.1, 0.3, 1.0, 3.0, 10.0]) # spans linear AND saturated regimes
c_lin = np.array([0.001, 0.002, 0.005, 0.01, 0.02]) # linear regime only
# Compare diagnostics in the two regimes (using true values as the
# nominal point so we are not chasing fit noise).
diag_sat = diagnose_identifiability(LangmuirExperiment(c_sat), {"qm": qm_true, "K": K_true})
diag_lin = diagnose_identifiability(LangmuirExperiment(c_lin), {"qm": qm_true, "K": K_true})
print(f"Saturated regime: rank={diag_sat.fim_rank}, cond={diag_sat.condition_number:.3g}")
print(f" VIFs: {diag_sat.vif}")
print()
print(f"Linear regime: rank={diag_lin.fim_rank}, cond={diag_lin.condition_number:.3g}")
print(f" VIFs: {diag_lin.vif}")
Saturated regime: rank=2, cond=10.1
VIFs: {'qm': 2.658123069377644, 'K': 2.6581230693776443}
Linear regime: rank=2, cond=1.06e+06
VIFs: {'qm': 40325.380387478624, 'K': 40325.38038747861}
The linear regime exposes a near-singular FIM (condition number \(\gg 100\)) and inflated VIFs — the classic signature. Now profile likelihood the \(K\) parameter in the well-designed regime to confirm a bounded CI:
np.random.seed(42)
exp_sat = LangmuirExperiment(c_sat)
data_sat = {
f"q_{i}": qm_true * K_true * c / (1 + K_true * c) + np.random.normal(0, 0.05)
for i, c in enumerate(c_sat)
}
est_sat = estimate_parameters(exp_sat, data_sat, initial_guess={"qm": 4.0, "K": 0.5})
print(est_sat.summary())
prof_K = profile_likelihood(exp_sat, data_sat, "K", initial_estimate=est_sat)
prof_qm = profile_likelihood(exp_sat, data_sat, "qm", initial_estimate=est_sat)
fig, (a1, a2) = plt.subplots(1, 2, figsize=(9, 3.5))
for ax, p, name, true in [(a1, prof_K, "K", K_true), (a2, prof_qm, "qm", qm_true)]:
ax.plot(p.theta_values, p.neg_log_lik - p.objective_hat, "o-", markersize=3)
ax.axhline(p.threshold - p.objective_hat, color="red", ls="--", label=r"$\chi^2_{1,0.95}$")
ax.axvline(true, color="green", ls=":", label="truth")
ax.set_title(f"{name} ({p.shape})")
ax.set_xlabel(name)
ax.legend()
a1.set_ylabel(r"$D(\theta) - D(\hat\theta)$")
plt.tight_layout()
Parameter Estimation Results
==================================================
qm = 5.00517 ± 0.06176 [4.809, 5.202]
K = 1.0285 ± 0.04408 [0.8882, 1.169]
Objective = 1.6786
N obs = 5
FIM det = 3.546e+05
FIM cond = 9.65
Worked example 3 — Batch reactor rate data#
A batch reactor obeying first-order kinetics with a non-zero asymptote:
Three parameters: initial concentration \(C_0\), rate constant \(k\), and asymptote \(C_\infty\).
Whether all three are identifiable from a given experiment depends on how the time grid is chosen relative to the system’s relaxation time \(1/k\). If sampling stops before reaching the asymptote, \(C_\infty\) trades off with \(C_0\) and the FIM becomes ill-conditioned [Yao et al., 2003].
class BatchReactorExperiment(Experiment):
def __init__(self, t_data, sigma=0.02):
self.t_data = list(t_data)
self.sigma = float(sigma)
def create_model(self, **kwargs):
m = dm.Model("batch_reactor")
C0 = m.continuous("C0", lb=0.1, ub=10.0)
k = m.continuous("k", lb=0.001, ub=10.0)
C_inf = m.continuous("C_inf", lb=0.0, ub=5.0)
responses = {
f"C_{i}": (C0 - C_inf) * dm.exp(-k * t) + C_inf for i, t in enumerate(self.t_data)
}
errors = {kname: self.sigma for kname in responses}
return ExperimentModel(m, {"C0": C0, "k": k, "C_inf": C_inf}, {}, responses, errors)
C0_true, k_true, Cinf_true = 5.0, 0.5, 1.0
# Long horizon: 5 half-lives -> all three identifiable.
t_long = np.linspace(0.0, 10.0, 11)
# Short horizon: < 1 half-life -> C_inf weakly identifiable.
t_short = np.linspace(0.0, 1.5, 8)
# Diagnose at the truth (no fitting noise).
diag_long = diagnose_identifiability(
BatchReactorExperiment(t_long), {"C0": C0_true, "k": k_true, "C_inf": Cinf_true}
)
diag_short = diagnose_identifiability(
BatchReactorExperiment(t_short), {"C0": C0_true, "k": k_true, "C_inf": Cinf_true}
)
print(f"Long horizon (t<=10): rank={diag_long.fim_rank}/3, cond={diag_long.condition_number:.3g}")
print(f" VIFs: {diag_long.vif}")
print(f" warnings: {diag_long.warnings or '(none)'}")
print()
print(
f"Short horizon (t<=1.5): rank={diag_short.fim_rank}/3, cond={diag_short.condition_number:.3g}"
)
print(f" VIFs: {diag_short.vif}")
print(" warnings:")
for w in diag_short.warnings:
print(f" - {w}")
Long horizon (t<=10): rank=3/3, cond=34.7
VIFs: {'C0': 1.4202105098138078, 'k': 2.9148442990018943, 'C_inf': 2.304553156223133}
warnings: (none)
Short horizon (t<=1.5): rank=3/3, cond=6.97e+03
VIFs: {'C0': 3.264772139993319, 'k': 193.07217844585176, 'C_inf': 172.31262856628257}
warnings:
- serious collinearity: condition index eta_3 = 31.1 > 30
- VIF[k] = 193 > 10
- VIF[C_inf] = 172 > 10
- |rho[k,C_inf]| = 0.994 > 0.95
The short-horizon design has a much larger condition number and one or more inflated VIFs — Yao’s ranking will recommend dropping the worst offender:
rank_short = estimability_rank(
BatchReactorExperiment(t_short),
{"C0": C0_true, "k": k_true, "C_inf": Cinf_true},
)
print("Yao ranking on the short-horizon experiment:")
top = rank_short.projected_norms[0]
for name, pn in zip(rank_short.ranking, rank_short.projected_norms):
bar = "=" * max(1, int(40 * pn / max(top, 1e-30)))
print(f" {name:>6s} |R_kk| = {pn:10.4g} {bar}")
print()
print(f"Recommended subset (cutoff 0.04): {rank_short.recommended_subset}")
print(f"Collinearity index: {rank_short.collinearity_index:.3g}")
Yao ranking on the short-horizon experiment:
C0 |R_kk| = 515.4 ========================================
k |R_kk| = 94.23 =======
C_inf |R_kk| = 3.653 =
Recommended subset (cutoff 0.04): ['C0', 'k']
Collinearity index: 2
Take-away. The three diagnostics
(diagnose_identifiability \(\to\) condition number / VIF;
estimability_rank \(\to\) projected norms;
profile_likelihood \(\to\) shape) form a coherent picture: when the
experiment’s design under-excites a parameter, all three flag it. A
better-designed experiment (here, longer sampling horizon for the
batch reactor; saturated-regime concentrations for Langmuir) restores
identifiability. This is the model-based design-of-experiments loop
in miniature [Villaverde, 2019].
Further reading#
See the Crucible wiki articles behind this notebook:
Parameter identifiability analysis — overview and workflow.
Algebraic model identifiability — the Jacobian-rank recipe.
Parameter estimability — Yao, Brun, Chu & Hahn methods [Chu and Hahn, 2007, Chu and Hahn, 2012].
Profile likelihood identifiability — Raue/Kreutz algorithm.