--- title: "Migrating from SAS HAZARD to TemporalHazard" format: html vignette: > %\VignetteIndexEntry{Migrating from SAS HAZARD to TemporalHazard} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} --- ```{r} #| include: false # R CMD build installs the package into R_BUILD_TEMPLIB. A build launched # through callr (rcmdcheck, devtools) passes Quarto's R session a profile # that resets .libPaths() without it, and every chunk skipped (#276). templib <- Sys.getenv("R_BUILD_TEMPLIB") if (dir.exists(templib)) .libPaths(c(templib, .libPaths())) has_pkg <- requireNamespace("TemporalHazard", quietly = TRUE) # On CI a package that cannot load must fail the build, not skip every # chunk and pass over nothing (#276). loadNamespace() says why it failed. if (nzchar(Sys.getenv("CI"))) loadNamespace("TemporalHazard") knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = has_pkg ) ``` ```{r} #| label: setup library(TemporalHazard) ``` ## Overview If you've been running multiphase hazard analyses in SAS HAZARD (the macro suite that wraps the C HAZARD binary written by Blackstone, Naftel, and Turner at UAB), this vignette is the bridge to running the same analyses in R. Every SAS `HAZARD` statement, every common macro option, and every output field maps to something in `TemporalHazard`; this document spells out the correspondences and flags the small handful of features the R port doesn't yet cover. The mapping is faithful, not literal. SAS HAZARD is a step-driven DSL: you write a `PROC HAZARD` block with separate `TIME`, `EVENT`, `PARMS`, `EARLY`, `CONSTANT`, and `SELECTION` statements, and the macro assembles them into a binary call. R is function-driven: you write a single `hazard()` call with named arguments. The conceptual pieces are identical; the syntactic ergonomics differ. Use this vignette to translate an existing SAS analysis, or as a reference when comparing the package's output against a SAS HAZARD reference fit for parity testing. The full formal argument mapping table is available programmatically and ships with the package. It is handy when you want to grep for a specific SAS parameter name and find its R equivalent without scrolling through the prose below: ```{r} #| label: mapping-table knitr::kable( TemporalHazard::hzr_argument_mapping(), caption = "Formal argument map: SAS HAZARD/C → hazard()", col.names = c( "SAS Statement", "Legacy Input", "C Concept", "R Parameter", "Required", "Expected Type", "Transform Rule", "Status", "Notes" ) ) ``` --- ## Statement-by-statement mapping A SAS HAZARD analysis is a stack of statements inside a `PROC HAZARD` block. We walk through them in roughly the order they appear in a typical analysis script, showing the SAS form and the corresponding R call. ### `PROC HAZARD DATA=` The procedure-level statement that names the dataset and sets global options. The `NOCOV` / `NOCOR` flags suppress covariance and correlation output; `CONDITION=` is a tolerance switch for the convergence check. ```sas PROC HAZARD DATA=AVCS NOCOV NOCOR CONDITION=14; ``` `DATA=` becomes `hazard()`'s `data` argument. None of the three options has a `control` equivalent. `NOCOV` and `NOCOR` only suppress printed output, and `hazard()` prints nothing until you ask for it with `summary()` or `vcov()`. `CONDITION=` sets a tolerance that `hazard()` does not take: it reports the Hessian's conditioning through `rcond` instead. So the statement translates to: ```{r} #| eval: false fit <- hazard(..., data = avcs) ``` The `control` settings a fit reads are `maxit` and `reltol` for every model, plus `n_starts`, `conserve`, `phase_share_tol` and `start_seed` for a multiphase one. If you pass `nocov`, `nocor` or `condition` anyway, `hazard()` warns that the option does nothing and fits the model unchanged. A name that nothing reads, such as a misspelling, draws the same warning, so a mistyped setting cannot silently leave the default in force. Additional `PROC HAZARD` options with no direct R equivalent yet are passed through `...` as named arguments and stored in `fit$legacy_args`. --- ### `TIME` Names the follow-up time variable. In SAS HAZARD this is a separate statement; in R it is the `time` argument of the vector interface shown here, or the first argument to `Surv()` if you use the formula interface. ```sas TIME INT_DEAD; ``` Maps directly to `time`: ```{r} #| eval: false fit <- hazard( time = avcs$INT_DEAD, ... ) ``` - Must be a non-negative numeric vector in the same time unit as the original analysis (months, in the AVC example). - Missing values are not permitted. --- ### `EVENT` Names the event-indicator variable. Like `TIME`, this is a separate statement in SAS HAZARD; in R it is the `status` argument shown here, or the second argument to `Surv()` in the formula interface. ```sas EVENT DEAD; ``` Maps to `status`: ```{r} #| eval: false fit <- hazard( status = avcs$DEAD, # 1 = event, 0 = censored ... ) ``` - `TemporalHazard` coerces `status` to `numeric`; logical vectors are also accepted. - The HAZARD convention uses `1` = occurred, `0` = censored/alive at last contact. --- ### `PARMS` Supplies starting values for the optimizer and flags which parameters to hold fixed (`FIXM`, `FIXMU`, etc.). The starting values matter. For the multiphase optimizer in particular, a poor starting point can park the fit at a local minimum well away from the global MLE. ```sas PARMS MUE=0.3504743 THALF=0.1905077 NU=1.437416 M=1 FIXM MUC=4.391673E-07; ``` Maps to the phases in `phases`, the start vector `theta`, and, for `FIXM`, `fixed = "m"` on the early phase: ```{r} #| eval: false fit <- hazard( ..., dist = "multiphase", phases = list( early = hzr_phase("cdf", t_half = 0.1905077, nu = 1.437416, m = 1, fixed = "m"), # FIXM: hold M at its start constant = hzr_phase("constant") ), theta = c(log(0.3504743), # MUE log(0.1905077), 1.437416, 1, # THALF, NU, M log(4.391673e-07)) # MUC ) ``` > **Note:** `theta` takes the `PARMS` values on the scale the optimizer > works on, which is not always the scale SAS prints. `MUE`, `MUC` and > `THALF` go in as logs; `NU` and `M` go in as written. A fixed parameter is > declared on the phase that owns it, so `FIXM` is `fixed = "m"` inside the > early `hzr_phase()`. `hazard()` does not read a `control$fix` entry: a fix > written there leaves `M` free, and `hazard()` warns that it never reads > `control$fix`. --- ### `EARLY` and `CONSTANT` covariate blocks These statements assign covariates to specific phases. In SAS HAZARD each block lists the covariates that affect that phase and their starting coefficients. Covariates can appear in multiple blocks with different starting values (same column, different phase-specific effect). ```sas EARLY AGE=-0.03205774, COM_IV=1.336675, MAL=0.6872028, OPMOS=-0.01963377, OP_AGE=0.0002086689, STATUS=0.5169533; CONSTANT INC_SURG=1.375285, ORIFICE=3.11765, STATUS=1.054988; ``` In `TemporalHazard` each block becomes the `formula` of the phase it names: the `EARLY` list goes on the early `hzr_phase("cdf", ...)`, the `CONSTANT` list on `hzr_phase("constant", ...)`. A covariate listed in both blocks, `STATUS` here, appears in both formulas and gets a separate coefficient in each phase, as it does in SAS. The model formula's right-hand side stays `~ 1` because every covariate belongs to a phase. The starting coefficients go into `theta` in phase order: the early phase's scale and shapes, then its covariates, then the constant phase's scale, then its covariates. In this chunk `avcs` is `data(avc)` with its missing `inc_surg` values filled, as in the worked example below. ```{r} #| eval: false fit <- hazard( Surv(int_dead, dead) ~ 1, data = avcs, dist = "multiphase", phases = list( early = hzr_phase("cdf", t_half = 0.1905077, nu = 1.437416, m = 1, formula = ~ age + com_iv + mal + opmos + op_age + status, fixed = "m"), constant = hzr_phase("constant", formula = ~ inc_surg + orifice + status) ), theta = c( log(0.3504743), log(0.1905077), 1.437416, 1, # MUE, THALF, NU, M -0.03205774, 1.336675, 0.6872028, # EARLY: AGE, COM_IV, MAL, -0.01963377, 0.0002086689, 0.5169533, # OPMOS, OP_AGE, STATUS log(4.391673e-07), # MUC 1.375285, 3.11765, 1.054988 # CONSTANT: INC_SURG, ), # ORIFICE, STATUS fit = TRUE ) ``` --- ### `SELECTION` Triggers stepwise covariate selection inside the hazard fit. `SLE` is the significance level for entry, `SLS` for staying. SAS HAZARD embedded this in the same `PROC HAZARD` call; in R the stepwise loop is a separate `hzr_stepwise()` function that wraps a fitted `hazard()` object. ```sas SELECTION SLE=0.2 SLS=0.1; ``` `TemporalHazard` now ships `hzr_stepwise()`, which implements forward, backward, and two-way stepwise selection with SAS-style SLENTRY / SLSTAY thresholds, phase-specific entry for multiphase models, and a MOVE oscillation guard. It now defaults to `criterion = "score"`, which reproduces SAS's `SELECTION` Q statistic (a score test at the current estimates, with no per-candidate refit except where the information cannot be formed), so entry decisions match SAS's directly. Passing `criterion = "wald"` selects the older refit-based path, which fits each candidate before testing it. See `?hzr_stepwise` for the full option list. --- ## SAS macro equivalents The legacy SAS distribution ships a suite of macros for nonparametric baselines, goodness-of-fit diagnostics, bootstrap inference, and variable calibration that sit alongside `PROC HAZARD`. Each has an R equivalent in `TemporalHazard`: | SAS macro | R function | Purpose | |:---|:---|:---| | `kaplan.sas` | `hzr_kaplan()` | KM survival with logit-transformed exact CL | | `nelsonl.sas` / `nelsont.sas` | `hzr_nelson()` | Nelson cumulative hazard with lognormal CL | | `hazplot.sas` | `hzr_gof()` | Parametric vs. KM overlay + CoE goodness-of-fit | | `deciles.hazard.sas` | `hzr_deciles()` | Decile-of-risk calibration + chi-square GOF | | `logit.sas` / `logitgr.sas` | `hzr_calibrate()` | Quantile-group calibration (logit / Gompertz / Cox link) | | `bootstrap.hazard.sas` | `hzr_bootstrap()` | Bootstrap resampling + summary | | `markov.sas` | `hzr_competing_risks()` | Aalen-Johansen competing-risks incidence | Minimal illustrations on the AVC dataset follow. The `vignette("inference-diagnostics")` walkthrough builds these into a full analysis workflow. ```{r} #| label: macro-demos-setup library(survival) data(avc) avc <- na.omit(avc) # Base parametric fit for downstream diagnostics fit <- hazard( Surv(int_dead, dead) ~ age + status + mal + com_iv, data = avc, dist = "weibull", theta = c(mu = 0.2, nu = 1, rep(0, 4)), fit = TRUE, control = list(maxit = 500) ) ``` ```{r} #| label: kaplan-demo km <- hzr_kaplan(time = avc$int_dead, status = avc$dead) head(km, 4) ``` ```{r} #| label: nelson-demo nel <- hzr_nelson(time = avc$int_dead, event = avc$dead) head(nel, 4) ``` ```{r} #| label: gof-demo head(hzr_gof(fit), 4) ``` ```{r} #| label: deciles-demo hzr_deciles(fit, time = 120, groups = 10) ``` ```{r} #| label: calibrate-demo hzr_calibrate(avc$age, avc$dead, groups = 10, link = "logit") ``` ```{r} #| label: bootstrap-demo set.seed(1) hzr_bootstrap(fit, n_boot = 20) # small for vignette build ``` ```{r} #| label: competing-demo data(valves) # Synthesize a competing-risks indicator for illustration. valves$ev <- with(valves, ifelse(dead == 1, 1L, ifelse(pve == 1, 2L, ifelse(reop == 1, 3L, 0L)))) head(hzr_competing_risks(valves$int_dead, valves$ev), 4) ``` --- ## Full worked example: AVC death after repair Statement-by-statement mapping is fine for reference, but the full gestalt only lands when you see a complete SAS HAZARD analysis and its R translation side by side. The example below is the final multivariable model from `examples/hm.death.AVC.sas` in the reference C repository: death after atrioventricular canal repair, two-phase model with covariates assigned to specific phases. It's the canonical "this is what a real SAS HAZARD analysis looks like" specimen. ### SAS (original) The original SAS code, exactly as it appears in the reference distribution. Notice how the statements stack inside a single `%HAZARD()` macro call. ```sas %HAZARD( PROC HAZARD DATA=AVCS P CONSERVE OUTHAZ=OUTEST CONDITION=14 QUASI; TIME INT_DEAD; EVENT DEAD; PARMS MUE=0.3504743 THALF=0.1905077 NU=1.437416 M=1 FIXM MUC=4.391673E-07; EARLY AGE=-0.03205774, COM_IV=1.336675, MAL=0.6872028, OPMOS=-0.01963377, OP_AGE=0.0002086689, STATUS=0.5169533; CONSTANT INC_SURG=1.375285, ORIFICE=3.11765, STATUS=1.054988; ); ``` ### R equivalent The same model in `TemporalHazard`, with one R argument for each SAS statement: - `TIME INT_DEAD` and `EVENT DEAD` become `Surv(int_dead, dead)` on the left of the formula. In `avc`, `dead` is 0/1, so it is a status as it stands. A job whose `EVENT` variable counts events also needs `weights`; see the translator section below. - `PARMS` becomes the phase shapes and the `theta` start vector, and `FIXM` becomes `fixed = "m"` on the early phase. - `EARLY` and `CONSTANT` become the `formula` of the matching `hzr_phase()`. - `CONSERVE` goes into `control` as `conserve = TRUE`. `CONDITION=14` and `QUASI` have no `hazard()` equivalent, so nothing is written for them: `CONDITION=` stops `PROC HAZARD`'s optimizer when its Hessian approximation becomes too ill-conditioned, and `QUASI` chooses its optimizer. `hazard()` offers no optimizer choice (it fits by BFGS, after a Nelder-Mead warm-up in some multiphase fits) and has no such stop. Passed anyway, as `control = list(condition = 14, method = "bfgs")`, they would only draw a warning that they do nothing. - `PROC STANDARD REPLACE` from the SAS `DATA` steps becomes one line that fills the missing `inc_surg` values with the column mean. Once you have translated one analysis this way, the pattern carries to every other. ```{r} #| label: avc-example #| warning: false data(avc) # reload: the diagnostics above dropped the incomplete rows avcs <- avc # PROC STANDARD REPLACE: fill missing INC_SURG with the column mean avcs$inc_surg[is.na(avcs$inc_surg)] <- mean(avcs$inc_surg, na.rm = TRUE) fit <- hazard( Surv(int_dead, dead) ~ 1, # TIME, EVENT data = avcs, dist = "multiphase", phases = list( early = hzr_phase("cdf", t_half = 0.1905077, nu = 1.437416, m = 1, formula = ~ age + com_iv + mal + opmos + op_age + status, # EARLY fixed = "m"), # FIXM constant = hzr_phase("constant", formula = ~ inc_surg + orifice + status) # CONSTANT ), theta = c( # PARMS, EARLY, CONSTANT log(0.3504743), log(0.1905077), 1.437416, 1, -0.03205774, 1.336675, 0.6872028, -0.01963377, 0.0002086689, 0.5169533, log(4.391673e-07), 1.375285, 3.11765, 1.054988 ), control = list(conserve = TRUE), # CONSERVE fit = TRUE ) summary(fit) ``` The `PARMS` and covariate values in this job are the estimates from an earlier SAS run of the same model, so both programs start at SAS's optimum. R agrees that it is one. The log-likelihood is -160.408, the value SAS prints for this model, and the standard errors agree with the SAS listing to within 0.31%: `AGE` is 0.009381 here against SAS's 0.009380, and the early-phase `log_mu` (SAS's `E0`) is 0.7527 against 0.7550. `summary()` notes an ill-conditioned Hessian (`rcond = 4.41e-12`); SAS reports the same condition for this fit as a log10 condition code of 11.12. Run interactively, the fit also raises `Hessian is ill-conditioned` warnings. The chunk hides them because `summary()` reports the same thing. --- ## Data preparation differences | SAS HAZARD | TemporalHazard R | |----------------------------------------|---------------------------------------------------| | `PROC STANDARD REPLACE` for missing | Replace with `mean(..., na.rm = TRUE)` manually | | Log transforms in `DATA` step | `transform()` or `dplyr::mutate()` | | Variable labels via `LABEL` | Column names of `x` matrix carry through | | `FILENAME` / `INFILE` / `INPUT` | `read.table()`, `read.csv()`, or `haven::read_sas()` | | SAS formats (date, currency, etc.) | Standard R numeric; apply `as.Date()` if needed | --- ## Output object `hazard()` returns a list of class `hazard`: | Slot | Contents | |-----------------------|--------------------------------------------------------------| | `$call` | Unevaluated `match.call()` for reproducibility | | `$spec$dist` | Baseline distribution label | | `$spec$control` | Control options passed at construction | | `$data$time` | Follow-up time vector | | `$data$status` | Event indicator (numeric) | | `$data$x` | Design matrix | | `$fit$theta` | Coefficient vector: estimates when `fit = TRUE`, the starting values when `fit = FALSE` | | `$fit$converged` | `TRUE`/`FALSE` from the optimizer; `NA` when `fit = FALSE` | | `$fit$objective` | Objective at convergence (the log-likelihood, except as noted below); `NA` when `fit = FALSE` | | `$legacy_args` | Named pass-through arguments for parity | > **Note:** `$fit$objective` above is a *number*: the value the optimizer > reached, which is the log-likelihood with the one exception described > below. `print()` and `summary()` label it > `log-lik` and `log_lik`, which are the names to prefer when writing about > it. The exception is a fit with `objective = "sas"` that read an > interval-censored row (one of positive weight, not dropped by a phase > design): there the value is `PROC HAZARD`'s objective, not a > log-likelihood, and `print()` labels it `SAS objective`, while > `summary()` returns `log_lik = NA` and the value as `objective_value`. The > field is distinct from the `objective` *argument* of `hazard()`, which > selects how it is calculated: `"likelihood"` (the default) or `"sas"`, > reproducing what `PROC HAZARD` accumulates for an interval-censored row. See `?hazard`; the > SAS objective exists to reproduce legacy runs, not for new analyses. --- ## Prediction In SAS HAZARD the `P` (predict / print) option on the `PROC HAZARD` line writes predicted survival, hazard, and cumulative-hazard tables to the output dataset. In R the same predictions come from `predict.hazard()` on the fitted object, with the requested quantity chosen via the `type=` argument. `predict.hazard()` currently supports four output types (`"linear_predictor"`, `"hazard"`, `"survival"`, and `"cumulative_hazard"`), covering every quantity the SAS `P` option produces: ```{r} #| eval: false #| label: predict-example # Linear predictor (X %*% theta) eta <- predict(fit, type = "linear_predictor") # Hazard scale (exp(eta)) hz <- predict(fit, type = "hazard") ``` The code chunk above shows only `"linear_predictor"` and `"hazard"`; `"survival"` and `"cumulative_hazard"` follow the same pattern. For multiphase fits pass `decompose = TRUE` to get per-phase cumulative hazard contributions instead of just the total. --- ## Automated translation with `hzr_translate_sas()` Everything above is a mapping you apply by hand. `hzr_translate_sas()` applies it for you: point it at a `.sas` file and it parses the `PROC HAZARD` / `PROC HAZPRED` blocks and emits a `.qmd` document whose chunks are the `hazard()` / `predict.hazard()` calls those blocks describe. It applies the statement-by-statement mapping this vignette documents, with differences of form you will see in its output: it writes the vector interface (`data =`, `time =`, `status =`) where this vignette writes a `Surv()` formula, it keeps the SAS job's upper-case column names, and it turns `EVENT` into a status and `weights` pair, covered below. The phases, `fixed = "m"`, `theta` and `control` it emits for the AVC model above are the ones the worked example uses, and on `avc` its call reaches the same estimates. It is useful for a one-off job, and essential for migrating a corpus of hundreds. ::: {.callout-warning} ## Experimental: what it will and will not translate The emitted document renders. The `hazard()` chunk binds its fit to a name and passes `fit = TRUE`, so the `predict()` chunks have something to predict from. What you should not read into that is "the translator works": it is a translation aid, and four limits are worth knowing before you point it at a corpus. - **A `SELECTION` statement is translated into a real screen, and the screen may choose a different model than `PROC HAZARD` did.** The job's candidates, per-variable flags (`/S`, `/I`, `/E`) and `SLE`/`SLS` thresholds are carried into an `hzr_stepwise()` call. `PROC HAZARD` uses approximate variances while selecting. In a screen that enters variables, the entry statistic here reproduces that approximation (except for a candidate refitted because its information is indefinite); in one that removes them, the Wald tests behind each removal use the full Hessian, so removals can differ. A `BACKWARD` screen only removes and a `NOSTEPWISE` one only enters, and the callout the emitted document carries above the chunk describes only what that job's screen does. Compare the selected model against the SAS listing before relying on it. Options with no faithful translation are still refused, each emitting a `stop()`: `FAST`, `MAXVARS`, `RESTRICT`, a per-variable `MOVE=` or `ORDER=`, and a variable held by `/I` in one phase but movable in another. On the reference corpus two jobs in four translate; the other two are a cross-phase `/I` and a `RESTRICT`. - **`LCENSOR` combined with `ICENSOR` is refused.** `hazard()`'s single `time_lower` argument carries the entry time for status 0/1 rows and the interval's lower bound for status 2 rows; one column cannot express both. Either statement alone translates. - **A job `PROC HAZARD` itself refuses, or one whose model this translation cannot emit, warns loudly and still fits.** A syntax error in a `PARMS` operand or on the `PROC` line (a value its lexer does not read as a number, a keyword outside the grammar, a template's `?`), a `SETG3` refusal, `FIXMNU1`, a non-zero `DELTA` on an active early phase, `FIXTAU` with no `TAU` written, and `FIXGE2`/`FIXGAE2` without `WEIBULL` each get a `warning()` chunk immediately above the fit, naming the cause, and a row in `$untranslated`. A fit standing in for a job `PROC HAZARD` does not run, or for a model it does not fit, is the answer this package most wants to avoid, so it is labelled where you cannot miss it rather than left to be inferred from the estimates. The warning is **per job**. A file holding several jobs emits one fit chunk each, and only the refused job's fit is preceded by a warning, so the document renders to completion and every other job runs unchanged. An earlier draft of this work stopped instead; Quarto then exited 1 and wrote no output document at all, including for the jobs *before* the refused one. Two limits on that. An **unknown** keyword is *recorded, not refused*: a `%repeat` call brings a DATA step's own keywords through the same parser, so refusing them would flag jobs that run. And an operand carrying a **macro reference** gets no syntax verdict at all, because SAS expands it before `PROC HAZARD` reads the statement, so the translation cannot know what arrived. It still warns on a job whose phases were built: the unread operand may be the one that sets a shape, and the emitted phase would otherwise carry SAS's default where the job wrote something else. One case still halts, and it halts more than the job. Where the value `PROC HAZARD` refuses is also one `hzr_phase()` will not build, the warning and the fit are both emitted as usual --- but `hzr_phase()` then rejects the value when the chunk actually runs, and an error in a chunk stops Quarto outright: it exits 1 and writes **no output document at all**, including for the jobs before this one. The console shows only `hzr_phase()`'s own message, because `knitr` collects warnings into the document that was never written. So the reason the job was refused does not reach you; only the symptom does. Delete that job from the file if you need the other jobs' results before you can correct it. - **Prediction grids the parser cannot resolve are refused whole**, and an unresolved `INHAZ=` stops the render on purpose. Both are covered under "What it doesn't translate" below. Note that resolving an external `INHAZ=` at all requires you to pass `librefs=`. - **Confidence limits from a loaded `INHAZ=` fit may stop the render.** A translated `PROC HAZPRED` block asks for `se.fit = TRUE` unless the SAS job says `NOCL`, and `predict()` on an `hzr_outhaz` object refuses `se.fit = TRUE` when `PROC HAZARD` estimated a late shape parameter on a composite scale (`log(GAMMA*ETA - 2)` and friends), which is the *generic* unconstrained three-phase case, not an exotic one. Point predictions are unaffected; drop `se.fit` from the emitted call to get them. Read `job$coverage` as a measure of how much of the SAS job the *parser recognised*, not of whether the result runs. A job can report full coverage with nothing in `$untranslated` and still error on render. The API and the emitted format will change as these limits close; the 1.2.2 entry of `NEWS.md` carries the issue numbers. ::: Take the `PROC HAZARD` block shipped as `system.file("extdata", "hz-example.sas", package = "TemporalHazard")` (de-identified from the public `hazard` examples): ```sas %HAZARD( PROC HAZARD DATA=AVCS CONSERVE OUTHAZ=EXAMPLES.HZDEATH QUASI CONDITION=14; EVENT DEAD; TIME INT_DEAD; PARMS MUE=0.2361727 THALF=0.1512095 NU=1.438652 M=1 FIXM MUC=0.0005436977; ); ``` ```{r} #| eval: false #| label: translate-example job <- hzr_translate_sas( system.file("extdata", "hz-example.sas", package = "TemporalHazard"), out_dir = "translated" ) ``` The emitted `.qmd` carries three chunks past its `setup`: a `data` chunk, a `status` chunk, and a `fit` chunk (reformatted here for width; the emitted document wraps at 60 columns like any `deparse()`d call, and the three are run together in one block below, each marked by a comment, rather than shown as the three separate chunks they are). `DATA` step translation is out of scope, so the `data` chunk is a guard, not a re-implementation. It checks that `AVCS` is already in scope and stops with a named error if not, rather than silently proceeding without data. The `status` chunk is where `EVENT DEAD` becomes something `hazard()` can take, and it is worth reading closely if you are checking a translation against the job it came from: ```{r} #| eval: false #| label: translate-output # chunk: data if (!exists("AVCS")) { stop("This job read ", "AVCS", " from a SAS DATA step, which ", "hzr_translate_sas() does not translate. Assign ", "AVCS", " before rendering.") } # chunk: status AVCS[[".hzr_status"]] <- with(AVCS, { if (any(is.na(DEAD) | DEAD < 0)) { stop("This job has rows where the EVENT count (DEAD) is missing or ", "negative. ...") # message trimmed here; see below } ifelse(DEAD > 0, 1, 0) }) # chunk: fit fit <- hazard( data = AVCS, time = INT_DEAD, status = .hzr_status, fit = TRUE, dist = "multiphase", phases = list( hzr_phase("cdf", t_half = 0.1512095, nu = 1.438652, m = 1, fixed = "m"), hzr_phase("constant") ), theta = c( log(0.2361727), log(0.1512095), 1.438652, 1, log(0.0005436977) ), weights = ifelse(DEAD > 0, DEAD, 1), control = list(conserve = TRUE) ) ``` `EVENT DEAD` names a **count**, not a flag. In the reference program a row with `DEAD = 2` is two events, and `setlik.c` weights its contribution accordingly. `hazard()` takes a 0/1 status and a separate `weights` argument, so the one SAS statement translates into two: the status derives from `DEAD > 0`, and the count carries into `weights`. Read `status = .hzr_status` and `weights = ifelse(DEAD > 0, DEAD, 1)` as a pair and they are your `EVENT` statement, not two new modelling choices. If you have translated this job before version 1.2.2, note that `DEAD = 2` used to map straight onto `status = 2` and be fitted as interval-censored, which is a different likelihood branch rather than an under-count. The guard in front of it refuses the rows SAS would have deleted. A missing or negative count is dropped outright by the reference program: `readc1.c` sets `mdel` on a missing count and `del` on a negative one, and `readobs.c` then skips the row and subtracts it from `Nobs`, so it contributes nothing to the likelihood. `hazard()` has no equivalent of that deletion, and translating such a row as censored would change the fit without saying so, so the chunk stops instead and asks you to drop the rows yourself. When it fires, subset to `!is.na(DEAD) & DEAD >= 0` and check your row count against the deletion tallies SAS printed for the job. The fit is bound to `fit` and `fit = TRUE` is emitted, which is what makes the document runnable: a later `predict()` chunk from a `PROC HAZPRED` block in the same job refers to that name. `CONSERVE` becomes `control$conserve = TRUE`, and the `PARMS` early-phase shape (`MUE`/`THALF`/`NU`/`M FIXM`) becomes a `hzr_phase("cdf", ...)` with `m` fixed. `FIXM` holds `m` only through `fixed = "m"` on the phase. `QUASI` and `CONDITION=14` are not in `control`, because `hazard()` reads neither: they are listed in `job$untranslated` with the reason. `QUASI` chooses `PROC HAZARD`'s optimizer, and `hazard()` has no choice to make: it always uses BFGS (a multiphase fit may run a Nelder-Mead warm-up first, and a stop that fails SAS's gradient test continues with `stats::nlm()`). `CONDITION=` stops `PROC HAZARD`'s optimizer as ill-conditioned once its Hessian approximation's condition estimate passes the limit, and `hazard()` has no such stop; it warns about the final Hessian instead. Earlier versions of the translator wrote both into `control`, where nothing read them. `dist = "multiphase"` is emitted whenever `phases` is, and `theta` carries the full interleaved multiphase start vector (early scale, then every early-phase shape parameter in phase order, then the constant/late scale), not just the phases' log(mu) starts. ### What it doesn't translate Constructs the parser can't resolve are never silently dropped: they land on the returned `hzr_sas_job` object and render as visible `UNTRANSLATED` callouts in the `.qmd` (where they do not stop the document outright, above), so a translated document shows its own gaps instead of quietly under-reporting the original job. Two gaps are common enough in production jobs to know about going in: - **A `PROC HAZPRED` grid is rebuilt from the job's own DATA steps, and some of those steps are only partly read.** The translator follows `DATA=` back the way SAS does: the last step that wrote the dataset before the `PROC HAZPRED`, then each `SET` input of that step, and so on. So a `DATA DESIGN` step's covariates, a `SET PREDICT DIGITAL` that stacks two grids, and a `YEARS = MONTHS/12` all reach the grid, and the emitted `predict()` runs at the variable the `TIME` statement names. Arithmetic, `LOG` and `EXP`, `DO` lists with constant bounds, `OUTPUT`, `DROP`, `KEEP` and `PROC SORT ... BY` are translated. A statement that sets a value the translator can't compute, such as an `IF ... THEN` assignment or a function it does not carry, is recorded as `UNTRANSLATED`, the grid chunk warns, and that variable is `NA` in the grid, so any prediction that uses it is `NA` too. The rest of the grid is still built. A step that decides which rows exist in a way it can't read (an `INFILE`, a `MERGE`, a subsetting `IF` or `IF ... THEN DELETE`, another procedure's `OUT=`, a `PROC SQL` `CREATE TABLE`, a loop bound taken from data) refuses the whole grid: the `predict()` chunks stop, so build the `newdata` grid by hand and pass it to `predict.hazard()`. The same happens when anything between a step the grid uses and the `PROC HAZPRED` could rewrite the grid out of sight: an `%INCLUDE`, a call of a macro the translator cannot read, open-code `%IF` or `%DO`, or `OPTIONS OBS=`. Such a call before the grid's first step keeps the grid and is recorded as `UNTRANSLATED`, since it may change session state the grid depends on. On the public corpus, 15 of 36 `PROC HAZPRED` blocks now get a grid, 7 of them with a warning; before 1.2.13, 13 did, several of them missing covariates or rows with nothing recorded. - **An unresolved `INHAZ=` fails the render on purpose.** If a `PROC HAZPRED` job's fitted-model dataset can't be found (neither from another translated job's `OUTHAZ=` nor from the `librefs=` argument you pass in), the document gets an `inhaz-unresolved` chunk ahead of the grid and `predict()` chunks whose whole body is a `stop()` naming the libref, so rendering fails loudly instead of producing predictions against a model it never loaded. ## Known limitations vs. SAS HAZARD A SAS HAZARD veteran migrating to `TemporalHazard` should be aware of the following scope limits. Detailed status per feature is tracked in `inst/dev/SAS-PARITY-GAP-ANALYSIS.md` and `inst/dev/DEVELOPMENT-PLAN.md`. ### Stepwise variable selection (`SELECTION` statement) - **Supported:** `hzr_stepwise()` implements forward / backward / two-way stepwise with SAS-style SLENTRY / SLSTAY thresholds, per-phase entry for multiphase, and a MOVE oscillation guard. - **Not yet supported:** FAST-screening (Lawless-Singhal approximate Wald updates). Much of what it exists to save is already gone, though: the default `criterion = "score"` tests each candidate at the current estimates without refitting it, except at a step where the information cannot be formed, where it refits and Wald-tests each candidate instead. `criterion = "wald"` refits each *entry* candidate, twice for a multiphase model, once from the default start and once from the current model's estimates, keeping the better. Removals never refit, under either criterion: a drop is tested on the current model's Wald p-value against `slstay`, which is what SAS does, and only the accepted drop is refitted. ### Per-phase time-varying windows - **Supported:** `time_windows` applies one piecewise-constant cut-point set across all covariates. - **Not yet supported:** distinct `EARLY`/`LATE` window sets per phase. Workaround: expand the design matrix manually before calling `hazard()`. ### Output datasets (`OUTEST=`, `OUTVCOV=`) - The R equivalents are `coef(fit)` and `vcov(fit)` in memory; there is no automatic dataset-export mode. `saveRDS()` them yourself if you need an on-disk artefact. ### Density / quantile prediction types - `predict.hazard()` covers `hazard`, `cumulative_hazard`, `survival`, and `linear_predictor`. Density and quantile (median survival) types are not wired up; derive them from `cumulative_hazard` and `survival` predictions if needed. ### Previously listed gaps that have since been closed For users migrating from older TemporalHazard versions or reading older SAS parity notes: - **`weights` on all distributions** (shipped v0.9.6). Weibull, exponential, log-logistic, log-normal, and multiphase all honour observation weights end-to-end (LL + analytic gradient + multiphase Conservation of Events). - **`Surv(start, stop, event)` with `start > 0`** (shipped v0.9.7). Counting-process rows contribute `H(stop) - H(start)` for Weibull and multiphase. The previous `hazard()` guard against non-zero starts is gone. - **Delta-method prediction confidence limits** (shipped v0.9.8). Use `predict(..., se.fit = TRUE, level = 0.95)` to get a data frame with `fit`, `se.fit`, `lower`, and `upper` per row. Closed-form Jacobian for Weibull and multiphase, `numDeriv::jacobian` fallback for exponential / log-logistic / log-normal. --- ## Rcpp acceleration note `TemporalHazard` is currently pure R. If profiling against large real datasets turns up bottlenecks in the likelihood kernel or the optimizer inner loop, those specific functions will be re-implemented with Rcpp. The public interface (`hazard()`, `predict.hazard()`, etc.) will not change. ## References {.unnumbered} Blackstone EH, Naftel DC, Turner ME Jr. The decomposition of time-varying hazard into phases, each incorporating a separate stream of concomitant information. *J Am Stat Assoc.* 1986;81(395):615–624. doi: [10.1080/01621459.1986.10478314](https://doi.org/10.1080/01621459.1986.10478314)