Modelling Norms for Speeded Tests and Count Data with the Conway-Maxwell-Poisson (CMP) Distribution

Wolfgang Lenhard & Alexandra Lenhard

2026-10-03

Introduction

In a speeded test, the raw score is the number of correctly processed units within a fixed time limit: symbols matched in a coding task, targets marked in a cancellation test such as the d2, words read correctly per minute in a fluency task. In many speeded tasks, items are deliberately kept relatively easy so that performance is substantially constrained by processing speed rather than by item difficulty.

Statistically, such scores are counts, \(y \in \{0, 1, 2, \dots\}\), and there are two standard count distributions. These however display severe limitations regarding speeded tests:

To model speeded tests, more flexibility is needed: In speeded tasks, underdispersion may arise when individual response rates are relatively regular, whereas heterogeneous response rates can generate overdispersion (see, e.g., Forthmann et al., 2024). Since the defining equidispersion property of the Poisson distribution is \[E(Y)=Var(Y)=\lambda \], a quick check is to compare the variance with the mean: at a mean of 40 raw points, a standard deviation clearly below 6.3 indicates under-dispersion, one clearly above it over-dispersion. Which of the two occurs depends on the test and the age group; it is an empirical question.

The Conway-Maxwell-Poisson (CMP) distribution relaxes this constraint. It generalises the Poisson distribution by a dispersion parameter \(\nu\) and thereby covers over-, equi- and under-dispersion in one distribution family. cnorm.cmp() builds a continuous norming model on this distribution: location and dispersion are both smooth functions of age, and norm scores follow from the percentile ranks of the fitted distribution. Additionally, many speeded tests are also bounded. They have a maximum number of items, e. g. the number of items in a test solved within a time constraint. For these scenarios, the distribution can be truncated at that maximum (max_score, see below). Thus, flexible count models for timed assessments are a tool for very prominent use cases in psychometrics.

The CMP Model

The CMP distribution (Conway & Maxwell, 1962; Shmueli et al., 2005; for regression see Sellers & Shmueli, 2010) introduces the \(\nu\) (‘nu’) parameter to allow dispersion to vary:

\[P(Y = y) = \frac{1}{Z(\mu, \nu)} \left(\frac{\mu^y}{y!}\right)^{\nu}, \qquad Z(\mu, \nu) = \sum_{j=0}^{\infty} \left(\frac{\mu^j}{j!}\right)^{\nu}, \qquad y = 0, 1, 2, \dots\]

with \(\mu > 0\) and \(\nu > 0\). For \(\nu = 1\), this is the Poisson distribution with mean \(\mu\). Values \(\nu > 1\) concentrate the distribution around its centre (under-dispersion), values \(\nu < 1\) spread it out (over-dispersion). \(\mu\) is approximately the mean, and the variance-to-mean ratio is approximately \(1/\nu\). Exact moments for a fitted model are available through predictMoments().

The figure shows three CMP distributions with the same \(\mu = 20\). Only \(\nu\) changes, and with it the spread.

When modelling over an explanatory variable like age, location and dispersion are modelled with polynomials. We use the logarithms of the parameters to stabilise estimation and to keep both parameters positive:

\[\log \mu(a) = \sum_{k=0}^{K_\mu} \beta_k a^k, \qquad \log \nu(a) = \sum_{k=0}^{K_\nu} \gamma_k a^k, \qquad a = \frac{\text{age} - \overline{\text{age}}}{\text{SD}(\text{age})}.\]

All coefficients are estimated simultaneously by maximum likelihood (L-BFGS-B with analytic gradients).

A First Look at the Data

For demonstration, we use the speed dataset. This is a synthetic dataset in the cNORM package, which resembles the highly speeded word-picture matching task from the ELFE-II reading comprehension test (Lenhard, Lenhard & Schneider, 2017). It has a maximum item number of k = 75 and a very strict time cutoff of 3 minutes. Each item is composed of a picture and four words. The child has to underline the correct alternative as quickly as possible. Within the test, errors happen rarely and the performance is mostly determined by decoding speed. The raw score (fluency) is the number of correctly solved items with the maximum attainable score of 75.

head(speed)
#>     age fluency
#> 1  8.77      45
#> 2  9.63      63
#> 3  9.92      15
#> 4  7.71       2
#> 5  8.62      35
#> 6 10.67      48
range(speed$fluency)
#> [1]  2 75

We first have a look at the dispersion index \(\text{Var}/\text{Mean}\), separated by age. If the model assumptions of a Poisson distribution holds, this index should stay roughly at 1 at any age (= equi-dispersion). Values below 1 represent under-dispersion, values above 1 over-dispersion.

by_year <- split(speed$fluency, floor(speed$age))
round(t(sapply(by_year, function(x)
  c(n = length(x), mean = mean(x), sd = sd(x), dispersion = var(x) / mean(x)))), 2)
#>      n  mean    sd dispersion
#> 7  288 25.05 10.96       4.80
#> 8  324 36.59 11.58       3.66
#> 9  369 46.82 13.53       3.91
#> 10 364 50.07 11.25       2.53

Fitting the Model

cnorm.cmp() fits the model. By default, it uses a cubic polynomial for the location (mu_degree = 3) and a quadratic polynomial for the dispersion (nu_degree = 2). We start by treating the scale as open-ended:

model.open <- cnorm.cmp(age = speed$age, score = speed$fluency,
  mu_degree = 3, nu_degree = 2, # default values
)
#> 'max_score' not set. Assuming positively unbounded (open-ended) range of score values. Please set 'max_score' if there is an upper limit.

The plot shows the fitted percentile curves together with the empirical percentiles (diamonds). The curves are continuity-corrected quantiles, i.e., the inverse of the mid-p percentile rank used for the norm scores below, which makes them directly comparable to the empirical values.

The mu_degree and nu_degree can be adjusted to apply a stronger smoothing (lower values) or to increase the closeness of the fit (higher values). There are two special settings of the dispersion worth knowing: nu_degree = 0 estimates a constant \(\nu\) over age, and nu_degree = NULL fixes it at the value given in nu, so that nu = 1 yields an ordinary Poisson regression. Keep the polynomial degrees as low as the data allow, since the dispersion is usually the less well-determined component. High values may result in overfitting. Internally, \(\nu\) is restricted to the range \(e^{-3}\) to \(e^{3}\) resulting in a range of \([0.05, 20]\). If its predictor reaches this limit at an observed age, the function issues a warning and the standard errors are not valid. In that case, reduce nu_degree.

Tests with a Ceiling: max_score

The speed test contains a finite number of targets, so scores above 75 are impossible. The open-ended CMP distribution is untruncated, however. It assigns a small probability to scores above the ceiling, and that probability is not available to the scores that can actually occur. As long as performance stays well below the maximum, the effect is negligible; in older age groups, where performance approaches the ceiling, it can bias the models. We can check this for the oldest age in the sample:

age_old <- max(speed$age)
tab <- normTable.cmp(model.open, ages = age_old, start = 0, end = 150)[[1]]
p_above <- 1 - tab$Pcum[tab$x == 75]

At age 11, the open-ended model assigns a probability of 2.2% to scores above 75. Where this probability is not negligible, percentile ranks near the top of the scale are not properly normalized. The CMP implementation in cNORM offers modelling truncated CMP distributions: If a maximum score \(M\) is supplied via max_score, the distribution is truncated at \(M\):

\[P(Y = y) = \frac{(\mu^y / y!)^{\nu}}{\sum_{j=0}^{M} (\mu^j / j!)^{\nu}}, \qquad y = 0, 1, \dots, M.\]

The norm score now stops at \(M\). Probabilities add up to one on the attainable scores, percentile ranks are normalized at the ceiling, and norm tables end at \(M\). Note that \(\mu\) and \(\nu\) now describe the truncated distribution: once the ceiling carries noticeable probability, \(\mu\) is no longer approximately the mean, and mean and variance should be obtained from predictMoments().

model.cmp <- cnorm.cmp(
  age = speed$age,
  score = speed$fluency,
  mu_degree = 3,
  nu_degree = 2,
  max_score = 75,
  plot = TRUE
)

When should max_score be used?

Checking the Model

summary() combines fit statistics, convergence information and parameter estimates. If age and score are supplied, it also reports R², RMSE and bias of the predicted against the manifest norm scores and a calibration table by age group:

summary(model.cmp, age = speed$age, score = speed$fluency)
#> Conway-Maxwell-Poisson Continuous Norming Model
#> -----------------------------------------------
#> Polynomial degrees:
#>   Location (log mu): 3 
#>   Dispersion (log nu): 2  
#>   Right-truncated at max_score = 75 
#> Number of observations: 1345 
#> Number of parameters: 7 
#> 
#> Model Fit:
#>   Log-likelihood: -5182.62 
#>   AIC: 10379.24 
#>   BIC: 10415.67 
#>   R-squared: 0.9637 
#>   RMSE: 1.9023 
#>   Bias: 0.0405 
#> 
#> Convergence:
#>   Converged: TRUE 
#>   Function evaluations: 29 
#>   Max |gradient|: 0.00634 
#>   Hessian condition number: 43.6 
#>   Message: Successful convergence 
#>   Optimizer message: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH 
#> 
#> Parameter Estimates:
#> Location (mu) parameters (log scale):
#>            Estimate Std. Error  z value  Pr(>|z|)
#> log(mu)_0  3.733959    0.01362 274.0972 0.000e+00
#> log(mu)_1  0.288299    0.02355  12.2424 1.845e-34
#> log(mu)_2 -0.096141    0.01158  -8.3011 1.032e-16
#> log(mu)_3 -0.002426    0.01276  -0.1902 8.492e-01
#> 
#> Dispersion (nu) parameters (log scale):
#>           Estimate Std. Error  z value  Pr(>|z|)
#> log(nu)_0 -1.33511    0.06838 -19.5242 6.832e-85
#> log(nu)_1  0.10658    0.04999   2.1321 3.300e-02
#> log(nu)_2 -0.05133    0.05406  -0.9494 3.424e-01
#> 
#> Calibration by age group (z-residuals; expected: mean 0, SD 1):
#>    age   n  mean_z  sd_z
#>   7.32 103 -0.0591 1.171
#>   7.58 103  0.0440 0.827
#>   7.90 104 -0.0443 0.989
#>   8.24 104  0.0886 1.050
#>   8.56  96 -0.0432 0.916
#>   8.86 111  0.0415 0.841
#>   9.19 101 -0.0919 1.007
#>   9.47 103  0.0216 1.036
#>   9.72 105  0.0654 1.190
#>   9.99 104  0.0691 1.116
#>  10.21 104 -0.2047 0.839
#>  10.48  99  0.0341 0.978
#>  10.79 108  0.0685 0.960

Some guidance on reading the output:

Automatic Model Selection via BIC

To balance developmental flexibility against parsimony, autoselect.cmp() evaluates a grid of polynomial degrees for \(\mu\) and \(\nu\) and selects the model with the lowest Bayesian Information Criterion (BIC). Candidates that converge and keep \(\nu\) inside its admissible range at all observed ages are preferred, and the warnings of the individual fits are recorded rather than discarded. Each fit starts from the solution of the next lower degree of \(\mu\), and the search runs in parallel across the degrees of \(\nu\).

# Grid: degrees 1 to 4 for mu, 0 to 2 for nu (0 = constant dispersion).
# Use min_nu = -1 to include the Poisson model (nu fixed at 1) as a candidate.
model.best <- autoselect.cmp(
  age = speed$age,
  score = speed$fluency,
  max_mu = 4,
  max_nu = 2,
  min_nu = 0,
  max_score = 75
)

# Evaluated models, sorted by BIC
head(model.best$selection$evaluated)

BIC is a guide, not a verdict. Always inspect the percentile plot and the calibration of the selected model.

Norm Scores

Norm scores can be predicted for individual cases or tabulated for specific ages.

Individual Predictions

Raw scores are discrete, so many persons share the same score. Assigning them the percentile rank \(P(Y \le x)\) would place all of them at the upper end of their tied block. cNORM therefore uses the mid-p percentile rank by default,

\[\text{PR}(x) = 100 \cdot \left[ P(Y < x) + \tfrac{1}{2}\, P(Y = x) \right],\]

and converts it into a norm score via the normal quantile function, \(z = \Phi^{-1}(\text{PR}/100)\), followed by the linear transformation to the chosen scale (T-scores by default, with \(M = 50\) and \(SD = 10\)). Norm scores are limited to \(\pm 3\) standard deviations unless a different range is given.

cases <- data.frame(
  age = c(7.5, 8.5, 9.5, 10.5),
  raw = c(18, 24, 30, 35)
)

cases$T_Score <- round(predict(model.cmp, age = cases$age, score = cases$raw), 1)
cases$PR <- round(predict(model.cmp, age = cases$age, score = cases$raw,
                          type = "percentile"), 1)
cases
#>    age raw T_Score   PR
#> 1  7.5  18    44.5 29.0
#> 2  8.5  24    39.8 15.3
#> 3  9.5  30    36.5  8.8
#> 4 10.5  35    36.9  9.6

The argument type selects the output: "norm" (default), "z" or "percentile". Percentile ranks are never limited by range.

Norm Tables

Norm tables are generated with normTable() (or normTable.cmp()). If a reliability coefficient is supplied, confidence intervals based on Kelley’s true score estimate (regression to the mean) are added:

tables <- normTable(
  c(8.0, 9.0),
  model.cmp,
  reliability = 0.88,
  CI = 0.90
)

head(tables[["8"]], 10)
#>    x           Px         Pcum  Percentile         z
#> 1  0 0.0001114256 0.0001114256 0.005571279 -3.000000
#> 2  1 0.0002356861 0.0003471117 0.022926862 -3.000000
#> 3  2 0.0004268340 0.0007739457 0.056052865 -3.000000
#> 4  3 0.0007059008 0.0014798464 0.112689604 -3.000000
#> 5  4 0.0010945733 0.0025744198 0.202713309 -2.873908
#> 6  5 0.0016145088 0.0041889285 0.338167414 -2.708277
#> 7  6 0.0022861312 0.0064750598 0.533199415 -2.553521
#> 8  7 0.0031272841 0.0096023439 0.803870183 -2.407154
#> 9  8 0.0041518785 0.0137542224 1.167828316 -2.267551
#> 10 9 0.0053686481 0.0191228705 1.643854647 -2.133581
#>        norm  lowerCI  upperCI lowerCI_PR upperCI_PR
#> 1  20.00000 18.25486 28.94514 0.07504378   1.762452
#> 2  20.00000 18.25486 28.94514 0.07504378   1.762452
#> 3  20.00000 18.25486 28.94514 0.07504378   1.762452
#> 4  20.00000 18.25486 28.94514 0.07504378   1.762452
#> 5  21.26092 19.36447 30.05475 0.10936255   2.304735
#> 6  22.91723 20.82202 31.51230 0.17625620   3.224552
#> 7  24.46479 22.18387 32.87416 0.27044787   4.339453
#> 8  25.92846 23.47191 34.16219 0.39912468   5.662178
#> 9  27.32449 24.70041 35.39069 0.57037873   7.201719
#> 10 28.66419 25.87935 36.56963 0.79312167   8.962997

The table contains:

For a model with max_score, the table ends at the maximum score.

Model-Implied Moments

predictMoments() evaluates mean, variance, skewness and (excess) kurtosis of the fitted score distribution at any age by exact summation of the probability mass function. For a model with max_score, the moments refer to the truncated distribution. We add the dispersion index to see how the spread relates to the Poisson benchmark:

moments <- predictMoments(model.cmp, age = seq(7, 11, by = 1), censor = FALSE)
moments$dispersion <- moments$variance / moments$mean
round(moments, 2)
#>   age  mean    sd variance skewness kurtosis dispersion
#> 1   7 19.50  9.77    95.36     0.59     0.34       4.89
#> 2   8 30.10 11.21   125.57     0.38     0.07       4.17
#> 3   9 41.79 12.08   145.85     0.16    -0.28       3.49
#> 4  10 49.64 11.96   143.11    -0.09    -0.48       2.88
#> 5  11 50.69 11.99   143.72    -0.14    -0.49       2.84

A dispersion index below 1 denotes under-dispersion, above 1 over-dispersion. Skewness and kurtosis show how the shape of the distribution changes across development. With censor = TRUE (the default), the distribution is instead censored at the empirical score range, which makes the moments comparable to those of the other model families.

Comparing Model Families

compare() places a CMP model next to alternative models, for instance the distribution-free Taylor polynomial model, and shows their percentile curves and fit statistics:

model.taylor <- cnorm(raw = speed$fluency, age = speed$age, plot = FALSE)

compare(
  model.taylor,
  model.cmp,
  age = speed$age,
  score = speed$fluency,
  title = "Taylor Polynomial vs. Conway-Maxwell-Poisson"
)
#> Retrieving norm scores, please stand by ...
#> 
#> Model Comparison Summary:
#> ------------------------
#>  Metric     Model1     Model2 Difference
#>      R2     0.9603     0.9569    -0.0034
#>    Bias     0.0145     0.0068    -0.0078
#>    RMSE     1.9971     2.0709     0.0738
#>     MAD     1.5342     1.5894     0.0552
#>     AIC  5767.0901 10379.2378  4612.1477
#>     BIC -5276.2840 10415.6669 15691.9508
#> 
#> Note: Difference = Model2 - Model1
#>       Fit indices are based on manifest and fitted norm scores (weighted if weights provided).
#>       Scale metrics are T scores (scaleSD = 10)
#>       AIC and BIC should only be used when comparing models of the same type.

Keep in mind that likelihood-based criteria (AIC, BIC) are comparable only between models of the same kind of distribution, such as the open-ended and the truncated CMP model above. They cannot be compared between discrete families (CMP, beta-binomial) and continuous ones (SHASH), nor with the distribution-free model, which has no likelihood. Across families, judge the models by the fit of the norm scores (R², RMSE, bias), by the visual inspection of the percentiles and by the calibration.

Choosing a Model Family

cNORM offers one semi-parametric, distribution-free approach and three parametric families. Each is tailored to a particular kind of test data.

Distribution-free (Taylor) Beta-binomial CMP SHASH
Score scale Continuous or discrete Integer, \(0 \le y \le n\) Integer, \(y \ge 0\) Continuous, \(y \in \mathbb{R}\)
Typical tests General ability and achievement Power and accuracy tests Speeded tests, rate tasks Continuous measures, response times
Upper bound Post-hoc clipping to \([minRaw, maxRaw]\) Natural bound (\(n\) items) Open-ended, or truncated via max_score None
Distributional shape Any (empirical) Binomial with beta mixing Over-, equi- or under-dispersed Skewed, light or heavy tails
Main strength Maximum flexibility Respects fixed item numbers Models the dispersion directly Handles negative and uncentred scores
Function cnorm() cnorm.betabinomial() cnorm.cmp() cnorm.shash()

Distribution-free models (cnorm(); Lenhard et al., 2018) are the choice when flexibility matters most: the data show irregular features that no standard family captures, you prefer not to assume a distribution, or the relationship between person location and age changes irregularly. Because bivariate polynomials can invert in extreme regions, always inspect checkConsistency() and plotPercentiles(). Parametric models cannot produce crossing percentile curves by construction.

Beta-binomial models (cnorm.betabinomial()) suit unspeeded power and accuracy tests in which the raw score is the number of correct responses among \(n\) items. They respect the fixed ceiling and handle floor and ceiling effects naturally. Sum scores of Rasch-type tests are approximately beta-binomial, which makes the model a good working description rather than an exact one.

CMP models (cnorm.cmp()) suit speeded tests whose raw score counts the units processed within a time limit, such as cancellation tasks, coding, symbol search or reading fluency. They model the dispersion as an age-dependent parameter instead of assuming it, covering under- and over-dispersion alike. If the test has a hard maximum score, set max_score. One limitation deserves mention: for large \(\mu\), the skewness of a CMP distribution is approximately \(1/\sqrt{\mu\nu}\), i.e., the shape is nearly determined by mean and variance. If skewness and kurtosis vary independently across age, SHASH or the distribution-free approach are more flexible.

SHASH models (cnorm.shash()) suit continuous raw scores: decimals, response times, physical measurements, or difference scores that can be negative. Four smooth age trajectories (location, scale, skewness and tail weight) describe the whole conditional distribution, which is parsimonious and flexible at the same time, and lets asymmetry and tail weight change across development.

References