0. Outline

This document estimates the incubation period of HFRS (Hemorrhagic Fever with Renal Syndrome) caused by Seoul virus, using data from Kulagin et al. (1962).

Section 2 fits three parametric distributions via MLE and selects the best-fitting family by AIC.

Section 3 is a Bayesian Estimation using Stan with HMC/NUTS, providing posterior distributions for parameters and a posterior predictive check.

1. Dataset

Incubation periods were recorded from 18 laboratory workers who visited the infected department exactly once, providing point observations with known exposure date.

Source: Kulagin, Fedorova & Ketiladze (1962), Journal of Microbiology, Epidemiology and Immunology, 33(10), 121–126. https://www.cabidigitallibrary.org/doi/full/10.5555/19632701583

inc <- c(rep(17, 3), 19, rep(21, 2), rep(22, 8), 23, rep(25, 2), 27)
cat("n =", length(inc), "\n")
## n = 18
cat("Range:", min(inc), "-", max(inc), "days\n")
## Range: 17 - 27 days
cat("Mean:", round(mean(inc), 1), "  Median:", median(inc), "\n")
## Mean: 21.6   Median: 22
ggplot(data.frame(x = inc), aes(x = "", y = x)) +
  geom_boxplot(fill = "#009E73", alpha = 0.4, width = 0.4, outlier.shape = NA) +
  geom_jitter(width = 0.08, alpha = 0.6, color = "#009E73", size = 2) +
  labs(title = "Incubation Period Distribution (n = 18)",
       x = NULL, y = "days") +
  theme_minimal()

The data are integer-valued and notably concentrated: 8 of 18 observations fall on day 22. This discrete, near-symmetric clustering shapes both the MLE and Bayesian results.

2. Maximum Likelihood Estimation

2.1 Fitting

Three parametric families (Lognormal, Weibull, Gamma) commonly used for incubation period modelling are fitted via MLE with L-BFGS-B optimization. Confidence intervals are derived from the observed Fisher information (Hessian inverse).

ll_lnorm <- function(params, data) {
  mu <- params[1]; sigma <- params[2]
  if (sigma <= 0) return(Inf)
  -sum(dlnorm(data, meanlog = mu, sdlog = sigma, log = TRUE))
}
fit_lnorm <- optim(c(mean(log(inc)), sd(log(inc))), ll_lnorm, data = inc,
                   method = "L-BFGS-B", lower = c(-Inf, 0.001), hessian = TRUE)
se_lnorm  <- sqrt(diag(solve(fit_lnorm$hessian)))
ci_lnorm  <- cbind(fit_lnorm$par - 1.96 * se_lnorm,
                   fit_lnorm$par + 1.96 * se_lnorm)
aic_lnorm <- 2 * 2 + 2 * fit_lnorm$value

ll_weib <- function(params, data) {
  shape <- params[1]; scale <- params[2]
  if (shape <= 0 || scale <= 0) return(Inf)
  -sum(dweibull(data, shape = shape, scale = scale, log = TRUE))
}
fit_weib <- optim(c(1, mean(inc)), ll_weib, data = inc,
                  method = "L-BFGS-B", lower = c(0.001, 0.001), hessian = TRUE)
se_weib  <- sqrt(diag(solve(fit_weib$hessian)))
ci_weib  <- cbind(fit_weib$par - 1.96 * se_weib,
                  fit_weib$par + 1.96 * se_weib)
aic_weib <- 2 * 2 + 2 * fit_weib$value

ll_gamma <- function(params, data) {
  shape <- params[1]; rate <- params[2]
  if (shape <= 0 || rate <= 0) return(Inf)
  -sum(dgamma(data, shape = shape, rate = rate, log = TRUE))
}
fit_gamma <- optim(c(1, 1 / mean(inc)), ll_gamma, data = inc,
                   method = "L-BFGS-B", lower = c(0.001, 0.001), hessian = TRUE)
se_gamma  <- sqrt(diag(solve(fit_gamma$hessian)))
ci_gamma  <- cbind(fit_gamma$par - 1.96 * se_gamma,
                   fit_gamma$par + 1.96 * se_gamma)
aic_gamma <- 2 * 2 + 2 * fit_gamma$value

2.2 Parameter Estimates

result <- data.frame(
  Distribution = c("Lognormal", "Weibull", "Gamma"),
  Parameters   = c("meanlog, sdlog", "shape, scale", "shape, rate"),
  Param1       = round(c(fit_lnorm$par[1], fit_weib$par[1], fit_gamma$par[1]), 4),
  CI1          = paste0("[", round(c(ci_lnorm[1,1], ci_weib[1,1], ci_gamma[1,1]), 3),
                         ", ", round(c(ci_lnorm[1,2], ci_weib[1,2], ci_gamma[1,2]), 3), "]"),
  Param2       = round(c(fit_lnorm$par[2], fit_weib$par[2], fit_gamma$par[2]), 4),
  CI2          = paste0("[", round(c(ci_lnorm[2,1], ci_weib[2,1], ci_gamma[2,1]), 3),
                         ", ", round(c(ci_lnorm[2,2], ci_weib[2,2], ci_gamma[2,2]), 3), "]"),
  AIC          = round(c(aic_lnorm, aic_weib, aic_gamma), 2)
)
kable(result, align = "llrrrrc",
      col.names = c("Distribution", "Parameters", "Param 1", "95% CI",
                    "Param 2", "95% CI", "AIC"))
Distribution Parameters Param 1 95% CI Param 2 95% CI AIC
Lognormal meanlog, sdlog 3.0628 [3.004, 3.121] 0.1264 [0.085, 0.168] 90.88
Weibull shape, scale 8.9965 [5.89, 12.103] 22.7147 [21.482, 23.948] 90.67
Gamma shape, rate 64.0606 [22.316, 105.805] 2.9719 [1.028, 4.916] 90.56

Interpretation: The model with the lowest AIC is preferred. A difference of ΔAIC < 2 indicates the models are statistically indistinguishable. With n = 18, all three distributions are likely to produce similar AIC values — in that case, lognormal is the conventional choice for incubation periods in epidemiology (multiplicative mechanism; consistent with prior literature on HFRS and other respiratory/zoonotic infections).

The CIs are Wald-type intervals derived from the Hessian; they are symmetric on the parameter scale and may be unreliable if the posterior is skewed — a limitation addressed by the Stan section below.

2.3 Distribution Curves (MLE)

freq <- seq(0, 50, length.out = 1000)

df_lnorm <- data.frame(
  x        = freq,
  density  = dlnorm(freq,  fit_lnorm$par[1], fit_lnorm$par[2]),
  ci_lower = dlnorm(freq,  ci_lnorm[1, 1],   ci_lnorm[2, 1]),
  ci_upper = dlnorm(freq,  ci_lnorm[1, 2],   ci_lnorm[2, 2])
)
plnorm <- ggplot(df_lnorm, aes(x = x)) +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "#009E73", alpha = 0.2) +
  geom_line(aes(y = density), color = "#009E73", linewidth = 1) +
  geom_rug(data = data.frame(x = inc), aes(x = x), inherit.aes = FALSE,
           sides = "b", alpha = 0.6) +
  labs(title = paste0("Lognormal  (AIC = ", round(aic_lnorm, 1), ")"),
       x = "days", y = "density") +
  theme_minimal()

df_weib <- data.frame(
  x        = freq,
  density  = dweibull(freq, fit_weib$par[1],  fit_weib$par[2]),
  ci_lower = dweibull(freq, ci_weib[1, 1],    ci_weib[2, 1]),
  ci_upper = dweibull(freq, ci_weib[1, 2],    ci_weib[2, 2])
)
pweib <- ggplot(df_weib, aes(x = x)) +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "#56B4E9", alpha = 0.2) +
  geom_line(aes(y = density), color = "#56B4E9", linewidth = 1) +
  geom_rug(data = data.frame(x = inc), aes(x = x), inherit.aes = FALSE,
           sides = "b", alpha = 0.6) +
  labs(title = paste0("Weibull  (AIC = ", round(aic_weib, 1), ")"),
       x = "days", y = "density") +
  theme_minimal()

df_gamma <- data.frame(
  x        = freq,
  density  = dgamma(freq,  fit_gamma$par[1], fit_gamma$par[2]),
  ci_lower = dgamma(freq,  ci_gamma[1, 1],   ci_gamma[2, 1]),
  ci_upper = dgamma(freq,  ci_gamma[1, 2],   ci_gamma[2, 2])
)
pgamma_p <- ggplot(df_gamma, aes(x = x)) +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "#E69F00", alpha = 0.2) +
  geom_line(aes(y = density), color = "#E69F00", linewidth = 1) +
  geom_rug(data = data.frame(x = inc), aes(x = x), inherit.aes = FALSE,
           sides = "b", alpha = 0.6) +
  labs(title = paste0("Gamma  (AIC = ", round(aic_gamma, 1), ")"),
       x = "days", y = "density") +
  theme_minimal()

df_combined <- data.frame(
  x            = rep(freq, 3),
  density      = c(dlnorm(freq,   fit_lnorm$par[1], fit_lnorm$par[2]),
                   dweibull(freq, fit_weib$par[1],  fit_weib$par[2]),
                   dgamma(freq,   fit_gamma$par[1], fit_gamma$par[2])),
  Distribution = rep(c("Lognormal", "Weibull", "Gamma"), each = length(freq))
)
pcom <- ggplot(df_combined, aes(x = x, y = density, color = Distribution)) +
  geom_line(linewidth = 1) +
  geom_rug(data = data.frame(x = inc), aes(x = x), inherit.aes = FALSE,
           sides = "b", alpha = 0.6) +
  scale_color_manual(values = c("#009E73", "#E69F00", "#56B4E9")) +
  labs(title = "Combined", x = "days", y = "density") +
  theme_minimal()

grid.arrange(plnorm, pweib, pgamma_p, pcom, ncol = 2)

Interpretation:

  • AIC comparison — all three AICs fall within a range of 0.3 (Lognormal 90.9, Weibull 90.7, Gamma 90.6), well below the ΔAIC < 2 threshold for statistical distinguishability. No distribution can be ruled out on fit alone; lognormal is retained as the primary model based on the multiplicative mechanism for incubation periods.
  • Peak and location — all three curves peak near day 20–22 and cover the same observation range (17–27 days), confirming that the choice of family does not meaningfully affect the estimated central tendency.
  • Tail behaviour — the combined panel reveals the main difference: Weibull has a steeper left rise and slightly earlier peak than Lognormal and Gamma, which are nearly indistinguishable. All three assign non-trivial density beyond day 27, reflecting uncertainty in the right tail given n = 18.
  • CI ribbon width — the Weibull ribbon is visibly wider and more asymmetric than the other two, particularly on the upper bound at the peak. This reflects larger parameter uncertainty relative to the data scale for the Weibull parameterisation. The ribbon is computed by varying each parameter independently at its Wald CI bounds, ignoring off-diagonal Hessian terms (parameter covariance); the posterior predictive band in Section 3.6 corrects for this by propagating joint uncertainty.

2.4 Q-Q Plot (Lognormal MLE)

n     <- length(inc)
probs <- (seq_len(n) - 0.5) / n
qq_df <- data.frame(
  theoretical = qlnorm(probs, fit_lnorm$par[1], fit_lnorm$par[2]),
  observed    = sort(inc)
)
ggplot(qq_df, aes(x = theoretical, y = observed)) +
  geom_point(color = "#009E73", alpha = 0.7, size = 2) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  labs(title = "Q-Q Plot — Lognormal (MLE)",
       x = "Theoretical Quantiles", y = "Observed Quantiles") +
  theme_minimal()

Interpretation: Points close to the diagonal indicate good fit. Systematic curvature or S-shapes suggest the distribution family is inappropriate. Tied observations (8 cases at day 22) produce the horizontal cluster of points — expected for this dataset and not indicative of model misfit.

3. Bayesian Estimation via Stan (HMC/NUTS)

3.1 Model Specification

The lognormal family is selected based on AIC and the multiplicative-growth mechanism for incubation periods (Sartwell 1950). The Stan model encodes the same likelihood and weakly informative priors as v1.

stan_code <- "
data {
  int<lower=0> N;
  vector<lower=0>[N] y;
}
parameters {
  real mu;
  real<lower=0> sigma;
}
model {
  // Weakly informative priors -- to be updated after literature review.
  // mu    ~ Normal(3, 10): centred near log(20) days; wide sd lets likelihood dominate.
  // sigma ~ HalfNormal(1): lower=0 constraint enforces the half-normal truncation.
  mu    ~ normal(3, 10);
  sigma ~ normal(0, 1);
  y     ~ lognormal(mu, sigma);
}
"

Priors:

  • mu ~ Normal(3, 10): centred near log(20) ≈ 3; SD = 10 is weakly informative
  • sigma ~ HalfNormal(1): Stan enforces sigma > 0 via the <lower=0> constraint, making normal(0, 1) equivalent to a half-normal prior
  • These are weakly informative, and would be updated after a literature review of HFRS incubation periods.

Sampler: Stan uses the No-U-Turn Sampler (NUTS), an adaptive variant of Hamiltonian Monte Carlo (HMC). Unlike random-walk Metropolis-Hastings, NUTS uses the gradient of the log-posterior to propose efficient moves across the parameter space, avoiding the need to specify or tune a proposal distribution.

3.2 Sampling

stan_data <- list(N = length(inc), y = as.numeric(inc))

fit_stan <- stan(
  model_code = stan_code,
  data       = stan_data,
  chains     = 4,
  iter       = 10000,
  warmup     = 5000,
  seed       = 42,
  refresh    = 0
)

4 chains × 5,000 post-warmup iterations = 20,000 posterior samples.

3.3 HMC Diagnostics

check_hmc_diagnostics(fit_stan)
## 
## Divergences:
## 0 of 20000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 20000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
stan_sum  <- summary(fit_stan, pars = c("mu", "sigma"))$summary
diag_df   <- data.frame(
  Parameter = c("mu (meanlog)", "sigma (sdlog)"),
  Mean      = round(stan_sum[, "mean"], 4),
  R_hat     = round(stan_sum[, "Rhat"], 4),
  ESS       = round(stan_sum[, "n_eff"])
)
kable(diag_df, col.names = c("Parameter", "Posterior Mean", "R-hat", "ESS"),
      row.names = FALSE)
Parameter Posterior Mean R-hat ESS
mu (meanlog) 3.0628 1.0002 12824
sigma (sdlog) 0.1402 1.0002 11742

Interpretation:

  • R-hat — Gelman-Rubin statistic comparing between-chain to within-chain variance. Values < 1.01 confirm convergence across all 4 chains.
  • ESS (n_eff) — effective sample size after accounting for autocorrelation. HMC typically achieves much higher ESS than random-walk MH for the same number of iterations, because gradient-guided proposals decorrelate more rapidly. ESS > 400 is the conventional minimum.
  • Divergent transitions — a Stan-specific diagnostic. Any divergence indicates the sampler encountered regions where the posterior geometry is challenging (e.g., funnel-shaped). Zero divergences is the target.
  • E-BFMI — Energy Bayesian Fraction of Missing Information; values < 0.2 suggest the sampler struggles to explore the posterior energy landscape. Values near 1 are ideal.

3.4 Trace Plots

chains_arr <- extract(fit_stan, pars = c("mu", "sigma"),
                      permuted = FALSE, inc_warmup = FALSE)

df_trace <- do.call(rbind, lapply(seq_len(dim(chains_arr)[2]), function(k)
  data.frame(
    iter  = seq_len(dim(chains_arr)[1]),
    mu    = chains_arr[, k, "mu"],
    sigma = chains_arr[, k, "sigma"],
    chain = factor(k)
  )
))

p_tr_mu <- ggplot(df_trace, aes(x = iter, y = mu, color = chain)) +
  geom_line(alpha = 0.6, linewidth = 0.3) +
  labs(title = "Trace — mu (meanlog) [post-warmup]",
       x = "post-warmup iteration", y = "mu") +
  theme_minimal() + theme(legend.position = "none")

p_tr_sigma <- ggplot(df_trace, aes(x = iter, y = sigma, color = chain)) +
  geom_line(alpha = 0.6, linewidth = 0.3) +
  labs(title = "Trace — sigma (sdlog) [post-warmup]",
       x = "post-warmup iteration", y = "sigma") +
  theme_minimal() + theme(legend.position = "none")

grid.arrange(p_tr_mu, p_tr_sigma, ncol = 1)

Interpretation: Only post-warmup iterations (5,000 per chain × 4 chains) are shown; the warmup phase is excluded to avoid scale distortion from initial transients. All three convergence criteria are satisfied:

  • Stationarity — both parameters show no upward or downward drift throughout the 5,000 iterations: mu remains stable near 3.05 and sigma near 0.15, with no systematic trend at any point in the chain.
  • Mixing — the oscillations are large, rapid, and uniformly dense across the entire trace. There are no flat segments where the chain stalls at a fixed value, indicating that nearly every HMC proposal is accepted and the sampler moves freely through the posterior. This is the most visible advantage of HMC over random-walk MH: gradient-guided proposals produce large, effective jumps with low autocorrelation and high ESS per iteration.
  • Overlap — the four chains are fully interleaved with no separation, making individual chains indistinguishable by colour. This is consistent with R-hat < 1.01 for both parameters.

3.5 Posterior Distributions

post_stan <- as.data.frame(extract(fit_stan, pars = c("mu", "sigma")))

p_post_mu <- ggplot(post_stan, aes(x = mu)) +
  geom_density(fill = "#009E73", alpha = 0.5) +
  geom_vline(xintercept = fit_lnorm$par[1], linetype = "dashed", color = "gray30") +
  labs(title = "Posterior — mu (meanlog)",
       subtitle = "Dashed: MLE estimate",
       x = "mu", y = "density") +
  theme_minimal()

p_post_sigma <- ggplot(post_stan, aes(x = sigma)) +
  geom_density(fill = "#009E73", alpha = 0.5) +
  geom_vline(xintercept = fit_lnorm$par[2], linetype = "dashed", color = "gray30") +
  labs(title = "Posterior — sigma (sdlog)",
       subtitle = "Dashed: MLE estimate",
       x = "sigma", y = "density") +
  theme_minimal()

grid.arrange(p_post_mu, p_post_sigma, ncol = 2)

Interpretation:

  • mu (meanlog) — the posterior is narrow (range approximately 2.93–3.20) and nearly symmetric, with the MLE dashed line falling almost exactly at the mode. The close agreement between MLE and posterior mode confirms that the weakly informative prior exerts negligible influence; the shape is driven entirely by the likelihood. The tight width reflects strong data constraint on the location parameter given the concentrated observations.
  • sigma (sdlog) — the posterior is right-skewed with a sharp left boundary near zero, a mode around 0.13, and a long right tail extending to ~0.30. The MLE falls close to the mode. The skewness is expected: sigma is constrained to be positive, and with only 18 observations the upper tail remains uncertain. The right tail represents plausible scenarios where the true dispersion is somewhat larger than the data suggest — a genuine reflection of small-sample uncertainty rather than prior influence.

3.6 Posterior Predictive Check

set.seed(1)
n_samp   <- 500
samp_idx <- sample(nrow(post_stan), n_samp)
pp_mat   <- sapply(samp_idx, function(i)
  dlnorm(freq, post_stan$mu[i], post_stan$sigma[i]))

df_pp <- data.frame(
  x     = freq,
  mean  = rowMeans(pp_mat),
  lower = apply(pp_mat, 1, quantile, probs = 0.025),
  upper = apply(pp_mat, 1, quantile, probs = 0.975)
)

ggplot(df_pp, aes(x = x)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "#009E73", alpha = 0.25) +
  geom_line(aes(y = mean), color = "#009E73", linewidth = 1) +
  geom_rug(data = data.frame(x = inc), aes(x = x), inherit.aes = FALSE,
           sides = "b", alpha = 0.8) +
  labs(title = "Posterior Predictive — Lognormal (Stan/HMC)",
       subtitle = paste0("Mean density +/- 95% interval from ", n_samp,
                         " posterior samples"),
       x = "days", y = "density") +
  theme_minimal()

Interpretation: The posterior predictive distribution averages the lognormal density over all plausible (mu, sigma) combinations sampled from the posterior, propagating joint parameter uncertainty. Key observations from the plot:

  • Coverage — all 18 rug marks fall within the shaded band and close to the peak region (days 17–27), confirming no systematic misalignment between the model and the data.
  • Shape — the mean density curve peaks near day 21–22 and shows mild right-skew, consistent with a lognormal with small sigma (≈ 0.14). The distribution is nearly symmetric, reflecting the tight empirical clustering of this dataset.
  • Band width — the 95% interval spans roughly 0.09–0.19 at the peak, indicating moderate uncertainty about the exact height of the density given n = 18. The band widens on the right tail where fewer observations constrain the fit. This reflects parameter uncertainty, not model misfit.
  • Right tail — the non-zero density extending beyond day 30 represents the model’s uncertainty about the upper tail; no observations fall there, but the lognormal family assigns positive probability to longer incubation periods consistent with the fitted sigma.

4. Final Estimates

ip_median <- exp(post_stan$mu)
ip_mean   <- exp(post_stan$mu + post_stan$sigma^2 / 2)

est_df <- data.frame(
  Quantity         = c("Median incubation period (days)",
                       "Mean incubation period (days)",
                       "mu (meanlog)",
                       "sigma (sdlog)"),
  Posterior_Median = round(c(median(ip_median), median(ip_mean),
                              median(post_stan$mu), median(post_stan$sigma)), 2),
  CI_Lower         = round(c(quantile(ip_median, 0.025), quantile(ip_mean, 0.025),
                              quantile(post_stan$mu, 0.025),
                              quantile(post_stan$sigma, 0.025)), 2),
  CI_Upper         = round(c(quantile(ip_median, 0.975), quantile(ip_mean, 0.975),
                              quantile(post_stan$mu, 0.975),
                              quantile(post_stan$sigma, 0.975)), 2)
)
kable(est_df, col.names = c("Quantity", "Posterior Median", "95% CrI Lower",
                             "95% CrI Upper"),
      row.names = FALSE)
Quantity Posterior Median 95% CrI Lower 95% CrI Upper
Median incubation period (days) 21.39 20.02 22.86
Mean incubation period (days) 21.60 20.25 23.16
mu (meanlog) 3.06 3.00 3.13
sigma (sdlog) 0.14 0.10 0.20

Interpretation:

The posterior median incubation period is 21.39 days (95% CrI: 20.02–22.86), and the posterior mean is 21.60 days (95% CrI: 20.25–23.16). The narrow gap between median and mean reflects the small estimated sigma (0.14), indicating that the lognormal distribution is only mildly right-skewed for this dataset — consistent with the concentrated raw data (8 of 18 cases on day 22).

The sigma posterior (median 0.14, 95% CrI: 0.10–0.20) is notably small, meaning the fitted lognormal is nearly symmetric on the log scale. This is driven by the tight empirical spread of observations (range 17–27 days) rather than by the prior, as the weakly informative HalfNormal(1) prior places negligible constraint in this region.

Median and mean incubation periods are derived by transforming posterior samples back to the day scale: median = exp(mu), mean = exp(mu + sigma²/2). The 95% CrI directly expresses the probability that the true value lies within the stated range given the data and priors — a stronger inferential statement than the frequentist CI from MLE, which has no direct probabilistic interpretation for a fixed parameter.

5. Credits

Data & Literature

  • Kulagin, S.V., Fedorova, N.I., & Ketiladze, E.S. (1962). Laboratory outbreak of hemorrhagic fever with renal syndrome. Journal of Microbiology, Epidemiology and Immunology, 33(10), 121–126.

  • Sartwell, P.E. (1950). The distribution of incubation periods of infectious disease. American Journal of Hygiene, 51(3), 310–318. (Foundational reference for the lognormal model of incubation periods; source of the multiplicative-mechanism rationale.)

Software

  • R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
  • Stan Development Team. RStan: the R interface to Stan.
  • Wickham, H. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.

AI Assistance

This document was developed with the assistance of Claude (Anthropic). AI was used to support code writing, interpretation drafting, and statistical explanation. All analytical decisions, prior specifications, and interpretations were reviewed and approved by the author.