--- title: "Modelling Norms for Speeded Tests and Count Data with the Conway-Maxwell-Poisson (CMP) Distribution" author: "Wolfgang Lenhard & Alexandra Lenhard" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Modelling Norms for Speeded Tests and Count Data with the Conway-Maxwell-Poisson (CMP) Distribution} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( tidy.opts = list(width.cutoff = 60), collapse = TRUE, comment = "#>" ) options(width = 60) library(cNORM) ``` # 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: - The **Poisson distribution** forces the variance to equal the mean (*equi-dispersion*). - The **negative binomial distribution** allows only a variance larger than the mean (*over-dispersion*). 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()`. ```{r pmf-figure, fig.height = 3.6, fig.width = 7, echo = FALSE} x <- 0:45 op <- par(mar = c(4, 4, 1, 1)) plot(x, dcmp(x, mu = 20, nu = 3), type = "l", lwd = 2, col = "#1b9e77", ylim = c(0, 0.17), xlab = "Raw score", ylab = "P(Y = y)", bty = "l") lines(x, dcmp(x, mu = 20, nu = 1), lwd = 2, col = "#7570b3") lines(x, dcmp(x, mu = 20, nu = 0.5), lwd = 2, col = "#d95f02") legend("topright", bty = "n", lwd = 2, col = c("#1b9e77", "#7570b3", "#d95f02"), legend = c(expression(nu == 3 ~ "(under-dispersed)"), expression(nu == 1 ~ "(Poisson)"), expression(nu == 0.5 ~ "(over-dispersed)"))) par(op) ``` 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. ```{r data} head(speed) range(speed$fluency) ``` 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. ```{r 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) ``` # 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: ```{r fig1, fig.height = 4, fig.width = 7} model.open <- cnorm.cmp(age = speed$age, score = speed$fluency, mu_degree = 3, nu_degree = 2, # default values ) ``` 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: ```{r ceiling-mass} 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 `r round(age_old, 1)`, the open-ended model assigns a probability of `r sprintf("%.1f%%", 100 * p_above)` 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()`. ```{r fig-ceiling, fig.height = 4, fig.width = 7} 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? - **Use it** whenever the scale has a hard upper limit that persons can reach or approach, typically the number of items or targets on the test sheet. Supply the *theoretical* maximum, not the highest observed score. Scores above `max_score` raise an error. - **Omit it** if there is no ceiling, or if it is practically out of reach within the time limit. The two models are then practically indistinguishable, and the open-ended one is simpler. - **Be critical** if a large share of the sample obtains the maximum score. For large $\mu$, the shape of a CMP distribution is nearly determined by its mean and variance, so strong ceiling effects may call for a closer look at the diagnostics below, or for a more flexible approach. # 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: ```{r} summary(model.cmp, age = speed$age, score = speed$fluency) ``` Some guidance on reading the output: - **Coefficients** are on the log scale, and $\log \mu_0$ and $\log \nu_0$ are the intercepts at the *mean age*. A positive $\log \nu_0$ means $\nu > 1$ and thus under-dispersion at the mean age; the higher-order coefficients describe how location and dispersion change with age. To see the dispersion at specific ages, use `predictMoments()` below. - **Convergence** should be successful with a small maximum gradient. A warning that $\log \nu$ reached its admissible range signals an over-parameterized dispersion model. - **Calibration**: If the model fits, the z-residuals have a mean close to 0 and a standard deviation close to 1 in every age group. The standard deviation can fall slightly below 1 for small counts, because the mid-p transform of a discrete variable has slightly less than unit variance. # 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$. ```{r, eval = FALSE} # 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. ```{r} 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 ``` 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: ```{r} tables <- normTable( c(8.0, 9.0), model.cmp, reliability = 0.88, CI = 0.90 ) head(tables[["8"]], 10) ``` The table contains: - **x**: raw score, - **Px**: point probability $P(Y = x)$, - **Pcum**: cumulative probability $P(Y \le x)$, - **Percentile**: mid-p percentile rank, - **z** and **norm**: z-score and norm score, limited to $\pm 3$ SD like in `predict()`, - **lowerCI** and **upperCI** (with **lowerCI_PR** and **upperCI_PR**): limits of the true-score confidence interval. 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: ```{r} moments <- predictMoments(model.cmp, age = seq(7, 11, by = 1), censor = FALSE) moments$dispersion <- moments$variance / moments$mean round(moments, 2) ``` 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: ```{r fig2, fig.height = 4, fig.width = 7} 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" ) ``` 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 * Conway, R. W., & Maxwell, W. L. (1962). A queuing model with state dependent service rates. *Journal of Industrial Engineering*, 12, 132-136. * Dunn, P. K., & Smyth, G. K. (1996). Randomized quantile residuals. *Journal of Computational and Graphical Statistics*, 5(3), 236-244. * Forthmann, B., Lenhard, W., Lenhard, A., & Förster, N. (2024). Flexible item response modeling for timed reading comprehension assessment. *The Journal of Experimental Education*, 1-17. https://doi.org/10.1080/00220973.2024.2367162 * Huang, A. (2017). Mean-parametrized Conway-Maxwell-Poisson regression models for dispersed counts. *Statistical Modelling*, 17(6), 359-380. * Lenhard, A., Lenhard, W., Suggate, S., & Segerer, R. (2018). A continuous solution to the norming problem. *Assessment*, 25(1), 112-125. https://doi.org/10.1177/1073191116656437 * Sellers, K. F., & Shmueli, G. (2010). A flexible regression model for count data. *The Annals of Applied Statistics*, 4(2), 943-961. * Shmueli, G., Minka, T. P., Kadane, J. B., Borle, S., & Boatwright, P. (2005). A useful distribution for fitting discrete data: revival of the Conway-Maxwell-Poisson distribution. *Journal of the Royal Statistical Society C*, 54(1), 127-142.