Skip to content

Cross-validation & model selection

Choosing the number of latent classes is an important modelling decision for any finite mixture. Information criteria (BIC, CAIC, adjusted BIC) are useful, and LCL reports all three. But these criteria can disagree, and can prove unstable over repeated sampling (see Preacher and Merkle [2012, Psychological Methods]). Out-of-sample log-likelihood under a panel-respecting split is preferable when you can spare the time (and shoulder the GPU expense).

cv_optimal_classes automates that procedure. Folds are drawn at the decision-maker level: an individual's entire choice history sits in exactly one fold so that a panel cannot leak between training and hold-out data.

Experimental

The cross-validation utility is functional but still labelled experimental. Expect occasional refinements as I use this functionality in my own work.

Running the sweep

Let's re-use the long-format Apollo frame (with the income_band column) and the formula-based LCLSpec from the estimation tutorial. Passing the same spec keeps the column mapping, the C(alt) constants, the C(income_band) membership terms, and the cost numeraire in one place; num_classes_list then sweeps the class count, overriding spec.classes. Each fold re-fits the formula from scratch, so the categorical base levels are learned per training fold. The search below evaluates two, three, four, and five classes with three-fold CV. To keep the example speedy on a single device, we trim the inner EM loop to twenty-five iterations; in practice you'd loosen this when the optimum is already obvious from a coarse sweep.

import lcl
from lcl import (
    ChoiceIds,
    FitOptions,
    LCLSpec,
    NegativeCoefficient,
    OptimizationOptions,
)

spec = LCLSpec(
    ids=ChoiceIds(alt="alt", case="qID", panel="ID", choice="choice"),
    utility_formula="choice ~ cost + time + C(alt)",
    membership_formula="~ C(income_band) + female",
    constraints={"cost": NegativeCoefficient()},
)

cv_results = lcl.cv_optimal_classes(
    df_long,
    spec,
    num_classes_list=[2, 3, 4, 5],
    folds=3,
    seed=42,
    fit_options=FitOptions(max_em_iter=25, num_devices=1),
    optimization_options=OptimizationOptions(maxiter=30),
)
print(cv_results)

shape: (4, 2)
┌─────────────┬──────────────┐
│ Num_Classes ┆ Avg_OOS_LL   │
│ ---         ┆ ---          │
│ i64         ┆ f64          │
╞═════════════╪══════════════╡
│ 2           ┆ -2178.832849 │
│ 3           ┆ -2148.768403 │
│ 4           ┆ -2147.733992 │
│ 5           ┆ -2148.063559 │
└─────────────┴──────────────┘

The out-of-sample log-likelihood rises sharply from two to three classes, edges up to a peak at four, and then slips back at five. On this richer specification four classes is the clear choice; the decline past four is the telltale sign of a fifth component fitting noise rather than signal.

Plotting the curve

Eye-balling a table of likelihoods is easy when there are four rows, but tougher when you're comparing a dozen candidate Ks. Vega-Altair generates a tidy interactive curve:

import altair as alt
import polars as pl

def plot_cv(cv_df: pl.DataFrame) -> alt.LayerChart:
    optimal = cv_df.filter(pl.col("Avg_OOS_LL") == pl.col("Avg_OOS_LL").max())

    line = (
        alt.Chart(cv_df)
           .mark_line(color="#3F2B47", size=3)
           .encode(
               x=alt.X("Num_Classes:O",
                       title="Number of latent classes",
                       axis=alt.Axis(labelAngle=0)),
               y=alt.Y("Avg_OOS_LL:Q",
                       title="Average out-of-sample log-likelihood",
                       scale=alt.Scale(zero=False)),
           )
    )
    points = line.mark_circle(size=80, color="#3F2B47", opacity=1)
    peak = (
        alt.Chart(optimal)
           .mark_circle(size=160, color="#E37449")
           .encode(x="Num_Classes:O", y="Avg_OOS_LL:Q",
                   tooltip=[
                       alt.Tooltip("Num_Classes:O", title="Classes"),
                       alt.Tooltip("Avg_OOS_LL:Q", title="OOS-LL", format=".2f"),
                   ])
    )
    label = peak.mark_text(align="left", baseline="middle",
                           dx=12, fontSize=12, fontWeight="bold",
                           color="#E37449").encode(text=alt.value("Best K"))

    return (line + points + peak + label).properties(
        title="Held-out log-likelihood by class count",
        width=600, height=380,
    ).configure_title(fontSize=16, anchor="start", offset=20)

plot_cv(cv_results).save("cv_plot.html")

Open cv_plot.html and you have an interactive plot of the curve with the peak highlighted. For a paper figure, swap .save("cv_plot.html") for .save("cv_plot.svg").

Practical notes

  • Stick with a small number of folds for initial screening. Three folds is plenty to see the qualitative shape of the curve. Bump to five or ten once you've narrowed the range.
  • Use the same FitOptions across folds. Differing iteration budgets across folds confound the comparison.
  • Inspect a fold that fails to converge. cv_optimal_classes swallows individual fold failures and records NaN; if you see one, refit that fold with lcl.fit (or results.diagnostics() on the refit) to see what went wrong.
  • Pair CV with the information criteria. When CV picks \(K^*\) and BIC picks \(K^* - 1\), the latter is often the right call for inference; the smaller model trades a small amount of fit for tighter standard errors.