Definition

Consider a system of \(m\) independent components. The lifetime of component \(j\) in the \(i\)-th system is \(T_{ij}\), with reliability (survival) function \(R_j(t) = P(T_{ij} > t)\), pdf \(f_j(t)\), and hazard function \[ h_j(t) = \frac{f_j(t)}{R_j(t)}. \]

In a series configuration (“weakest link”), the system fails as soon as any component fails: \[ T_i = \min(T_{i1}, \ldots, T_{im}). \]

Assumption (independence). Component lifetimes \(T_{i1}, \ldots, T_{im}\) are mutually independent given the parameter vector \(\theta\).

System reliability

Additive hazards

The additive hazard decomposition is the central structural result for series systems. It means that each component contributes independently to the instantaneous system failure risk at every time \(t\).

The system pdf follows immediately: \[ f(t \mid \theta) = h(t \mid \theta) \cdot R(t \mid \theta) = \left[\sum_{j=1}^m h_j(t)\right] \prod_{l=1}^m R_l(t). \]

model_exp <- exp_series_md_c1_c2_c3()

# Component failure rates
theta_exp <- c(1.0, 1.1, 0.95)
m <- length(theta_exp)

# Extract hazard closures
h1 <- component_hazard(model_exp, 1)
h2 <- component_hazard(model_exp, 2)
h3 <- component_hazard(model_exp, 3)

t_grid <- seq(0.01, 3, length.out = 200)
h_sys <- h1(t_grid, theta_exp) + h2(t_grid, theta_exp) + h3(t_grid, theta_exp)

cat("System hazard (constant):", h_sys[1], "\n")
#> System hazard (constant): 3.05
cat("Sum of component rates:  ", sum(theta_exp), "\n")
#> Sum of component rates:   3.05

Definition

Let \(K_i\) denote the index of the component that caused the \(i\)-th system failure: \[ K_i = \arg\min_{j \in \{1,\ldots,m\}} T_{ij}. \]

\(K_i\) is a discrete random variable on \(\{1, \ldots, m\}\) that, in general, depends on the failure time \(T_i\).

Conditional cause probability

Intuition. At any instant \(t\), the component with the highest hazard is the most likely cause. The cause probability is the component’s share of the total instantaneous risk — a proportional allocation.

# Exponential: constant cause probabilities
ccp_exp <- conditional_cause_probability(model_exp)
probs_exp <- ccp_exp(t_grid, theta_exp)

# Heterogeneous Weibull: time-varying
theta_wei <- c(0.7, 200,   # DFR electronics
               1.0, 150,   # CFR seals
               2.0, 100)   # IFR bearing
model_wei <- wei_series_md_c1_c2_c3(
  shapes = theta_wei[seq(1, 6, 2)],
  scales = theta_wei[seq(2, 6, 2)]
)
ccp_wei <- conditional_cause_probability(model_wei)

t_wei <- seq(0.5, 120, length.out = 200)
probs_wei <- ccp_wei(t_wei, theta_wei)

oldpar <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
on.exit(par(oldpar))

# Exponential
plot(t_grid, probs_exp[, 1], type = "l", col = "steelblue", lwd = 2,
     ylim = c(0, 0.6), xlab = "Time", ylab = "P(K = j | T = t)",
     main = "Exponential (constant)")
lines(t_grid, probs_exp[, 2], col = "forestgreen", lwd = 2, lty = 2)
lines(t_grid, probs_exp[, 3], col = "firebrick", lwd = 2, lty = 3)
legend("right", paste("Comp", 1:3),
       col = c("steelblue", "forestgreen", "firebrick"),
       lty = 1:3, lwd = 2, cex = 0.8)

# Heterogeneous Weibull
plot(t_wei, probs_wei[, 1], type = "l", col = "steelblue", lwd = 2,
     ylim = c(0, 1), xlab = "Time", ylab = "P(K = j | T = t)",
     main = "Heterogeneous Weibull (time-varying)")
lines(t_wei, probs_wei[, 2], col = "forestgreen", lwd = 2, lty = 2)
lines(t_wei, probs_wei[, 3], col = "firebrick", lwd = 2, lty = 3)
legend("right", c("DFR (k=0.7)", "CFR (k=1.0)", "IFR (k=2.0)"),
       col = c("steelblue", "forestgreen", "firebrick"),
       lty = 1:3, lwd = 2, cex = 0.8)

Four observation types

The package supports four types, each arising from a different monitoring scheme:

Here \(h_c(t) = \sum_{j \in c_i} h_j(t)\) is the candidate-set hazard.

Masking

When \(|c_i| = 1\), the failed component is known exactly. When \(|c_i| > 1\), the cause of failure is masked: the diagnostic identifies a set of possible culprits but not the specific one.

Observe functors

The package provides composable observation functors for generating data under different monitoring schemes:

# Continuous monitoring with right-censoring at tau
obs_right <- observe_right_censor(tau = 5)
obs_right(3.2)   # failure before tau -> exact
#> $t
#> [1] 3.2
#> 
#> $omega
#> [1] "exact"
#> 
#> $t_upper
#> [1] NA
obs_right(7.1)   # survival past tau  -> right-censored
#> $t
#> [1] 5
#> 
#> $omega
#> [1] "right"
#> 
#> $t_upper
#> [1] NA

# Single inspection at tau (left-censoring)
obs_left <- observe_left_censor(tau = 5)
obs_left(3.2)    # failed before inspection -> left-censored
#> $t
#> [1] 5
#> 
#> $omega
#> [1] "left"
#> 
#> $t_upper
#> [1] NA
obs_left(7.1)    # surviving at inspection  -> right-censored
#> $t
#> [1] 5
#> 
#> $omega
#> [1] "right"
#> 
#> $t_upper
#> [1] NA

# Periodic inspections every delta time units
obs_periodic <- observe_periodic(delta = 1, tau = 5)
obs_periodic(3.2)  # failure in (3, 4) -> interval-censored
#> $t
#> [1] 3
#> 
#> $omega
#> [1] "interval"
#> 
#> $t_upper
#> [1] 4

# Mixture: random assignment to schemes
obs_mixed <- observe_mixture(
  observe_right_censor(tau = 5),
  observe_left_censor(tau = 3),
  observe_periodic(delta = 0.5, tau = 5),
  weights = c(0.5, 0.2, 0.3)
)

All rdata() methods accept an observe argument to generate data under a specified monitoring scheme:

gen_exp <- rdata(model_exp)

set.seed(42)
df_demo <- gen_exp(theta_exp, n = 500, p = 0.3,
                   observe = observe_mixture(
                     observe_right_censor(tau = 5),
                     observe_left_censor(tau = 3),
                     weights = c(0.7, 0.3)
                   ))
cat("Observation type counts:\n")
#> Observation type counts:
print(table(df_demo$omega))
#> 
#> exact  left 
#>   344   156

The three conditions

C1 (containment). The candidate set always contains the true failed component: \[ P(K_i \in c_i) = 1. \] This is a minimal accuracy requirement on the diagnostic. Without C1, the candidate set could be misleading, and the likelihood would need to model diagnostic errors explicitly.

Example: An on-board diagnostic (OBD) system reports a set of possibly-failed modules. C1 requires that the truly failed module is always in this set — the OBD may include false positives but never a false negative.

C2 (symmetric masking). The masking probability is the same for all components in the candidate set: \[ P(c_i \mid K_i = j, T_i = t_i, \theta) = P(c_i \mid K_i = j', T_i = t_i, \theta) \quad \text{for all } j, j' \in c_i. \] Intuitively, the diagnostic cannot distinguish among the candidates — if it could, it would narrow the candidate set.

Example: A line-replaceable unit (LRU) grouping replaces all components in a subsystem together. The grouping is determined by physical layout, not by which component failed, so C2 holds by construction.

C3 (parameter independence). The masking probabilities do not depend on the lifetime parameters \(\theta\): \[ P(c_i \mid K_i = j, T_i = t_i) \text{ is not a function of } \theta. \]

Example: Diagnostic accuracy depends on sensor calibration and technician skill, not on component failure rates. If the diagnostic equipment was calibrated before the study, C3 holds.

Deriving the likelihood

We derive the likelihood step by step. Consider observation \(i\) with failure time \(t_i\) and candidate set \(c_i\).

Step 1: Full joint density. The observed-data density includes a sum over the latent cause \(K_i\): \[ L_i(\theta) = \sum_{j=1}^m f_{K,T}(j, t_i) \cdot P(c_i \mid K_i = j, T_i = t_i, \theta). \]

Step 2: Apply C1. Since \(P(K_i \in c_i) = 1\), contributions from \(j \notin c_i\) are zero: \[ L_i(\theta) = \sum_{j \in c_i} f_{K,T}(j, t_i) \cdot P(c_i \mid K_i = j, T_i = t_i, \theta). \]

Step 3: Apply C2. Symmetric masking means \(P(c_i \mid K_i = j, T_i = t_i, \theta) = \beta_i\) for all \(j \in c_i\): \[ L_i(\theta) = \beta_i \sum_{j \in c_i} f_{K,T}(j, t_i). \]

Step 4: Apply C3. Since \(\beta_i\) does not depend on \(\theta\), it is a constant with respect to the parameters and drops from the maximization: \[ L_i(\theta) \propto \sum_{j \in c_i} f_{K,T}(j, t_i) = \sum_{j \in c_i} h_j(t_i) \cdot R(t_i) = R(t_i) \cdot h_c(t_i). \]

The exact-failure likelihood contribution is therefore: \[\begin{equation} \label{eq:exact_contrib} \ell_i^{\text{exact}}(\theta) = \log R(t_i \mid \theta) + \log h_c(t_i \mid \theta), \end{equation}\] where \(h_c(t) = \sum_{j \in c_i} h_j(t)\) is the candidate-set hazard.

Censored observations

Right-censored. The system survived past \(t_i\). No failure was observed, so there is no candidate set and the contribution is simply: \[ \ell_i^{\text{right}}(\theta) = \log R(t_i \mid \theta). \]

Left-censored. The system was found failed at inspection time \(\tau_i\), but the exact failure time is unknown. Applying C1–C3 to the integrated likelihood: \[ \ell_i^{\text{left}}(\theta) = \log \int_0^{\tau_i} h_c(u) \, R(u) \, du. \]

Interval-censored. The failure occurred in \((a_i, b_i)\): \[ \ell_i^{\text{interval}}(\theta) = \log \int_{a_i}^{b_i} h_c(u) \, R(u) \, du. \]

Combined log-likelihood

The full log-likelihood decomposes as: \[\begin{equation} \label{eq:full_loglik} \ell(\theta) = \ell_E + \ell_R + \ell_L + \ell_I, \end{equation}\] where \[\begin{align} \ell_E &= \sum_{i : \omega_i = \text{exact}} \bigl[\log R(t_i) + \log h_{c_i}(t_i)\bigr], \\ \ell_R &= \sum_{i : \omega_i = \text{right}} \log R(t_i), \\ \ell_L &= \sum_{i : \omega_i = \text{left}} \log \int_0^{\tau_i} h_{c_i}(u) R(u) \, du, \\ \ell_I &= \sum_{i : \omega_i = \text{interval}} \log \int_{a_i}^{b_i} h_{c_i}(u) R(u) \, du. \end{align}\]

ll_exp <- loglik(model_exp)

# Evaluate at true parameters on the mixed-censoring demo data
ll_val <- ll_exp(df_demo, theta_exp)
cat("Log-likelihood at true theta:", round(ll_val, 4), "\n")
#> Log-likelihood at true theta: -345.3

# Perturb parameters and compare
theta_bad <- c(2, 2, 2)
cat("Log-likelihood at wrong theta:", round(ll_exp(df_demo, theta_bad), 4), "\n")
#> Log-likelihood at wrong theta: -470.8

Setup

theta <- c(1.0, 1.1, 0.95, 1.15, 1.1)
m <- length(theta)
model <- exp_series_md_c1_c2_c3()

cat("True rates:", theta, "\n")
#> True rates: 1 1.1 0.95 1.15 1.1
cat("System rate:", sum(theta), "\n")
#> System rate: 5.3

Specializing the likelihood

For exponential components with \(h_j(t) = \lambda_j\) and \(R_j(t) = e^{-\lambda_j t}\):

\[\begin{align} \text{Exact:} \quad & \log\!\left(\sum_{j \in c_i} \lambda_j\right) - \left(\sum_{j=1}^m \lambda_j\right) t_i, \\ \text{Right:} \quad & -\lambda_{\text{sys}} t_i, \\ \text{Left:} \quad & \log \lambda_c + \log(1 - e^{-\lambda_{\text{sys}} \tau_i}) - \log \lambda_{\text{sys}}, \\ \text{Interval:} \quad & \log \lambda_c - \lambda_{\text{sys}} a_i + \log(1 - e^{-\lambda_{\text{sys}}(b_i - a_i)}) - \log \lambda_{\text{sys}}, \end{align}\] where \(\lambda_c = \sum_{j \in c_i} \lambda_j\) and \(\lambda_{\text{sys}} = \sum_j \lambda_j\).

All four observation types have closed-form log-likelihood, score, and Hessian — the exponential model is unique in this respect.

Data generation and fitting

gen <- rdata(model)
set.seed(7231)

df <- gen(theta, n = 2000, p = 0.3,
          observe = observe_right_censor(tau = 3))

cat("Observation types:\n")
#> Observation types:
print(table(df$omega))
#> 
#> exact 
#>  2000

solver <- fit(model)
theta0 <- rep(1, m)
estimate <- solver(df, par = theta0, method = "Nelder-Mead")

recovery <- data.frame(
  Component = 1:m,
  True = theta,
  MLE = estimate$par,
  SE = sqrt(diag(estimate$vcov)),
  Rel_Error_Pct = 100 * (estimate$par - theta) / theta
)

knitr::kable(recovery, digits = 4,
             caption = "Exponential MLE: parameter recovery",
             col.names = c("Component", "True", "MLE", "SE", "Rel. Error %"))

Score and Hessian verification

The analytical score should vanish at the MLE. We also verify it against numerical differentiation:

ll_fn <- loglik(model)
scr_fn <- score(model)
hess_fn <- hess_loglik(model)

scr_at_mle <- scr_fn(df, estimate$par)
cat("Score at MLE:", round(scr_at_mle, 4), "\n")
#> Score at MLE: -0.0494 0.1326 -0.0488 0.0388 0.1303

scr_num <- numDeriv::grad(function(th) ll_fn(df, th), estimate$par)
cat("Max |analytical - numerical| score:",
    formatC(max(abs(scr_at_mle - scr_num)), format = "e", digits = 2), "\n")
#> Max |analytical - numerical| score: 1.14e-06

# Hessian eigenvalues (should all be negative for concavity)
H <- hess_fn(df, estimate$par)
cat("Hessian eigenvalues:", round(eigen(H)$values, 2), "\n")
#> Hessian eigenvalues: -142.6 -153.3 -167.4 -174.5 -370.4

Fisher information and confidence intervals

The observed Fisher information is \(I(\hat\theta) = -H(\hat\theta)\). Asymptotic \((1-\alpha)\)% Wald confidence intervals are: \[ \hat\theta_j \pm z_{1-\alpha/2} \sqrt{I(\hat\theta)^{-1}_{jj}}. \]

alpha <- 0.05
z <- qnorm(1 - alpha / 2)
se <- sqrt(diag(estimate$vcov))

ci_df <- data.frame(
  Component = 1:m,
  Lower = estimate$par - z * se,
  MLE = estimate$par,
  Upper = estimate$par + z * se,
  True = theta,
  Covered = (estimate$par - z * se <= theta) & (theta <= estimate$par + z * se)
)

knitr::kable(ci_df, digits = 4,
             caption = "95% Wald confidence intervals",
             col.names = c("Comp.", "Lower", "MLE", "Upper", "True", "Covered?"))

Key properties

Property 1: The system lifetime is Weibull. Because all shapes are equal, the system reliability simplifies: \[ R(t) = \prod_{j=1}^m \exp\!\left(-\left(\frac{t}{\beta_j}\right)^k\right) = \exp\!\left(-\left(\frac{t}{\beta_{\text{sys}}}\right)^k\right), \] where \(\beta_{\text{sys}} = \bigl(\sum_{j=1}^m \beta_j^{-k}\bigr)^{-1/k}\). This closure under the minimum is the structural advantage of the homogeneous model.

Property 2: Constant cause probabilities. The time dependence \(t^{k-1}\) cancels in the hazard ratio: \[ P(K = j \mid T = t) = \frac{\beta_j^{-k}}{\sum_l \beta_l^{-k}} \eqqcolon w_j. \] Just as in the exponential case, the conditional cause probability does not depend on \(t\).

Property 3: Closed-form censored contributions. Left- and interval-censored terms factor as \(\log w_{c_i} + \log [\text{Weibull CDF difference}]\) — no numerical integration needed.

R code: Setup and hazard visualization

theta_hom <- c(k = 1.5, beta1 = 100, beta2 = 150, beta3 = 200)
k <- theta_hom[1]
scales <- theta_hom[-1]
m_hom <- length(scales)

model_hom <- wei_series_homogeneous_md_c1_c2_c3()

# System scale
beta_sys <- wei_series_system_scale(k, scales)
cat("System scale:", round(beta_sys, 2), "\n")
#> System scale: 65.24

# Plot component hazards
t_grid <- seq(0.1, 200, length.out = 300)
cols <- c("steelblue", "forestgreen", "firebrick")

plot(NULL, xlim = c(0, 200), ylim = c(0, 0.06),
     xlab = "Time", ylab = "Hazard h_j(t)",
     main = "Homogeneous Weibull: component hazards (k = 1.5)")
for (j in seq_len(m_hom)) {
  h_j <- component_hazard(model_hom, j)
  lines(t_grid, h_j(t_grid, theta_hom), col = cols[j], lwd = 2)
}
legend("topleft", paste0("Comp ", 1:m_hom, " (beta=", scales, ")"),
       col = cols, lwd = 2, cex = 0.9)
grid()

All hazards are increasing (\(k > 1\)), sharing the same shape but differing in magnitude. Component 1 (smallest \(\beta\)) dominates at all times.

Cause probabilities

# Analytical cause weights
w <- scales^(-k) / sum(scales^(-k))
names(w) <- paste0("Component ", 1:m_hom)
cat("Time-invariant cause weights:\n")
#> Time-invariant cause weights:
print(round(w, 4))
#> Component 1 Component 2 Component 3 
#>      0.5269      0.2868      0.1863

# Verify with package function (pass scales so model knows m)
ccp_fn <- conditional_cause_probability(
  wei_series_homogeneous_md_c1_c2_c3(scales = scales)
)
probs <- ccp_fn(c(10, 50, 100, 150), theta_hom)
knitr::kable(probs, digits = 4,
             caption = "P(K=j | T=t) at four time points (should be constant)",
             col.names = paste0("Comp ", 1:m_hom))

Data generation and MLE

gen_hom <- rdata(model_hom)
set.seed(2024)

df_hom <- gen_hom(theta_hom, n = 500, p = 0.3,
                  observe = observe_periodic(delta = 20, tau = 250))

cat("Observation types:\n")
#> Observation types:
print(table(df_hom$omega))
#> 
#> interval 
#>      500

solver_hom <- fit(model_hom)
theta0_hom <- c(1.2, 110, 140, 180)

est_hom <- solver_hom(df_hom, par = theta0_hom, method = "Nelder-Mead")

recovery_hom <- data.frame(
  Parameter = c("k", paste0("beta_", 1:m_hom)),
  True = theta_hom,
  MLE = est_hom$par,
  SE = sqrt(diag(est_hom$vcov)),
  Rel_Error_Pct = 100 * (est_hom$par - theta_hom) / theta_hom
)

knitr::kable(recovery_hom, digits = 3,
             caption = paste0("Homogeneous Weibull MLE (n = 500)"),
             col.names = c("Parameter", "True", "MLE", "SE", "Rel. Error %"))

Differences from the homogeneous model

The additional flexibility comes at a cost: the system lifetime \(T = \min(T_1, \ldots, T_m)\) is no longer Weibull when shapes differ, and left-censored or interval-censored likelihood contributions require numerical integration (stats::integrate).

R code: Mixed failure regimes

We construct a 3-component system with DFR (infant mortality), CFR (random), and IFR (wear-out) components:

theta_het <- c(0.7, 200,   # DFR electronics
               1.0, 150,   # CFR seals
               2.0, 100)   # IFR bearing
m_het <- length(theta_het) / 2

model_het <- wei_series_md_c1_c2_c3(
  shapes = theta_het[seq(1, 6, 2)],
  scales = theta_het[seq(2, 6, 2)]
)

# Component hazard functions
t_het <- seq(0.5, 120, length.out = 300)

h_fns <- lapply(1:m_het, function(j) component_hazard(model_het, j))
h_vals <- sapply(h_fns, function(hf) hf(t_het, theta_het))

plot(t_het, h_vals[, 1], type = "l", col = "steelblue", lwd = 2,
     ylim = c(0, 0.05),
     xlab = "Time", ylab = "Hazard h_j(t)",
     main = "Heterogeneous Weibull: mixed failure regimes")
lines(t_het, h_vals[, 2], col = "forestgreen", lwd = 2, lty = 2)
lines(t_het, h_vals[, 3], col = "firebrick", lwd = 2, lty = 3)
legend("topright",
       c("Comp 1: DFR (k=0.7)", "Comp 2: CFR (k=1.0)", "Comp 3: IFR (k=2.0)"),
       col = c("steelblue", "forestgreen", "firebrick"),
       lty = 1:3, lwd = 2, cex = 0.85)
grid()

Time-varying cause probabilities

ccp_het <- conditional_cause_probability(model_het)
probs_het <- ccp_het(t_het, theta_het)

plot(t_het, probs_het[, 1], type = "l", col = "steelblue", lwd = 2,
     ylim = c(0, 1), xlab = "Time", ylab = "P(K = j | T = t)",
     main = "Heterogeneous Weibull: time-varying cause probabilities")
lines(t_het, probs_het[, 2], col = "forestgreen", lwd = 2, lty = 2)
lines(t_het, probs_het[, 3], col = "firebrick", lwd = 2, lty = 3)
legend("right",
       c("DFR (k=0.7)", "CFR (k=1.0)", "IFR (k=2.0)"),
       col = c("steelblue", "forestgreen", "firebrick"),
       lty = 1:3, lwd = 2, cex = 0.85)
grid()

The crossover pattern is characteristic of bathtub-curve behavior at the system level: early failures are dominated by infant-mortality mechanisms, while late failures are dominated by wear-out.

Data generation and MLE

gen_het <- rdata(model_het)
set.seed(7231)

df_het <- gen_het(theta_het, n = 300, tau = 120, p = 0.3,
                  observe = observe_right_censor(tau = 120))

cat("Observation types:\n")
#> Observation types:
print(table(df_het$omega))
#> 
#> exact right 
#>   286    14

We use right-censoring for speed — exact and right-censored observations have closed-form likelihood contributions. Left- and interval-censored observations require numerical integration per observation per optimization iteration.

solver_het <- fit(model_het)
theta0_het <- rep(c(1, 150), m_het)

est_het <- solver_het(df_het, par = theta0_het, method = "Nelder-Mead")

par_names <- paste0(
  rep(c("k", "beta"), m_het), "_",
  rep(1:m_het, each = 2)
)

recovery_het <- data.frame(
  Parameter = par_names,
  True = theta_het,
  MLE = est_het$par,
  SE = sqrt(diag(est_het$vcov)),
  Rel_Error_Pct = 100 * (est_het$par - theta_het) / theta_het
)

knitr::kable(recovery_het, digits = 3,
             caption = paste0("Heterogeneous Weibull MLE (n = 300)"),
             col.names = c("Parameter", "True", "MLE", "SE", "Rel. Error %"))

Model comparison: heterogeneous vs homogeneous

When the true DGP has heterogeneous shapes, the homogeneous model is misspecified. We compare the two on the same data:

# Refit on larger dataset for sharper comparison
set.seed(999)
df_comp <- gen_het(theta_het, n = 800, tau = 150, p = 0.2,
                   observe = observe_right_censor(tau = 150))

# Heterogeneous (6 params)
fit_het <- solver_het(df_comp, par = theta0_het, method = "Nelder-Mead")

# Homogeneous (4 params)
model_hom2 <- wei_series_homogeneous_md_c1_c2_c3()
solver_hom2 <- fit(model_hom2)
fit_hom2 <- solver_hom2(df_comp, par = c(1, 150, 150, 150),
                        method = "Nelder-Mead")

# Likelihood ratio test
LRT <- 2 * (fit_het$loglik - fit_hom2$loglik)
df_lrt <- 6 - 4
p_val <- pchisq(LRT, df = df_lrt, lower.tail = FALSE)

comp_df <- data.frame(
  Model = c("Heterogeneous (2m = 6)", "Homogeneous (m+1 = 4)"),
  Params = c(6, 4),
  LogLik = c(fit_het$loglik, fit_hom2$loglik),
  AIC = c(-2 * fit_het$loglik + 2 * 6, -2 * fit_hom2$loglik + 2 * 4)
)

knitr::kable(comp_df, digits = 2,
             caption = "Model comparison",
             col.names = c("Model", "# Params", "Log-Lik", "AIC"))

cat(sprintf("\nLRT statistic: %.2f (df = %d, p = %.4f)\n",
            LRT, df_lrt, p_val))
#> 
#> LRT statistic: 155.88 (df = 2, p = 0.0000)

Both AIC and the likelihood ratio test favor the heterogeneous model when the true DGP has unequal shapes. The common-shape estimate from the homogeneous fit represents a compromise value between 0.7, 1.0, and 2.0, leading to biased estimates of the scale parameters.