PolyGenius
Contents

evaluate$compare

Contrast PRS models against a reference model

evaluate$compare() asks how much better or worse each model predicts an outcome than one reference model does, as the paired difference in one performance metric.

Usage

evaluate$compare(
  data,
  outcomes,
  scores.layer = X,
  split.by = NULL,
  metric = NULL,
  reference.model = NULL,
  bootstrap = 1000,
  conf.level = 0.95,
  p.adjust.method = "BH"
)

Arguments

ArgumentDescription
dataA [PolyGeniusStudy](/reference/polygeniusstudy/). Holds the score layer and the sample columns the other arguments name. It is not modified.
outcomesUnquoted outcome expression resolved from data: a bare column, an expression such as status == "case", an [outcome()](/reference/outcome/) descriptor, or a c() or list() of these whose names label the outcomes. The type is binary or continuous, inferred from the values unless the descriptor sets it. A [surv()](/reference/surv/) outcome aborts.
scores.layerUnquoted name of an existing score layer, default X. The layer is used as supplied: nothing is scored or standardised here.
split.byUnquoted expression naming factor, character or logical sample columns, or NULL (default). Adds rows for each level beside the overall stratum = "all" rows. A level named "all" aborts.
metricCharacter vector with at most one value per outcome type, or NULL (default). The metric each contrast takes the difference of. Binary outcomes take "auc" (the default) or "r2.liability". "r2.liability" needs the outcome's prevalence, set with outcome(prevalence =), and aborts without it. Continuous outcomes take "r.squared" (the default). Each value applies to the outcomes of its own type. A type with no value takes its default. A value for a type absent from the call is ignored.
reference.modelNULL (default), one model name, or a character vector of model names named by outcome. NULL takes, per outcome, the model with the highest metric in the "all" stratum. NA values are skipped, and ties go to the first model in layer order. Each model is ranked on its own complete cases. That model is the reference in every split.by level, and $metadata records it as data-selected. A named vector may leave outcomes out; those take the NULL rule. A name matching no outcome, or a model not in the layer, aborts.
bootstrapNon-negative integer, default 1000. Replicates behind the percentile intervals of delta.r2.liability and delta.r.squared. 0 returns estimates without bootstrap intervals. Resampling draws from the caller's random number stream. The DeLong interval of delta.auc does not use it.
conf.levelNumeric scalar in (0, 1), default 0.95.
p.adjust.methodOne of stats::p.adjust.methods, default "BH". Applied within analysis, outcome, metric and stratum.

Value

A PolyGeniusEvaluation following the schema-evaluation schema. $results has the columns listed in evaluate, one row per outcome, non-reference model and stratum, with analysis = "compare". metric is delta.auc, delta.r2.liability or delta.r.squared. estimate is the model's metric minus the reference's, and model.ref names the reference. n and n.events count the pair's complete cases. Only delta.auc rows carry a p-value. There are no artifacts.

$diagnostics$messages is always present, recording scores with the reference's ranks, a stratum with no complete cases or a single class, an outcome left without a reference, a failed DeLong test (NA row), a zero DeLong variance (interval and p-value set to NA), a non-finite metric or DeLong bound (set to NA), each warning from pROC, and the number of failed bootstrap replicates.

$metadata holds the outcome labels, scores.layer, the split.by columns, bootstrap, conf.level and p.adjust.method, and per outcome, named by label, metric, reference.model and reference.selected (TRUE when the reference was data-selected).

Details

Each outcome has one reference model. Every other model is contrasted against it as delta.<metric>: the model's metric minus the reference's, with model.ref naming the reference. The reference itself gets no row. Models are contrasted only against the reference, never with each other. The shared rules on outcomes, strata, intervals and multiple testing are in evaluate.

A contrast uses the complete cases of its pair: the outcome and both models' scores. The score is the only predictor. Both R-squared metrics are the unadjusted R-squared of lm(y ~ score), and r2.liability moves it to the liability scale as evaluate$performance() describes.

delta.auc uses the paired DeLong test (pROC::roc.test()). Its interval is the analytic DeLong interval and its p-value is the DeLong p-value. delta.r2.liability and delta.r.squared use a paired bootstrap: each replicate draws one resample and evaluates both models on it. Their interval is the percentile interval, and they carry no p-value.

A model whose scores have identical ranks to the reference's gets a delta.auc of 0 with lower = upper = 0, an NA p-value, and a diagnostic. A stratum in which a binary outcome has only cases or only controls gives NA rows and a diagnostic. Under reference.model = NULL, an outcome whose metric is NA for every model has no reference, so it gets no rows and a diagnostic. With a named reference, a cell whose metric is NA gives an NA row and a diagnostic. Fewer than two models in the layer aborts.

Interpretation

A reference selected on the same data is the best model there partly by chance. Contrasts against it are biased away from 0, and their intervals and p-values are not adjusted for that selection. This holds for the default reference.model = NULL and for a name picked from the same data. An honest comparison names a reference chosen on independent data, such as a tuning split ranked with evaluate$select().

Examples

# Against the model with the highest AUC in these data.
cmp <- evaluate$compare(study, outcomes = c(case = status == "case"))

# Against a reference chosen on a disjoint tuning split.
ranked <- evaluate$select(
  evaluate$performance(tune, outcomes = c(case = status == "case"))
)
best <- ranked$model[ranked$rank == 1][1]
cmp  <- evaluate$compare(test,
  outcomes        = c(case = status == "case"),
  reference.model = c(case = best))

# Liability scale, within each sex.
cmp <- evaluate$compare(study,
  outcomes = outcome(status, event = "case", prevalence = 0.05),
  metric   = "r2.liability",
  split.by = sex)

See Also

evaluate$select() to choose a reference on independent data. evaluate$profile() runs this together with the other components. visualize$evaluate$compare() plots the contrasts.

Other evaluate-components: evaluate, evaluate.association(), evaluate.incremental(), evaluate.performance(), evaluate.profile(), evaluate.redundancy(), evaluate.select(), evaluate.stratification()

References

DeLong, DeLong and Clarke-Pearson (1988). Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach. Biometrics 44(3), 837-845. doi:10.2307/2531595

Aliases: evaluate.compare, evaluate$compare