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. EQ-VAS norms

The VAS is modelled through the health states. A cumulative logit on the observed VAS values gives the VAS for each state, with one step per level of each dimension (an indicator that the state is at that level or worse), plus age and sex for what they add beyond the dimensions:

\[ \operatorname{logit} P(\text{VAS} \le v_j \mid s, a, x) = \theta_j - \Big(\sum_{k,l} \gamma_{kl}\, 1[s_k \ge l] + f(a) + \beta x\Big) \]

The norm at age $a$ and sex $x$ is the mixture over the 3,125 states, $P(\text{VAS} = v \mid a, x) = \sum_s P(s \mid a, x)\, P(\text{VAS} = v \mid s, a, x)$, and its mean, SD and percentiles are computed from that distribution. $f(a)$ is linear below 2,000 respondents with a VAS and a natural spline from there. On three surveys this matched or beat modelling the VAS directly at every sample size tested (report/vas.html).

9. Norms by age alone (sexes pooled)

A norm without sex is the average of the distributions for men and women, weighted by the share of women at each age, $\pi(a)$:

\[ P(s \mid a) = (1 - \pi(a))\, P(s \mid a, \text{male}) + \pi(a)\, P(s \mid a, \text{female}) \]

The averaging is done on the distribution of states (or of VAS values), and every statistic is then computed from the averaged distribution. Statistics such as the SD and percentiles are never averaged directly. $\pi(a)$ is a logistic spline in age fitted with the same weights as the model, so in a weighted fit it is the share of women in the population. A share from a register can be used instead. In the population this identity is exact. On held-out data in three surveys it was as accurate as fitting a model without sex, and one fit then gives both the sex-specific and the combined norms.

The observed points shown with a pooled norm are mixed in the same proportion. The weight of each respondent is scaled by $\pi(a) / s(a)$ for women and $(1 - \pi(a)) / (1 - s(a))$ for men, where $s(a)$ is the fitted share of women at that age in the data the model was fitted to (weighted or not, as the fit). With the fitted share the weights are unchanged; with no women only the men are left.

10. Inspecting the copula

The Inspect tab draws the model behind the norms for two or three dimensions. Each dimension is read from a latent continuous position through four thresholds, $\tau_k = G^{-1}(F(k))$, where $F$ is the cumulative level probability at the age and sex shown and $G$ the normal or t distribution of the family of the copula. The copula is the shape of the joint cloud of latent positions, and a health state is a box whose probability is the mass of the copula in the box. With two dimensions the tab shows the latent density with the threshold planes; with three, simulated people (the same draws at every age, so they move as the age changes). The draws are a Sobol sequence by default: the points are spread evenly, so the share of people in each box is close to its probability (at 2,000 people the largest box error is about a third of that from random draws); random draws are an option. The box probabilities beside the figure are computed exactly (bivariate and trivariate normal or t rectangle probabilities), not from the drawing.

By default men and women are shown side by side in one figure: the same age, the same view (turning or zooming one scene turns the other) and, with three dimensions, the same simulated draws, each placed by the correlation and thresholds of its own sex. Any difference between the two scenes is therefore a difference between the sexes in the model.

Pooled, for both sexes together, men and women have different thresholds, so their latent pictures do not share a box grid. The combined picture uses the thresholds of the combined margins, and maps the latent of each sex onto it with a monotone spline through the threshold points on the probability scale. The level bands of each sex map exactly onto the combined bands, so every box has exactly the share-weighted probability, and the drawn density is the share-weighted sum of the two mapped densities. When a fit has a mixture or correction layer, the figures show the copula beneath it.

11. 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

Model specification

Or upload your own CSV

Active dataset

Column mapping

Data preview

Upload population targets

Balance check

Variable matching

Model description:

Share of Women by Age

Marginal Level Distributions by Age (Observed vs Predicted)

Mean Utility by Age

P(Full Health) by Age

SD of Utility by Age

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:

Top 20 Predicted States -- Both sexes pooled

Norms by age band, sexes pooled

Model Descriptions

R Code: Reproduce Predictions

Mean EQ-VAS by Age

EQ-VAS Distribution

VAS norms by age band

Highlight

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

Export the fitted norms

Load a distribution

Loaded distribution

Covariates recorded

Level probabilities

Correlations

Text rendering

1. Choose a distribution

2. Settings

What is being compared

Headline comparison

Per-dimension tests

Whole distribution

Per-dimension tests

Whole distribution

Probability of any problem

Level distributions

Level distributions

Distribution of index values

Health states with the largest differences

Health states with the largest differences

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. EQ-VAS norms

9. Norms by age alone (sexes pooled)

10. Inspecting the copula

11. References and further reading