Fit and rank a set of candidate causal models

compare() fits every candidate drm_dag() in a drm_model_set() with drm_sem() (one node per response, using a per-node or shared family), runs dsep()/fisher_c() on each, and ranks the candidates by an information criterion built on Fisher's C. This is drmSEM's analogue of phylopath's phylo_path() (van der Bijl 2018).

Usage

compare(object, data, family = NULL, criterion = c("CBIC", "CICc"), ...)

# S3 method for class 'drm_model_set'
compare(object, data, family = NULL, criterion = c("CBIC", "CICc"), ...)

Arguments

object: A drm_model_set.
data: A data frame supplied to every node fit.
family: Optional. A single drmTMB family applied to every node, or a named list of families keyed by node (response) name. NULL uses the drm_node() default (Gaussian) for every node.
criterion: Ranking criterion. "CBIC" is the default BIC-style parsimony criterion; "CICc" gives the phylopath-style support ranking.
...: Further arguments passed on to each drm_sem() fit.

Value

A drm_compare object: a data frame with one row per candidate and columns model, fisher_c, df, p.value, k, n, CICc, dCICc, wCICc, CBIC, dCBIC, wCBIC, and weight, sorted by the selected criterion. weight is the weight for the selected criterion, so average() follows the same ranking rule. The fitted drm_sem objects and per-model d-separation tables are carried in attributes for best()/average().

Details

The default ranking statistic is $$\mathrm{CBIC} = C + k\log(n),$$ where $C$ is Fisher's C statistic from the model's d-separation test, $k$ is the number of estimated fixed-effect coefficients across all nodes (counted via paths()), and $n$ is the number of observations. CBIC uses a BIC-style penalty and is the default because it is more conservative about extra paths. criterion = "CICc" instead ranks by $$\mathrm{CICc} = C + 2k\,\frac{n}{n - k - 1},$$ $\Delta$CICc is the difference from the best (lowest-CICc) model and the CICc weight is $\exp(-\tfrac12\,\Delta\mathrm{CICc})$, normalised to sum to one. As $n \to \infty$ the correction term tends to $2k$, so CICc reduces to the usual $C + 2k$.

Citation and novelty. CICc is the C-statistic information criterion of Shipley (2013), used as the default in phylopath (van der Bijl 2018). CBIC is introduced here as a BIC-style analogue: it combines Shipley's C statistic with the BIC penalty $k\log(n)$ (Schwarz 1978). As far as we can find, no published criterion of this exact form has been named or defaulted to in the piecewise-SEM / d-separation literature; users who want to stay on the established AIC-style ranking should set criterion = "CICc". Information-criterion weights for multimodel inference follow Burnham & Anderson (2002).

References

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

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

Burnham KP, Anderson DR (2002). Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach, 2nd edition. Springer, New York.

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

Lefcheck JS (2016). “piecewiseSEM: Piecewise Structural Equation Modelling in R for Ecology, Evolution, and Systematics.” Methods in Ecology and Evolution, 7(5), 573–579. doi:10.1111/2041-210X.12512 .

Examples

if (FALSE) { # \dontrun{
models <- drm_model_set(
  direct   = drm_dag(fitness ~ temp + size),
  mediated = drm_dag(size ~ temp, fitness ~ size + temp)
)
cmp <- compare(models, data = dat,
               family = list(size = stats::gaussian(),
                             fitness = stats::gaussian()))
cmp
best(cmp)
average(cmp)
} # }