Overview

The package fits and predicts the full joint distribution of the 3,125 EQ-5D-5L health states as a function of age and sex. The model has three independent components, each estimated separately:

Per-dimension marginal models -- probabilities $P(\text{level}=l \mid \text{age}, \text{sex})$ for each dimension $k$ and level $l \in \{1,2,3,4,5\}$.
Gaussian copula -- a 5x5 correlation matrix R linking the dimensions, possibly stratified by age band x sex.
Mixture P(11111) correction (optional) -- a separate logistic model for the probability of full health.

Predictions evaluate the implied joint distribution by Sobol QMC integration of the multivariate-normal rectangle for each of the 3,125 states, then derive summary statistics (mean utility, SD, P(11111), quantiles, top-state rankings, per-dimension distributions) under any EQ-5D-5L value set carried by eq5dsuite.

1. Marginal models

For dimension k (one of MO, SC, UA, PD, AD), an ordinal regression predicts the probability of each severity level given age and sex. Age enters as a natural cubic spline with inner knots at 30, 45, 60, 75 years and boundary knots at the configured range (default [18, 100]).

Two model families are tried, with and without an age-by-sex interaction:

Proportional-odds (cumulative logit, MASS::polr): \[ \log \frac{P(L \le l)}{P(L > l)} = \alpha_l - \big( f(\text{age}) + \beta \cdot \text{female} + g(\text{age}) \cdot \text{female} \big) \] where $f, g$ are linear combinations of NS basis functions. Assumes the covariate effect is the same across all severity thresholds (parallel-regression assumption).
Multinomial logit (nnet::multinom): \[ \log \frac{P(L=l)}{P(L=1)} = \alpha_l + f_l(\text{age}) + \beta_l \cdot \text{female} + g_l(\text{age}) \cdot \text{female} \] Relaxes the parallel-regression assumption -- each level gets its own covariate effects.

For each dimension, the four candidates (polr x with/without interaction; multinom x with/without interaction) are fit and the best is selected by AIC.

Package code:

fit_marginal_models(data, dims = c("MO","SC","UA","PD","AD"),
                    weights = NULL,
                    boundary_knots = c(18, 100))
# returns list of `eqnorms_marginal` objects, one per dim, with $model,
# $model_name, $aic, $candidate_aics.

predict_marginal(object, age, female)
# returns an n x 5 matrix of level probabilities.

2. Gaussian copula

Given the marginals, the dependence between dimensions is captured by a Gaussian copula. Conceptually, each respondent has 5 latent continuous "health" scores $Z = (Z_1, Z_2, Z_3, Z_4, Z_5) \sim N(\mathbf{0}, R)$ where $R$ is a 5x5 correlation matrix; the observed level for dimension $k$ is the bin of $Z_k$ under thresholds $\tau_{k,1} < \tau_{k,2} < \tau_{k,3} < \tau_{k,4}$ chosen so the marginals match the per-dimension model:

\[ \tau_{k,j}(\text{age}, \text{sex}) = \Phi^{-1}\big( P(L_k \le j \mid \text{age}, \text{sex}) \big). \]

The probability of any health state $(s_1, ..., s_5)$ is then the 5-dimensional MVN rectangle probability:

\[ P(\text{state} = (s_1, ..., s_5) \mid \text{age}, \text{sex}) = \int_{a_1}^{b_1} \cdots \int_{a_5}^{b_5} \phi_5(\mathbf{z}; R)\, d\mathbf{z} \]

with $a_k = \tau_{k, s_k - 1}$ and $b_k = \tau_{k, s_k}$, where $\tau_{k,0} = -\infty$ and $\tau_{k,5} = +\infty$.

2.1 Estimating R (polychoric correlation)

The latent correlations $R_{ij}$ are estimated by polychoric correlation -- the maximum-likelihood (or two-step Olsson) estimator for bivariate normal correlation when only the bin-discretized levels are observed. Concretely, for each pair (i, j) we build the weighted 5x5 cross-tabulation:

\[ n_{lm} = \sum_{r : \text{dim}_i(r) = l,\, \text{dim}_j(r) = m} w_r \]

and feed it to polycor::polychor(), which finds the bivariate-normal threshold model $(\tau_i, \tau_j, \rho_{ij})$ that best reproduces those cell counts. Two methods are supported:

twoStep (default) -- thresholds from marginals, $\rho_{ij}$ by 1D optimisation of the bivariate-normal likelihood given those thresholds. Fast.
ML -- joint maximum-likelihood over thresholds and correlation. ~30x slower; gives almost identical $\rho_{ij}$ on sample sizes > 500.

The full pairwise matrix is assembled (10 unique off-diagonals for 5 dims) and clamped to positive-definiteness via eigenvalue thresholding if necessary.

Why this matters: Earlier releases used a "probit-midpoint" Pearson approximation in which every respondent at level $l$ for dimension $k$ was assigned the same latent score $\Phi^{-1}((F_{k,l-1} + F_{k,l})/2)$. Because within-level variation collapsed to zero, that estimator systematically attenuated correlations. Empirically, on Norwegian 2025 the mean off-diagonal $\overline{|R_{ij}|}$ was 0.33 under probit-midpoint vs 0.61 under polychoric, leading to predicted SDs ~0.027 below observed (see SD_diagnostic_report.html). Switching to polychoric closes the SD gap to ~0.001. Set options(eqnorms.copula_method = "midpoint") to recover the legacy behaviour.

Package code:

# Internal: weighted polychoric R for 5 EQ-5D dims
.polychor_R(data, weights, method = "twoStep")

# Stratified variants (independence, pooled, by_sex, by_age, by_age_sex)
fit_copula(data, marginals, variant = "by_sex", weights, method)
# stores $R as a matrix (pooled/independence) or a named list of matrices.

2.2 Stratification variants

The 5x5 R can be estimated globally, or stratified. Each variant adds parameters; users pick by trading off bias vs variance:

Variant	R varies by	unique R	params
independence	--	1 (identity)	0
pooled	--	1	10
by_sex	Sex	2	20
by_age	Age band (4)	4	40
by_age_sex	Age band x sex	8	80
continuous	Polynomial in age, per sex	smooth	~60

3. Joint distribution and Sobol QMC integration

Given the marginals $P$ (a 5x5 matrix of level probabilities at (age, sex)) and the copula correlation $R$ at that stratum, the package evaluates all 3,125 rectangle probabilities by Sobol Quasi-Monte Carlo integration of $N(0, R)$.

The procedure is:

Generate $N$ Sobol points $u_1, ..., u_N \in [0,1]^5$ (deterministic, low-discrepancy).
Convert to standard normals via the inverse CDF: $w_i = \Phi^{-1}(u_i)$.
Apply Cholesky: $z_i = w_i \cdot U$ where $R = U^T U$. Each row $z_i$ is one draw from $N(0, R)$.
Bin each draw into one of the 3,125 cells using the thresholds $\tau_{k,j}$, and count.
Probability of state $s = (s_1, ..., s_5)$ is $\text{count}(s)/N$, optionally smoothed.

One Sobol draw is reused across all 3,125 cells (the cells partition $\mathbb{R}^5$, so each draw lands in exactly one). At $N = 262144$ the maximum per-state error vs the deterministic mvtnorm integration is <5x10^-4, while a single 3,125-state rectangle takes ~25 ms vs ~4.5 s for mvtnorm.

Optional Laplace smoothing:

\[ \widehat{p}_s = \frac{c_s + \alpha}{N + \alpha K} \]

with $K = 3125$. Set $\alpha > 0$ when downstream code needs strictly positive probabilities (e.g. for $\log$ or division); $\alpha = 0.01$ shifts displayed statistics by far less than the reporting precision.

Package code:

# Low-level kernel (one (age, sex) stratum):
compute_joint_qmc_kernel(P, R, N = 262144L, alpha = 0)

# High-level (handles age sweep + mixture + value set + many statistics):
predict(fit, ages = 18:100, sexes = c(0L, 1L),
        value_set = "NO",
        what = c("mean_utility","sd_utility","p_full","quantiles"))

4. Mixture P(11111) correction (deprecated default)

Earlier package versions (≤ 0.6.2) fit an optional logistic model for $P(\text{state}=11111 \mid \text{age}, \text{sex})$ and rescaled the copula joint to pin the corner. This was originally introduced to compensate for the probit-midpoint copula method (used before v0.5.0), which systematically attenuated correlations and under-predicted P(11111) by ~0.05.

Since v0.5.0 the package uses polychoric correlations, which fit the corner mass per stratum directly. Diagnostic comparison on the Norway 2025 norms (8 strata, weighted) shows the mixture is now slightly less accurate than the pure copula on every absolute error metric: |P(11111) err| pure 0.014 vs mixture 0.018, |mean utility err| 0.012 vs 0.013, |SD err| 0.008 vs 0.009. The mechanism: the mixture's GLM smooths over age, while the polychoric R adapts to each (age band, sex) stratum.

From v0.6.3 the default is mixture = FALSE. The code remains available for backward-compatible re-runs (fit_norms(..., mixture = TRUE)) but is no longer drawn in the dashboard.

5. From the joint distribution to reported statistics

5.1 Mean and SD of utility

\[ E[U \mid \text{age, sex}] = \sum_{s=1}^{3125} p_s \cdot u(s, \text{value set}) \] \[ \mathrm{SD}[U] = \sqrt{ \sum_s p_s (u(s) - E[U])^2 } \]

For additive value sets (Norway, UK, most published 5L sets), the pure-copula $E[U]$ depends only on the marginals -- it is invariant under the choice of copula R. SD on the other hand depends on the joint, so it is sensitive to R.

Note: if a fit was built with mixture = TRUE (no longer the default), the joint loses exact marginal consistency and $E[U]$ acquires a copula-dependent term $E[U]_{\text{mix}} = p_{\text{full}} + \frac{1 - p_{\text{full}}}{1 - p_{11111}^{\text{copula}}} (E[U]_{\text{copula}} - p_{11111}^{\text{copula}})$. The default pure-copula path avoids this side-effect.

5.2 P(full health)

$P(s = 11111)$ is the 5D corner probability. Higher correlations make co-occurrence at level 1 more likely, so this metric is sensitive to R.

5.3 Quantiles

The package computes weighted quantiles of utility over the discrete distribution $\{(u(s), p_s) : s = 1, ..., 3125\}$ -- sort utilities, accumulate probabilities, find crossing points by linear interpolation.

5.4 Top-N states

Sort the 3,125 cells by probability, take the top N. Observed proportions in the matching stratum (broad age band x sex) are joined in for comparison.

5.5 Per-dimension level distributions

Marginalising the joint:

\[ P(\text{dim}_k = l \mid \text{age, sex}) = \sum_{s : s_k = l} p_s \]

which should recover the marginal-model prediction up to QMC sampling noise -- the package exposes "consistency checks" that flag the difference.

6. Value sets and the dashboard

The fit object stores predictions and observed summaries scored under a single value set (recorded as fit$value_set). The dashboard's Value set selector can override this at any time; when the selected value set differs from the fit's, the dashboard rescores fit$data$EQ5D_index and recomputes the observed summaries on the fly so the observed line in every plot uses the same utility units as the predicted curves. Predictions also re-derive under the active value set. Caching is per (fit, value_set) so switching is fast after the first hit.

The active value set is displayed under the dashboard's "Value set" selector ("Active value set: ...") so it is never ambiguous.

7. Performance: precomputed dashboard curves

The Fitted-norms tab needs predicted curves (mean utility, SD, P(11111) per age x sex x copula variant; quantile bands for the active variant; marginal level distributions per dim and sex). These are deterministic given (fit, value_set), so the package precomputes them once and stores the result on the fit as fit$dashboard_curves (about 40 KB). Plot helpers consume the cache when available and fall back to live predict() calls only when the user picks a value set or age range outside the cached one. Effect: opening a pre-fit dashboard is ~5 s of plotly rendering instead of ~30 s of plotly + prediction.

8. References and further reading

Olsson U. (1979). Maximum likelihood estimation of the polychoric correlation coefficient. Psychometrika 44, 443-460.
Genz A., Bretz F. (2002). Comparison of methods for the computation of multivariate t probabilities. JCGS 11, 950-971.
Sklar A. (1959). Fonctions de repartition a n dimensions et leurs marges. Publ. Inst. Statist. Univ. Paris 8, 229-231.
Garratt A.M. et al. (2024). Norwegian EQ-5D-5L value set. Quality of Life Research. (Norwegian "NO" value set used by default throughout the demonstration data.)

The full source for every formula here is in the package R/ folder; in particular, see R/copula.R (correlation estimation), R/joint_qmc.R (Sobol kernel), R/marginals.R (per-dim models), R/mixture.R (P(11111) correction), and R/predict_norms.R (the dispatcher used by the dashboard).

Pre-stored datasets

Or upload your own CSV

Active dataset

Column mapping

Data preview

Upload population targets

Balance check

Variable matching

Copula R (Male):

Copula R (Female):

Model description:

Marginal Level Distributions by Age (Observed vs Predicted)

Mean Utility by Age

P(Full Health) by Copula Model

SD of Utility by Copula Model

Mean Utility by Copula Model

Model Comparison: Predictive Accuracy by Copula Variant

Dimension-Level Distributions

Utility Distribution

Top 20 Predicted States -- Males

Top 20 Predicted States -- Females

Model Validation -- Males

Consistency Checks:

Model Validation -- Females

Consistency Checks:

Model Descriptions

R Code: Reproduce Predictions

Fit overview

Marginal models (per dimension, AIC search)

Marginal-model coefficients (best model per dim)

Copula variants stored on this fit

Correlation matrices

Mixture (P(full health)) model

Overview

1. Marginal models

2. Gaussian copula

2.1 Estimating R (polychoric correlation)

2.2 Stratification variants

3. Joint distribution and Sobol QMC integration

4. Mixture P(11111) correction (deprecated default)

5. From the joint distribution to reported statistics

5.1 Mean and SD of utility

5.2 P(full health)

5.3 Quantiles

5.4 Top-N states

5.5 Per-dimension level distributions

6. Value sets and the dashboard

7. Performance: precomputed dashboard curves

8. References and further reading

1. Package tarball

2. Standalone analysis script

3. Fit model file

4. Predicted proportions for all 3,125 health states