Model selection with CBIC • drmSEM

Use model selection when you have a small set of biologically plausible DAGs and want to ask which graph is best supported after each endogenous node has been fitted with drmTMB. The workflow is confirmatory: you name the candidate DAGs, fit the same data under each candidate, and compare their d-separation evidence. It is not an all-subsets search over every possible graph.

drmSEM now uses CBIC as the default ranking criterion. CBIC is a BIC-style criterion built on Fisher’s C,

$\mathrm{CBIC} = C + k\log(n),$

where $C$ is Fisher’s C statistic from the d-separation test (Shipley 2000, 2009), $k$ is the number of fitted fixed-effect paths, and $n$ is the number of observations. The stronger penalty matters because a harmless over-fitted DAG can have a very small C statistic. CBIC asks a stricter question: is the extra path worth its cost?

Where this construction comes from. The information-criterion family for Shipley’s C statistic was introduced by Shipley (2013) as $\mathrm{CIC} = C + 2k$ , in direct analogy to AIC, with the small-sample correction $\mathrm{CICc} = C + 2k\,n/(n - k - 1)$ ; phylopath (Bijl 2018) defaults to CICc. AIC-style penalties are well known to be inconsistent for true-model recovery — their selection probability for the true model asymptotes to a fixed value below one, even when the truth is in the candidate set and the sample grows (Schwarz 1978; Burnham and Anderson 2002). Combining Shipley’s C with the BIC penalty $k\log(n)$ gives a consistent variant; as far as we can find, this combination has not been named or defaulted to before, so the construction is introduced here as a drmSEM choice. The simulation evidence supporting it is reported in validation.Rmd.

Build a candidate set

Suppose the scientific question is whether temperature acts on fitness directly or through body size. The candidate set might contain a direct model, a mediated model, and a deliberately over-fitted model with both routes.

models <- drm_model_set(
  direct = drm_dag(
    size ~ temp,
    fitness ~ temp
  ),
  mediated = drm_dag(
    size ~ temp,
    fitness ~ size
  ),
  direct_plus_mediated = drm_dag(
    size ~ temp,
    fitness ~ temp + size
  )
)

Each drm_dag() is an unfitted candidate graph. The formulas describe the structural parents of each endogenous node. A fair comparison keeps the response set and data set fixed across candidates, then changes the arrows that represent the competing biological explanations.

Fit and rank the models

compare() fits every candidate with drm_sem(), runs the d-separation test, combines the independence-claim p-values into Fisher’s C, and adds both CBIC and CICc scores.

cmp <- compare(
  models,
  data = dat,
  family = list(
    size = stats::gaussian(),
    fitness = drmTMB::nbinom2()
  )
)

cmp

The printed table contains both criteria:

CBIC, dCBIC, and wCBIC are the default ranking score, delta, and weight.
CICc, dCICc, and wCICc are retained for the phylopath-style support ranking.
weight is the weight for the selected criterion. With the default compare(), it is the same as wCBIC.

The top model is the row with the lowest CBIC. best() returns that fitted SEM.

top <- best(cmp)
paths(top)

avg <- average(cmp)
avg

average() returns criterion-weighted standardized paths. With the default ranking, those are CBIC-weighted paths. A path that appears in only some candidates is averaged conditionally over the models that contain it, with weight_sum showing how much model weight supported that path.

When to ask for CICc

CICc remains available:

cmp_cicc <- compare(
  models,
  data = dat,
  family = list(
    size = stats::gaussian(),
    fitness = drmTMB::nbinom2()
  ),
  criterion = "CICc"
)

Use CICc when you want a phylopath-style support ranking or a sensitivity check. Do not read CICc as the default true-DAG recovery claim. In the validation study below, CICc often retained an over-fitted rival that added a harmless direct path, while CBIC selected the generating DAG.

What the validation showed

The model-selection validation study generated data from a true chain x -> m -> y, compared that truth against an over-fitted rival and a missing-edge rival, and repeated the comparison 300 times at each sample size. The current cached results are:

Model-selection recovery in the cached validation study.
criterion	n	truth selection	mean truth weight	reps	pass
CBIC	300	0.927	0.817	300	TRUE
CBIC	1000	0.973	0.883	300	TRUE
CICc	300	0.610	0.502	300	FALSE
CICc	1000	0.573	0.495	300	FALSE

CBIC passed the recovery gate in this grid: it selected the true DAG in 92.7% of replicates at n = 300 and 97.3% at n = 1000, and it never selected the missing-edge model. CICc also avoided the missing-edge model, but it selected the true DAG only 61.0% and 57.3% of the time because the over-fitted rival kept substantial support. This is a scoped validation claim: it supports CBIC as the default on the tested three-candidate chain, not as a guarantee over arbitrary DAG sets or untested families. That is why CBIC is the default.

A practical reading rule

Read the comparison table in this order:

Define the candidate DAGs before ranking them. compare() is a confirmatory comparison tool, not a graph-discovery routine.
Check whether any candidate has a poor d-separation result. A model with a violated missing-edge claim is usually not a good scientific explanation.
Use CBIC and wCBIC to choose the parsimonious candidate.
Inspect CICc as a support sensitivity check. If CICc gives a harmless over-fitted model non-trivial support, report that as uncertainty about a weak extra path, not as evidence that the extra path is needed.
Use average() when the substantive paths are more important than the identity of a single top model.

This keeps the workflow honest: CBIC answers the “which DAG should I carry forward?” question, while CICc remains a useful view of near-ties and over-fitted alternatives.

References

Bijl, Wouter van der. 2018. “phylopath: Easy Phylogenetic Path Analysis in R.” PeerJ 6: e4718. https://doi.org/10.7717/peerj.4718.

Burnham, Kenneth P., and David R. Anderson. 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. 2nd ed. Springer.

Schwarz, Gideon. 1978. “Estimating the Dimension of a Model.” The Annals of Statistics 6 (2): 461–64. https://doi.org/10.1214/aos/1176344136.

Shipley, Bill. 2000. “A New Inferential Test for Path Models Based on Directed Acyclic Graphs.” Structural Equation Modeling 7 (2): 206–18. https://doi.org/10.1207/S15328007SEM0702_4.

Shipley, Bill. 2009. “Confirmatory Path Analysis in a Generalized Multilevel Context.” Ecology 90 (2): 363–68. https://doi.org/10.1890/08-1034.1.

Shipley, Bill. 2013. “The AIC Model Selection Method Applied to Path Analytic Models Compared Using a d-Sep Test.” Ecology 94 (3): 560–64. https://doi.org/10.1890/12-0976.1.