Identifiability Diagnostics

Installation

install.packages("pak")                   # development version
pak::pak("gcol33/degen")

What is identifiability?

A parameter is identifiable if its value can be uniquely determined from the data. When parameters are non-identifiable, the likelihood function is flat in some direction, meaning infinitely many parameter values produce the same likelihood.

Types of identifiability problems

Structural non-identifiability

Some parameters cannot be identified regardless of sample size:

# Only a + b is identifiable, not a and b separately
nonid_spec <- model_spec(
  loglik_fn = function(y, a, b) {
    sum(dnorm(y, mean = a + b, sd = 1, log = TRUE))
  },
  par_names = c("a", "b"),
  par_bounds = list(a = c(-10, 10), b = c(-10, 10)),
  name = "Non-identifiable sum"
)

Weak identifiability

Parameters may be theoretically identifiable but poorly determined in practice:

# With limited data, variance-covariance can be hard to estimate
bivariate_spec <- model_spec(
  loglik_fn = function(y, mu1, mu2, sigma1, sigma2, rho) {
    # Simplified bivariate normal (ignoring rho for this example)
    ll <- dnorm(y[, 1], mu1, sigma1, log = TRUE) +
          dnorm(y[, 2], mu2, sigma2, log = TRUE)
    sum(ll)
  },
  par_names = c("mu1", "mu2", "sigma1", "sigma2", "rho"),
  par_bounds = list(
    mu1 = c(-10, 10), mu2 = c(-10, 10),
    sigma1 = c(0.01, 10), sigma2 = c(0.01, 10),
    rho = c(-0.99, 0.99)
  ),
  name = "Bivariate normal"
)

Basic identifiability check

identifiability_check() combines Fisher information analysis with profile likelihood diagnostics:

set.seed(42)
y <- rnorm(100, mean = 5, sd = 1)

# Check the non-identifiable model
result <- identifiability_check(nonid_spec, y, par = c(a = 2, b = 3))
print(result)
#> <identifiability_result>
#> Overall: Model has NON-IDENTIFIABLE parameters
#> 
#> Parameter status:
#>   a: NON-IDENTIFIABLE
#>   b: NON-IDENTIFIABLE
#> 
#> Non-identified directions:
#>   -0.71*a + 0.71*b
#> 
#> Identified functions:
#>   0.71*a + 0.71*b
#> 
#> Condition number: 3.1e+12
#> Rank: 1 / 2

Understanding the output

The result shows:

Parameter classifications: Each parameter is classified as identified, weakly identified, or non-identifiable
Condition number: High values (> 10^6) suggest ill-conditioning
Rank deficiency: Indicates how many parameters are non-identifiable

Fisher information analysis

Computing Fisher information

The Fisher information matrix quantifies how much information the data provide about parameters:

info <- fisher_information(nonid_spec, y, par = c(a = 2, b = 3))
print(info)
#> <fisher_info> (observed) at par = (a=2, b=3)
#> Condition number: 3.1e+12
#> Rank: 1 / 2 (RANK DEFICIENT)
#> Eigenvalues: 200, -6.38e-11

Eigenvalue analysis

Eigenvalues reveal identifiability structure:

# Get eigenvalues
info_eigenvalues(info)
#> [1]  2.000000e+02 -6.377121e-11

# Condition number (ratio of largest to smallest eigenvalue)
info_condition(info)
#> [1] 3.136211e+12

# Rank (number of identifiable directions)
info_rank(info)
#> [1] 1

Finding non-identifiable directions

# Direction(s) along which likelihood is flat
null_directions(info)
#>       null_1
#> a -0.7071068
#> b  0.7071068

For the sum model, the null direction is (1, -1), meaning a and b can trade off while keeping a + b constant.

Visualizing identifiability

Eigenvalue spectrum

plot(info)

Profile likelihood

Profile likelihood examines each parameter individually:

# Profile likelihood for parameter 'a'
profile <- profile_likelihood(nonid_spec, y, par = c(a = 2, b = 3),
                              which_par = "a", n_points = 20)
plot_profile_likelihood(profile)

Expected vs. observed information

Observed information

Computed from the Hessian of the log-likelihood at a specific point:

obs_info <- fisher_information(nonid_spec, y, par = c(a = 2, b = 3),
                                type = "observed")

Expected information

Computed as the average outer product of the score across observations:

exp_info <- fisher_information(nonid_spec, y, par = c(a = 2, b = 3),
                                type = "expected")

Practical identifiability workflow

Step 1: Initial check

# Define your model
model <- model_spec(
  loglik_fn = function(y, alpha, beta, gamma) {
    mu <- alpha + beta * (1 - exp(-gamma * seq_along(y)))
    sum(dnorm(y, mean = mu, sd = 1, log = TRUE))
  },
  par_names = c("alpha", "beta", "gamma"),
  par_bounds = list(
    alpha = c(-10, 10),
    beta = c(0.01, 10),
    gamma = c(0.01, 5)
  ),
  name = "Growth model"
)

# Generate test data
set.seed(123)
n <- 50
true_par <- c(alpha = 1, beta = 5, gamma = 0.5)
mu_true <- true_par["alpha"] + true_par["beta"] * (1 - exp(-true_par["gamma"] * (1:n)))
y_growth <- rnorm(n, mean = mu_true, sd = 1)

Step 2: Run diagnostics

# Check identifiability at estimated values
check <- identifiability_check(model, y_growth, par = true_par)
print(check)
#> <identifiability_result>
#> Overall: Model has POORLY IDENTIFIED parameters
#> 
#> Parameter status:
#>   alpha: NON-IDENTIFIABLE
#>   beta: NON-IDENTIFIABLE
#>   gamma: identified
#> 
#> Non-identified directions:
#>   0.71*alpha -0.71*beta
#> 
#> Condition number: 1.4e+03
#> Rank: 3 / 3

Step 3: Examine Fisher information

info <- fisher_information(model, y_growth, par = true_par)
cat("Condition number:", info_condition(info), "\n")
#> Condition number: 1402.097
cat("Eigenvalues:", round(info_eigenvalues(info), 4), "\n")
#> Eigenvalues: 112.4737 60.2475 0.0802

Step 4: Profile likelihood (if needed)

When parameters appear weakly identifiable, profile likelihood provides deeper insight:

# Profile each parameter
for (p in c("alpha", "beta", "gamma")) {
  profile <- profile_likelihood(model, y_growth, par = true_par,
                                which_par = p, n_points = 30)
  plot_profile_likelihood(profile)
}

Addressing identifiability issues

Reparameterization

Transform to identifiable parameters:

# Instead of estimating a and b separately,
# estimate their sum (which IS identifiable)
reparam_spec <- model_spec(
  loglik_fn = function(y, sum_ab) {
    sum(dnorm(y, mean = sum_ab, sd = 1, log = TRUE))
  },
  par_names = "sum_ab",
  par_bounds = list(sum_ab = c(-20, 20)),
  name = "Reparameterized"
)

check_reparam <- identifiability_check(reparam_spec, y, par = c(sum_ab = 5))
print(check_reparam)
#> <identifiability_result>
#> Overall: Model appears WELL-IDENTIFIED
#> 
#> Parameter status:
#>   sum_ab: identified
#> 
#> Condition number: 1
#> Rank: 1 / 1

Fixing parameters

Fix non-identifiable parameters to interpretable values:

# Fix b = 0 and only estimate a
fixed_spec <- model_spec(
  loglik_fn = function(y, a) {
    sum(dnorm(y, mean = a + 0, sd = 1, log = TRUE))  # b fixed at 0
  },
  par_names = "a",
  name = "Fixed b"
)

Adding constraints

Impose constraints based on domain knowledge:

# Constrain a > b (not directly enforceable, but useful conceptually)
# Often implemented via transformation
constrained_spec <- model_spec(
  loglik_fn = function(y, a, diff) {
    b <- a - exp(diff)  # b = a - exp(diff), so b < a always
    sum(dnorm(y, mean = a + b, sd = 1, log = TRUE))
  },
  par_names = c("a", "diff"),
  name = "Constrained"
)

Model diagnostics

Full diagnostic report

diag <- diagnose_model(model, y_growth, par = true_par)
print(diag)

Gradient checking

Verify that numerical derivatives are accurate:

check_gradient(model, y_growth, par = true_par)

Interpreting condition numbers

Condition Number	Interpretation
< 100	Well-conditioned, all parameters identifiable
100 - 10^4	Moderate conditioning, some weak identifiability
10^4 - 10^8	Poorly conditioned, likely identifiability issues
> 10^8	Severely ill-conditioned, structural non-identifiability

Best practices

Always check identifiability before drawing conclusions from fitted parameters
Use multiple starting points to detect if optimization finds different solutions
Examine profile likelihoods for parameters with high uncertainty
Consider the condition number as a summary diagnostic
Report null directions when non-identifiability is detected - they often have meaningful interpretation

Session Info

sessionInfo()
#> R version 4.6.1 (2026-06-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26300)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: Europe/Luxembourg
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] degen_0.15.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] digest_0.6.39       desc_1.4.3          R6_2.6.1           
#>  [4] numDeriv_2016.8-1.1 fastmap_1.2.0       xfun_0.60          
#>  [7] cachem_1.1.0        knitr_1.52          htmltools_0.5.9    
#> [10] rmarkdown_2.32      lifecycle_1.0.5     cli_3.6.6          
#> [13] svglite_2.2.2       sass_0.4.10         pkgdown_2.2.1      
#> [16] textshaping_1.0.5   jquerylib_0.1.4     systemfonts_1.3.2  
#> [19] compiler_4.6.1      tools_4.6.1         bslib_0.12.0       
#> [22] evaluate_1.0.5      Rcpp_1.1.2          yaml_2.3.12        
#> [25] otel_0.2.0          jsonlite_2.0.0      rlang_1.3.0        
#> [28] fs_2.1.0            htmlwidgets_1.6.4

Gilles Colling

2026-11-02