Stepwise covariate selection for a parametric hazard model

Run forward, backward, or two-way stepwise selection on an existing hazard fit using score (Q) statistics, Wald p-values, or AIC deltas as the entry / retention criterion. Phase-specific entry is supported for multiphase models: a covariate can enter one phase and not another.

Usage

hzr_stepwise(
  fit,
  scope = NULL,
  data,
  direction = c("both", "forward", "backward"),
  criterion = c("score", "wald", "aic"),
  slentry = 0.3,
  slstay = 0.2,
  max_steps = 50L,
  max_move = 4L,
  force_in = character(),
  force_out = character(),
  trace = TRUE,
  ...
)

# S3 method for class 'hzr_stepwise'
print(x, ...)

# S3 method for class 'hzr_stepwise'
summary(object, ...)

# S3 method for class 'summary.hzr_stepwise'
print(x, ...)

# S3 method for class 'hzr_stepwise'
as.data.frame(x, ...)

Arguments

fit: A fitted hazard object built via the formula = Surv(...) ~ predictors, data = df interface.
scope: Candidate set. NULL (default) uses every data-frame column not already in the model for every phase. A candidate enters as the column it names, whatever its name: a column literally named age:mal enters as that column, not the interaction, and a column "age " beside age enters as itself (#449). For single-distribution fits, pass a one-sided formula (~ age + nyha) or a character vector of names. Each name in a character scope is looked up, not parsed: a name that is exactly a column of data is that column, otherwise a name that is exactly a term label as terms() writes it ("`_X1`", "log(age)", "age:mal") is that term, and any other name is ignored with a warning naming it. The column is looked up first, so when data has a column literally named age:mal, "age:mal" puts that column in the scope rather than the interaction. A candidate enters as the column or term its name resolved to, so a bare "_X1" enters the column _X1 (#449, #441), and the score reads the values of that column however the name was written (#438). For multiphase fits, pass a named list of one-sided formulas keyed by phase, naming each phase once. scope lists what may enter; a drop considers every term in the model except force_in and terms frozen by max_move at an earlier iteration (see the The frozen set section). A two-sided formula is an error, since its left-hand side would never be a candidate, and so is a non-empty scope under direction = "backward", which does not read it. An empty scope (~ 1, character(), or a list of NULLs and ~ 1s) is accepted there.
data: Data frame the base fit was built on. Required for refits. Its rows must be the fit's rows in the same order: the score criterion reads each candidate by position, so a sorted or reordered frame is refused when the fit stored the frame it was given. On the vector interface, and for any fit with weights, every criterion's refits pair data with vectors stored in the fit's row order, so there the frame is checked whatever the criterion. Columns added after the fit, such as derived candidates, are allowed. A vector-interface fit made without data =, or with a data of another row count used only to look names up, stores no frame to compare, and a data holding only columns added after the fit shares none; and shared columns with duplicate rows cannot show rows reordered among those duplicates. For those the order is checked against a column of data holding the fit's event times, and refused if they are out of order. Rows reordered within a tie leave that column unchanged, so when the times have ties the fit's other per-row inputs (status, interval bounds, weights and covariate design columns) are looked for in data too, each under the column it came from: the column the call named (status = dead, status = d$dead, or the event in Surv(time, event)), or a covariate's own name. Rows reordered within a tie are refused. The order is accepted when the inputs found tell every row apart, or when all of them are found: rows that could still be swapped are then identical in everything the fit reads, and the screen's answer is the same. With no column of event times, or tied times the other inputs cannot resolve, the order cannot be checked, and a warning of class hzr_score_rows_unverified says so and names the inputs to add. Column names must be unique: every check reads data by name, so a data with duplicated names is refused.
direction: Search strategy: one of "both" (default), "forward", or "backward". Controls whether variables may only enter, only leave, or both. See the Selection direction and criterion section.
criterion: Entry / retention rule: one of "score" (default), "wald", or "aic". "score" and "wald" both apply SAS-style p-value thresholds (slentry / slstay) but score entry candidates differently, and can therefore select different variable sets; "score" reproduces C/SAS HAZARD and needs no per-candidate refit, except at a step whose information cannot be formed, where it falls back to refitting and Wald-testing each candidate. "aic" adds or drops whenever it lowers the AIC. "score" and "aic" refuse a base fit that did not converge, since both measure candidates against it, and so does "wald" unless direction = "forward", since a removal is tested on the base's own estimates. See the Selection direction and criterion section.
slentry: Entry p-value threshold for the score / Wald criteria. Default 0.30 matches SAS SLENTRY.
slstay: Retention p-value threshold for the score / Wald criteria. Default 0.20 matches SAS SLSTAY.
max_steps: Hard cap on total accepted actions. Emits a warning() if hit. Default 50.
max_move: Per-variable oscillation cap. When a variable has entered + exited more than max_move times it is frozen for the remainder of the run, in the state it ends that iteration in. Default 4. See the The frozen set section.
force_in: Character vector of variables that must remain in the model. Such variables are still scored and reported in the selection trace, but are never dropped. Each name is looked up, not parsed, once, when the screen starts: a name that is exactly a column of data is that column, so the bare "_X1" pins the column _X1 although terms() labels it `_X1`, and "TRUE" pins a column named TRUE. The exception is a column no model formula can hold, "." or "": it can never become a model term, so it cannot be pinned, and naming it warns and lists it in unresolved. Otherwise a name that is exactly a term label of the model or scope is that term, so "`_X1`" and "age:mal" work too. Any other name, "age " with a trailing space when there is no such column, say, matches nothing and is ignored with a warning naming it. The column is looked up first: when data has a column literally named age:mal, "age:mal" resolves to that column and not to the interaction.
force_out: Character vector of variables that may never be considered as candidates. Names are looked up as for force_in: a column of data first, then a term label of the model or scope, and a warning for a name that is neither.
trace: Logical; print step-by-step progress to the console. Default TRUE.
...: Passed to every candidate refit. Only control (e.g. control = list(maxit = 500)) and an objective equal to the base fit's are accepted. Any other name is an error: a misspelling such as slentyr would be stored by hazard() without being read, and every other hazard() argument (the response, weights, time_windows, dist, theta, phases, fit and so on) is set by the refit itself, from the base model and data, so that each candidate is compared with the model it extends. The print(), summary() and as.data.frame() methods ignore ....
x: An hzr_stepwise object.
object: An hzr_stepwise object.

Value

An object of class c("hzr_stepwise", "hazard"), the final fit augmented with:

steps: Data frame with one row per accepted / frozen action; see Details.
scope: Record of the candidate scope, plus force_in, force_out, and the frozen set. Four fields answer four different questions about the pins, and none substitutes for another: force_in and force_out are the names as given, including any that pinned nothing; force_in_resolved and force_out_resolved are what those names resolved to, as a column or term label, so a bare "_X1" appears there as "`_X1`"; unresolved is which of them resolved to nothing; and frozen is what the loop held in that the caller never named. The resolved fields carry the resolved names only: the frozen set is not merged into them, because frozen already records it and merging would list variables nobody asked for (#451). Resolving is not applying: a name that resolves to a variable the model does not contain is recorded here and still pins nothing, because force_in only keeps a variable that is already in. The resolved fields are not aligned element for element with the as-given ones, which stay longer by every name that resolved to nothing. A variable in frozen is in the final model if it was in the model when it was frozen, and out of it otherwise; see the The frozen set section. unresolved is a list with elements force_in, force_out and scope, each the names that could not be used and were therefore ignored (character() when none were). That is usually a name matching neither a column of data nor a term label; for force_in and force_out it also covers a name that IS a column but which no model formula can hold, such as "." or "", since such a column can never become a model term and so can never be pinned. The trace, and so print() and summary(), carries a line for each non-empty one, and a screen whose character scope was emptied this way says so where it stops. candidates, force_in and force_out here are the arguments as given, so a character candidates still lists the names that were ignored, as force_in and force_out do; read unresolved for which those were (#451).
criteria: Named list of the threshold / direction settings actually applied, plus n_uncomputable_scores (how many candidate scores were NA, counted once per step: an entry the score test could not score under criterion = "score"; an entry under criterion = "wald", or a removal under any criterion, whose Wald statistic could not be computed for want of a variance, wald_no_variance; or an entry under criterion = "aic" whose fit had no finite objective, nonfinite, or was fitted on different rows from the current model because the candidate is missing on some, rows_differ; or an entry under criterion = "aic" or "wald", or rescued by the Wald fallback under "score", whose refit ended below the current model's log-likelihood, which it contains, so that the refit cannot have converged, loglik_below_base), uncomputable_reasons (a named integer vector of why), wald_untested_removals and wald_untested_entries (the "var" / "var@phase" tokens of variables kept in, or left out, on a step whose Wald test for them could not be computed; a variable tested, or refused as loglik_below_base, at a later step is not listed. Entries are listed under criterion = "wald" only: under "score" an entry no test could reach is reported by its reason, such as fallback_no_variance) and stopped_uncomputable (TRUE when the last iteration had candidates for entry or for removal and could test none of them). Read uncomputable_reasons before treating an unscored candidate as a bad one: information_indefinite marks candidates whose effect is too large for the score test's approximation at zero, which are typically the strongest variables on offer rather than degenerate ones. Candidates with that cause, or with coefficient_diverging, information_nonpositive or nuisance_singular, are refit and tested by Wald automatically, counted in n_wald_fallbacks. nuisance_singular is a fault of the score test at the current model, whose information matrix could not be formed or inverted, so it applies to every candidate at that step and each of them is refitted. A candidate still reaches uncomputable_reasons when that refit fails, or when its cause is any other, for which no refit is attempted. Read uncomputable_reasons for which one it was in any given run. For every criterion it also carries refit_failures (the "var" / "var@phase" tokens of candidate moves whose refit errored, failed to converge, or was refused because the move would not change the model – a drop that removes no design column, #320), refit_failure_reasons (why each one failed or was refused: the refit's error message, that it did not converge, or that the move changes nothing; named by the same tokens), n_refit_failures, and stopped_refit_failed (TRUE when the run ended on an iteration in which a refit failed or a move was refused. A refit failure is a screen that could not test its candidates, rather than one that tested them and liked none; a refusal is determinate – the move was tested and would have left the model as it was. Read refit_failure_reasons for which it was). Check it before reading a zero-row steps as an honest null result.
trace_msg: Character vector of the trace lines, captured regardless of the trace flag.
elapsed: difftime from start to finish.
final_call: The call that produced this result.

print.hzr_stepwise returns x invisibly.

summary.hzr_stepwise returns a summary.hzr_stepwise object (extends summary.hazard) with $stepwise_steps and $stepwise_trace appended.

print.summary.hzr_stepwise returns x invisibly.

as.data.frame.hzr_stepwise returns the $steps data frame.

Details

The steps data frame has columns:

step_num: Integer sequence starting at 1.
action: "enter", "drop", or "frozen".
variable: The term entered, dropped or frozen, as the model's terms() label: a non-syntactic column _X1 is recorded as `_X1` on every row, whether the scope named it "_X1", "`_X1`" or in a formula.
phase: Phase name (multiphase) or NA_character_.
criterion: The criterion actually applied to this step: "score", "wald", or "aic". Under criterion = "score" the drop rows read "wald", because score is entry-only.
score: Winning score used for the decision.
stat, df: Test statistic and degrees of freedom.
stat_type: What stat is on this row, and so which reference distribution recomputes its p-value: "score_q" (chi-square on df), "wald_z" (standard normal) or "wald_chisq" (chi-square on df). df alone does not distinguish them; a scalar Wald is reported as a z, not as its square, so it and a score Q are both recorded at df = 1 while calling for different distributions. It also identifies the rows the Wald fallback rescued, but only among entry rows: under criterion = "score" those are the rows with action == "enter" & stat_type == "wald_z". Drop rows are always Wald-tested under that criterion (removal follows SAS and is tested on the current model's Wald p-value), so they read "wald_z" whether or not the fallback ever fired. See $criteria$n_wald_fallbacks for the count.
p_value, delta_aic: Always populated when computable, regardless of the active criterion.
logLik, aic, n_coef: Goodness-of-fit diagnostics of the model after this step. Under objective = "sas", on a step whose model reads an interval-censored row (one of positive weight, not dropped by a phase design), logLik holds PROC HAZARD's objective, not a log-likelihood, and aic is computed from it. A step can change which rows are read (see n_rows), so on a step that starts or stops reading such rows delta_logLik is the difference of two different quantities. The trace's final line labels the final model's value, and criterion = "aic" warns once (class hzr_stepwise_sas_objective) as soon as a step's model reads such a row.
n_rows: Number of rows in the fit's data after this step (rows given weight 0 are counted). A multiphase fit drops every row where a variable in the model is missing, so entering a variable with missing values shrinks the sample every later step is tested on, and dropping it grows the sample back. A step that changes the rows of positive weight raises a warning of class hzr_stepwise_sample_changed naming the step, the variable and the row counts before and after. Which variables enter is not changed by this: compare n_rows across steps to see how many rows each test rested on.

Selection direction and criterion

Two arguments shape the search. direction decides which moves are allowed at each step; criterion decides how a candidate move is scored and whether it is accepted.

direction = "forward": Start from the base model and only add variables; the best eligible candidate enters each step until none clears the entry rule. Variables never leave once in.
direction = "backward": Start from the base model, which must already hold every candidate, and only drop variables; the weakest term leaves each step until all survivors clear the retention rule. scope is not read, so a non-empty one is an error: protect terms with force_in.
direction = "both" (default): Two-way stepwise: on every iteration, whether or not a variable entered, every term in the model, the base model's included, is re-tested and may be dropped unless it is in force_in or was frozen by max_move at an earlier iteration (see the The frozen set section). scope limits what may enter, not what may leave. This is the SAS SELECTION = STEPWISE strategy. max_move caps how often a single variable may oscillate before it is frozen.

criterion = "score" (default)

Accept moves on SAS-style significance thresholds, using the score (Q) statistic of the candidate coefficient; this reproduces C/SAS HAZARD's SELECTION statistic. Q is evaluated at the current model's MLE with the candidate's coefficient pinned at zero, so no candidate refit is needed: the reduced-model information is inverted once per step and reused across every candidate. Only the winner is refit. A candidate enters if its p-value is below slentry.

Score is an entry criterion; the drop path never refit per candidate in the first place, so removals are tested on the current model's Wald p-value against slstay, as SAS does.

For single-distribution fits, the score criterion computes the observed information numerically via the suggested numDeriv package and errors with a clear message if it is not installed; a multiphase fit uses the analytic Hessian instead and does not need it.

Following SAS, the variance used during selection is approximate: shaping-parameter covariances are ignored. This affects selection only; final-model standard errors are unchanged and still come from the full Hessian. Candidates must be single-column numeric main-effect terms; a factor is rejected with an error rather than skipped.

criterion = "wald"

Accept moves on SAS-style significance thresholds, using the Wald $\chi^2$ of the affected coefficient(s): a candidate enters if its p-value is below slentry, and a term is dropped if its p-value rises above slstay. Entry candidates are scored from a refit that adds the candidate (so its new coefficient can be tested); drop candidates are scored from the current model's Wald p-values without a per-candidate refit, and a single refit is run only after a drop is chosen. This was the default before version 1.2.0. It differs algorithmically from C/SAS HAZARD, so the two criteria can take different step paths (and select different variable sets) even when they converge to a similar final model.

criterion = "aic"

Accept any move with $\Delta\mathrm{AIC} < 0$ (a strictly better penalised fit), ignoring slentry / slstay. Entry candidates use the actual $\Delta\mathrm{AIC}$ from the candidate refit; drop candidates use a Wald-to-likelihood-ratio approximation, $\Delta\mathrm{AIC} \approx W - 2\,\mathrm{df}$, computed from the current model without a per-candidate refit (the chosen drop is refit afterwards). Use this for a non-significance-based, information-criterion search.

The frozen set

A variable that has moved (entered or left) more than max_move times is frozen: it is held in the state it is in, in the model or out of it, for the rest of the run, and listed in $scope$frozen.

The freeze takes effect at the end of the iteration in which the variable reached the cap. A two-way iteration (direction = "both") makes a forward step and then a backward step, so a variable that reaches the cap by entering and is dropped by the same iteration's backward step is frozen out. Its "frozen" row in $steps follows the "drop". A variable is therefore in the final model exactly when it was in the model when it was frozen, and $scope$frozen agrees with the final model.

Earlier versions froze such a variable as it entered and then dropped it anyway, so $scope$frozen could name a variable the final model did not contain (#580).

This is not PROC HAZARD's MOVE rule. PROC HAZARD makes one move per step, a removal before any entry; it counts a variable's exits (and, under NOSTEPWISE, its entries too), separately for each phase; and a variable at its limit can neither leave nor enter. Outside NOSTEPWISE it therefore freezes a variable only as it leaves. max_move counts entries and exits alike, by name across phases, which is why the SAS translator records MOVE= rather than mapping it. Adopting PROC HAZARD's rule would change which variables are selected, and is deferred; see https://github.com/ehrlinger/TemporalHazard/issues/378 and https://github.com/ehrlinger/TemporalHazard/issues/379.

Examples

data(avc)
avc <- na.omit(avc)
base <- hazard(survival::Surv(int_dead, dead) ~ age,
               data = avc, dist = "weibull", fit = TRUE,
               theta = c(mu = 0.01, nu = 0.5, 0))
# \donttest{
sw <- hzr_stepwise(base, scope = ~ age + mal,
                   data = avc, direction = "forward",
                   control = list(n_starts = 1))
#> Warning: 'control' element(s) with no effect on this dist = "weibull" fit, ignored: control$n_starts (it applies only to dist = "multiphase").
#> Stepwise selection (direction = forward, criterion = score, slentry = 0.30, slstay = 0.20)
#> 
#> Step 1: ENTER  mal   (p = 0.001)
#> (no further action after 1 step)
#> 
#> Final model: 2 covariates, logLik = -218.60, AIC = 445.19
print(sw)
#> Stepwise selection (direction = forward, criterion = score, slentry = 0.30, slstay = 0.20)
#> 
#> Step 1: ENTER  mal   (p = 0.001)
#> (no further action after 1 step)
#> 
#> Final model: 2 covariates, logLik = -218.60, AIC = 445.19
# }