Model assessment
Residuals and reporting
A residual for ordered outcomes should preserve the direction of discrepancy without imposing numerical distances between categories. The probability-scale residual of Li and Shepherd (2012) compares the fitted probability below the observed category with the fitted probability above it. Its value is positive when the observed category is high relative to the fitted distribution, and negative when it is low.
probability_residual <- residuals(fit_br)
head(cbind(probability_scale=probability_residual,
residuals(fit_br,type="cumulative")))
#> probability_scale Low Medium
#> 1 -0.62111854 0.6211185 0.3344078
#> 2 0.04447362 -0.3788815 0.3344078
#> 3 0.66559216 -0.3788815 -0.6655922
#> 4 -0.74244554 0.7424455 0.4690656
#> 5 -0.21151111 -0.2575545 0.4690656
#> 6 0.53093444 -0.2575545 -0.5309344
by(seq_len(nrow(housing)),housing$Infl,function(i)
weighted.mean(probability_residual[i],housing$Freq[i],na.rm=TRUE))
#> housing$Infl: Low
#> [1] -0.0002456822
#> ------------------------------------------------------------
#> housing$Infl: Medium
#> [1] 0.0001079569
#> ------------------------------------------------------------
#> housing$Infl: High
#> [1] 0.0006339001The residuals lie between minus one and one. For a count-matrix response, the probability-scale residual is the mean of the individual residuals within that row. Cumulative residuals compare observed and fitted proportions through each threshold, which can help locate discrepancies in different parts of the outcome distribution. Zero-frequency rows have missing residuals. These are descriptive checks: ordinal probability-scale residuals are generally neither normal nor uniform, and cumulative residual columns are dependent.
For reporting and further analysis, optional broom
methods retain the distinction between common slopes and thresholds.
Exponentiation applies to the slopes and their intervals, yielding
common odds ratios; threshold estimates and coefficient standard errors
remain on the logit scale.
if(requireNamespace("broom",quietly=TRUE)) {
broom::tidy(fit_br,component="location",conf.int=TRUE,exponentiate=TRUE)
broom::glance(fit_br)
broom::augment(fit_br,newdata=profiles,se_fit=TRUE)
}
#> Infl Type Cont .fitted.Low .fitted.Medium .fitted.High .class
#> 1 Low Tower Low 0.3788815 0.2867107 0.3344078 Low
#> 4 Medium Tower Low 0.2575545 0.2733800 0.4690656 High
#> 7 High Tower Low 0.1445395 0.2108406 0.6446199 High
#> 10 Low Apartment Low 0.5189481 0.2598099 0.2212420 Low
#> .se.fit.Low .se.fit.Medium .se.fit.High
#> 1 0.02930516 0.01171912 0.02786562
#> 4 0.02326513 0.01273978 0.02956705
#> 7 0.01764028 0.01605866 0.03127093
#> 10 0.02783526 0.01312037 0.01996866The .fitted matrix column from augment()
contains one probability per category, and .class gives the
category with the largest fitted probability. Training-data augmentation
also includes probability-scale residuals. Supplied covariance matrices
can be used in summary(), broom::tidy(), and
average probability contrasts, while the prediction standard errors in
this example use expected Fisher information.
Prediction scores, calibration, and graphical slope comparisons
Good numerical convergence establishes that the estimating equations were solved, but the fitted probabilities also need to describe the outcome distribution adequately. Prediction scores summarize discrepancies, while calibration plots indicate where those discrepancies occur. Assessing predictions on observations excluded from fitting helps distinguish predictive performance from agreement with the training sample.
The housing data contain frequency rows rather than individual records. We randomly divide the residents represented by each row between training and evaluation samples, retaining the same predictor coding. The evaluation frequencies are used only after fitting the training model.
set.seed(61105)
training_frequency <- rbinom(nrow(housing),housing$Freq,.7)
evaluation_frequency <- housing$Freq-training_frequency
training_fit <- fastglm_ordinal(
Sat~Infl+Type+Cont,data=housing,weights=training_frequency,
firth=TRUE,tol=1e-10)
assessment <- fastglm_ordinal_assess(
training_fit,newdata=housing,y=housing$Sat,
weights=evaluation_frequency,bins=5)
assessment$scores
#> frequency log_loss rps rps_normalized
#> 1 545 1.038522 0.4246117 0.2123059The logarithmic loss penalizes assigning small probabilities to observed categories. The ranked probability score sums the squared discrepancies between predicted cumulative probabilities and the individual indicators of an outcome at or below each threshold. Both are proper prediction scores: their expected values are minimized by the true outcome probabilities; see Gneiting and Raftery (2007). Smaller values indicate better predictions for the same outcome scale and evaluation observations. The normalized ranked probability score divides by the number of thresholds, placing its possible values between zero and one. For grouped counts, the calculation averages individual scores; squaring the difference from a grouped observed proportion would measure a different quantity.
For each threshold, the calibration plot compares the mean predicted
probability with the observed proportion at or below that threshold. The
bins are equal-width intervals of predicted probability, chosen without
reference to the observed outcomes. Points near the diagonal indicate
agreement within those bins. Sparse bins can vary considerably, and the
frequency column in assessment$calibration records the
number of observations represented. Empty bins have unavailable observed
proportions. These plots are descriptive and can conceal discrepancies
within a bin; they do not establish that the full conditional
probability model is correct.
The proportional-odds model also assumes that a predictor has the same slope at each threshold. A descriptive comparison fits a separate binary model for an outcome above each threshold, retaining the original predictors and offsets. As illustrated in Harrell’s ordinal regression chapter, differences between estimates at different cutoffs can arise from sampling variation even when proportional odds holds. The auxiliary models use mean bias reduction to permit finite estimates in many sparse or separated samples.
assumption_check <- fastglm_ordinal_assess(fit_br,proportional.odds=TRUE)
head(assumption_check$proportional.odds[
,c("category","term","estimate","common.estimate","converged")])
#> category term estimate common.estimate converged
#> 1 Low InflMedium 0.5981626 0.5644198 TRUE
#> 2 Low InflHigh 1.2356223 1.2837886 TRUE
#> 3 Low TypeApartment -0.6124204 -0.5701275 TRUE
#> 4 Low TypeAtrium -0.1813150 -0.3649909 TRUE
#> 5 Low TypeTerrace -1.0758702 -1.0867312 TRUE
#> 6 Low ContHigh 0.4282444 0.3589079 TRUE
plot(assumption_check,type="proportional.odds")
Positive auxiliary slopes describe shifts toward higher outcomes, matching the ordinal model’s sign convention. Dashed horizontal lines mark the common slopes from the ordinal fit, and the vertical intervals are separate 95% normal intervals for the auxiliary slopes. Variation across thresholds can suggest that the common-slope restriction deserves further investigation. The estimates use overlapping observations and are correlated, so their displayed uncertainty cannot be interpreted as independent tests of equal slopes. The comparison is a graphical diagnostic, with failed threshold fits retained in the returned table; it does not provide a formal proportional-odds test.
Testing the common-slope restriction
We can assess the common-slope restriction with
fastglm_po_test(). The score measures how the log
likelihood would initially change if the slopes were allowed to differ
across thresholds. Its expected variation is calculated from the
category probabilities under the fitted proportional-odds model, with
adjustment for the estimated thresholds and common slopes. This
calculation uses expected Fisher information, obtained by averaging over
possible outcomes analytically. It requires neither simulated outcomes
nor resampling to compute a p-value.
po_test <- fastglm_po_test(fit_br)
po_test
#> Expected-information proportional-odds score test
#> Terms: Infl, Type, Cont
#> Null fit: maximum likelihood
#> Sample: 1681 observations represented by 72 rows
#>
#> Statistic df p-value
#> 8.729 6 0.1894
#>
#> Note: First-order chi-square approximation.By default, the test refits the proportional-odds null by maximum likelihood, including when the original model used mean bias reduction. The original fitted model and its predictions remain unchanged. The printed table identifies the null estimation method and uses a first-order chi-square approximation. Expected and observed information can give different score statistics in finite samples. The SAS/STAT description of the proportional-odds score test, on pages 4246 to 4248, allows either form; the implementation here uses expected information throughout.
The same test can be included with prediction scores by calling
fastglm_ordinal_assess(fit_br,po.test="score"). Printing
that assessment displays its scores and test in separate tables.
Printing its proportional.odds.test component displays the
test alone. For a report or subsequent calculation,
as.data.frame(po_test) and
broom::tidy(po_test) provide numeric columns with the full
computed precision; display rounding does not change those values.
For \(s\) predictor columns and \(K\) categories, the joint test has \(s(K-2)\) degrees of freedom. A small p-value gives evidence against the common-slope model in the specified direction of departure. A large p-value does not establish proportional odds, particularly when the sample contains little information about a category or predictor. The graphical comparison remains useful for understanding which threshold-specific associations differ and whether those differences are substantively important.
We can also select formula terms. The following test allows the slopes associated with housing type to differ across thresholds while retaining common slopes for perceived influence and contact. Selecting a factor includes all of its model-matrix columns.
type_test <- fastglm_po_test(fit_br,terms="Type")
type_test
#> Expected-information proportional-odds score test
#> Terms: Type
#> Null fit: maximum likelihood
#> Sample: 1681 observations represented by 72 rows
#>
#> Statistic df p-value
#> 5.908 3 0.1162
#>
#> Note: First-order chi-square approximation.
as.data.frame(type_test)
#> statistic df p.value null.estimation available
#> 1 5.907505 3 0.116198 ML TRUEThe joint p-value is 0.189, and the housing-type p-value is 0.116. Neither comparison gives strong evidence against its specified common-slope restriction in these data. This finding should be considered alongside the plotted associations and the prediction assessment, since limited power can also produce a large p-value.
This selected-term test still assumes that the other terms satisfy
proportional odds. It does not fit a model with unrestricted
threshold-specific effects for those terms. For a matrix fit,
terms selects named predictor columns. Within
fastglm_ordinal_assess(), the corresponding argument is
po.terms. Frequency weights, grouped category counts, and
location offsets follow the same conventions as fitting, without
expanding grouped observations.
When several terms are of interest, by.term=TRUE
provides both the joint test and a separate test for each selected term.
The calculations share the null fit and its information, which avoids
repeated fitting and passes through the predictor data.
po_terms <- fastglm_po_test(fit_br,by.term=TRUE)
po_terms
#> Expected-information proportional-odds score test
#> Null fit: maximum likelihood
#>
#> Hypothesis Statistic df p-value Available
#> Joint 8.729 6 0.1894 TRUE
#> Infl 1.100 2 0.5770 TRUE
#> Type 5.908 3 0.1162 TRUE
#> Cont 2.435 1 0.1186 TRUE
#>
#> Note: First-order chi-square approximation. Term-specific p-values are
#> unadjusted for multiple testing.
as.data.frame(po_terms)
#> hypothesis statistic df p.value null.estimation available
#> 1 Joint 8.729302 6 0.1893840 ML TRUE
#> 2 Infl 1.099725 2 0.5770291 ML TRUE
#> 3 Type 5.907505 3 0.1161980 ML TRUE
#> 4 Cont 2.435222 1 0.1186371 ML TRUEEach term-specific alternative relaxes only that term’s common-slope
restriction. The reported p-values are unadjusted for multiple testing;
the joint test provides a single assessment of all selected
restrictions. Individual results remain accessible through
po_terms$tests, including explanations for unavailable
comparisons.
The chi-square approximation requires independent observations and a
correctly specified, identifiable null model with positive category
probabilities. It can be inaccurate in small or sparse samples. If a
declared category is empty, the design is deficient, the ML null does
not converge, or the expected information cannot be used reliably,
available is false and the statistic and p-value are
missing; reason and detail explain the
failure. A finite bias-reduced fit therefore does not guarantee that
this ML-based test is available. The function does not silently remove
categories or switch tests.
When assessment outcomes are supplied through newdata
and y, the test refits its null to those outcomes using the
complete rows with positive frequency. It then assesses proportional
odds in that evaluation population. This differs from the prediction
scores above, which assess predictions from the original fitted
coefficients. Neither calculation establishes independence or correct
specification of omitted predictors.
Separation and a bias-reduced null
When a predictor group contains only the highest outcome, increasing its slope can keep improving the likelihood without reaching a finite maximum. This is a form of separation. Additional iterations cannot produce a finite ML estimate in such a sample. The test checks for this configuration in predictors with two observed values and reports it directly.
exposure <- matrix(c(0,1),ncol=1,dimnames=list(NULL,"exposed"))
sparse_counts <- rbind(c(12,18,20),c(0,0,4))
colnames(sparse_counts) <- c("low","middle","high")
sparse_fit <- fastglm_ordinal(exposure,sparse_counts,firth=TRUE)
fastglm_po_test(sparse_fit)
#> Expected-information proportional-odds score test
#> Terms: exposed
#> Null fit: maximum likelihood
#> Sample: 54 observations represented by 2 rows
#> Status: Unavailable
#> Reason: Predictor 'exposed' at 1 occurs only in endpoint category 'high'. The
#> ML null has no finite solution. A separately requested meanBR null may be
#> used; its test remains an approximation.
#>
#> Note: First-order chi-square approximation.An explicitly requested mean-bias-reduced null can provide finite estimates under separation. The test then evaluates the ordinary score at that null and removes the contribution associated with estimating its thresholds and common slopes. This removal is essential because the ordinary nuisance score is nonzero at a bias-reduced estimate.
fastglm_po_test(sparse_fit,null.estimation="meanBR")
#> Expected-information proportional-odds score test
#> Terms: exposed
#> Null fit: mean bias reduction
#> Sample: 54 observations represented by 2 rows
#>
#> Statistic df p-value
#> 0.001896 1 0.9653
#>
#> Note: Ordinary efficient score at a mean-bias-reduced null; first-order
#> chi-square approximation.This option has the same first-order chi-square limit under a regular
interior null, but it can be conservative in small samples. It is an
ordinary efficient-score calculation at a bias-reduced null, rather than
the larger-model adjusted-score statistic discussed by Kosmidis (2014, Section
12). A higher-order accuracy claim would require a separate
justification. The choice of null estimator should be specified as part
of the analysis; the function never switches estimators automatically.
Within fastglm_ordinal_assess(), use
po.estimation="meanBR" together with
po.test="score".
Bias reduction does not supply information about a predictor level
that is entirely absent. Empty declared outcome categories and deficient
designs therefore remain explicit unsupported cases for these tests.
Failed null fits retain their convergence details. The printed
explanation can be exported with
as.data.frame(result,details=TRUE).
