Migrating from SAS HAZARD to TemporalHazard

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:

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"
  )
)

Formal argument map: SAS HAZARD/C → hazard()

SAS Statement	Legacy Input	C Concept	R Parameter	Required	Expected Type	Transform Rule	Status	Notes
HAZARD	TIME variable	obs time array	time	TRUE	numeric vector	pass through as numeric	implemented	Core observation time input.
HAZARD	EVENT/censor variable	event indicator array	status	TRUE	numeric/logical vector	coerce to numeric 0/1	implemented	Event indicator currently retained as numeric in object\(data\)status.
HAZARD	X covariate block	design matrix	x	FALSE	numeric matrix or data.frame	data.frame -> data.matrix	implemented	Future versions will support richer design encoding helpers.
HAZARD	initial parameters	parameter vector	theta	FALSE	numeric vector	length must equal ncol(x) when x is present	implemented	Used by predict.hazard as coefficient vector.
HAZARD	baseline distribution	phase distribution selector	dist	FALSE	character scalar	normalized lower-case label	implemented	Current default is ‘weibull’; more options planned.
HAZARD	control options	optimizer/control struct	control	FALSE	named list	stored in spec$control	implemented	Control list is stored and reserved for optimizer parity.
HAZARD	additional legacy options	misc legacy switches	…	FALSE	named arguments	stored in legacy_args for parity	implemented	Supports legacy-style pass-through options during migration.
TIME	t	time vector	time	TRUE	numeric vector	pass through	implemented	Canonical SAS migration uses TIME= mapping.
EVENT	status	event vector	status	TRUE	numeric/logical vector	coerce to numeric	implemented	Canonical SAS migration uses EVENT= mapping.
PARMS	theta0	starting coef	theta	FALSE	numeric vector	map PARMS/INITIAL to theta	planned	SAS PARMS syntax parser not yet implemented.
DIST	dist	dist selector	dist	FALSE	character scalar	map DIST= to dist	implemented	SAS DIST keyword maps directly to dist.
HAZARD	phases (3-phase structure)	3-phase Early/Const/Late	phases (list of hzr_phase())	FALSE	list of hzr_phase objects	list(early=hzr_phase(‘cdf’,…), constant=hzr_phase(‘constant’), late=hzr_phase(‘g3’,…))	implemented	Use dist=‘multiphase’ with phases argument. N-phase generalization of legacy 3-phase model.
HAZARD	MU_1, MU_2, MU_3	per-phase scale factors	mu (via exp(log_mu) in theta)	FALSE	numeric (per-phase)	exp(alpha_j) in internal parameterization; estimated on log scale	implemented	Each phase has its own scale mu_j(x) = exp(alpha_j + x*beta_j). Starting value via hzr_phase().
G1	THALF / RHO (early)	early half-life	hzr_phase(t_half=)	FALSE	positive scalar	maps directly to hzr_phase(t_half=) starting value	implemented	Half-life: time at which G(t_half) = 0.5. Same concept as SAS RHO/THALF.
G1	NU (early)	early time exponent	hzr_phase(nu=)	FALSE	numeric scalar	maps directly to hzr_phase(nu=) starting value	implemented	Time exponent controlling rate dynamics. Same parameter name as SAS early NU.
G1	M (early)	early shape	hzr_phase(m=)	FALSE	numeric scalar	maps directly to hzr_phase(m=) starting value	implemented	Shape exponent controlling distributional form. Same parameter name as SAS early M.
G1	DELTA (early)	early time transform	(not implemented)	FALSE	numeric scalar	time transform B(t) = (exp(delta*t)-1)/delta; drops out at DELTA = 0, non-zero DELTA not implemented	planned	NOT IMPLEMENTED, not absorbed. The C DELTA controls B(t) = (exp(delta*t)-1)/delta, which enters rho, the time argument and the density Jacobian separately; R computes the delta=0 branch of each. A job with DELTA != 0 is refused or flagged, never fitted silently.
G2	G2 constant phase	constant hazard rate phase	hzr_phase(‘constant’)	FALSE	hzr_phase(‘constant’)	hzr_phase(‘constant’) with no shape parameters	implemented	Flat background rate. No shape parameters estimated. SAS G2 equivalent.
G3	TAU (late)	late G3 scale	hzr_phase(‘g3’, tau=)	FALSE	positive scalar	maps directly to hzr_phase(‘g3’, tau=) for late phase	implemented	Late-phase G3 scale parameter. Maps directly to hzr_phase(‘g3’, tau=).
G3	GAMMA (late)	late G3 time exponent	hzr_phase(‘g3’, gamma=)	FALSE	numeric scalar	maps directly to hzr_phase(‘g3’, gamma=) for late phase	implemented	Late-phase G3 time exponent. Maps directly to hzr_phase(‘g3’, gamma=).
G3	ALPHA (late)	late G3 shape	hzr_phase(‘g3’, alpha=)	FALSE	numeric scalar	maps directly to hzr_phase(‘g3’, alpha=) for late phase	implemented	Late-phase G3 shape parameter. alpha=0 gives exponential case. Maps directly to hzr_phase(‘g3’, alpha=).
G3	ETA (late)	late G3 outer exponent	hzr_phase(‘g3’, eta=)	FALSE	numeric scalar	maps directly to hzr_phase(‘g3’, eta=) for late phase	implemented	Late-phase G3 outer exponent. Maps directly to hzr_phase(‘g3’, eta=).

---

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.

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:

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.

TIME INT_DEAD;

Maps directly to time:

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.

EVENT DEAD;

Maps to status:

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.

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:

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).

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.

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.

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.

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)
)

km <- hzr_kaplan(time = avc$int_dead, status = avc$dead)
head(km, 4)
#> Kaplan-Meier estimate with logit confidence limits
#> Events: 14  | Time points: 4 
#> Survival range: 0.9541 to 0.9836 
#> RMST at last event: 0.016 
#> 
#>         time n_risk n_event n_censor survival std_err cl_lower cl_upper cumhaz
#>  0.001368954    305       5        0   0.9836  0.0073   0.9612   0.9932 0.0165
#>  0.002737907    300       2        0   0.9770  0.0086   0.9527   0.9890 0.0232
#>  0.008213721    298       2        0   0.9705  0.0097   0.9443   0.9846 0.0300
#>  0.016427440    296       5        0   0.9541  0.0120   0.9240   0.9726 0.0470
#>   hazard density   life
#>  12.0744 11.9752 0.0014
#>   4.8862  4.7901 0.0027
#>   1.2298  1.1975 0.0081
#>   2.0741  1.9959 0.0160

nel <- hzr_nelson(time = avc$int_dead, event = avc$dead)
head(nel, 4)
#> Nelson cumulative hazard estimate with lognormal CL
#> Events: 14  | Time points: 4 
#> 
#>         time n_risk n_event weight_sum cumhaz std_err cl_lower cl_upper  hazard
#>  0.001368954    305       5          5 0.0164  0.0033   0.0109   0.0237 11.9752
#>  0.002737907    300       2          2 0.0231  0.0047   0.0153   0.0335  4.8699
#>  0.008213721    298       2          2 0.0298  0.0057   0.0201   0.0425  1.2256
#>  0.016427440    296       5          5 0.0467  0.0067   0.0350   0.0610  2.0565
#>  cum_events
#>           5
#>           7
#>           9
#>          14

head(hzr_gof(fit), 4)
#> Goodness-of-fit: observed vs. expected events
#> Distribution: weibull  | n = 305 
#> 
#> Total observed events: 68 
#> Total expected events: 68 
#> Final residual (E - O): 0 
#> Conservation ratio (E/O): 1 
#> 
#> Use plot columns: time, km_surv, par_surv, cum_observed, cum_expected, residual

hzr_deciles(fit, time = 120, groups = 10)
#> Decile-of-risk calibration (risk grouped at time = 120 )
#> 305 subjects, all included.
#> 10 groups, 68 observed events, 68 expected
#> 
#>  group  n events expected observed_rate expected_rate chi_sq p_value
#>      1 31      0     1.37        0.0000        0.0441 1.3700  0.2420
#>      2 30      1     1.64        0.0333        0.0546 0.2480  0.6190
#>      3 31      1     2.53        0.0323        0.0816 0.9260  0.3360
#>      4 30      2     3.49        0.0667        0.1160 0.6350  0.4250
#>      5 31      7     4.58        0.2260        0.1480 1.2800  0.2570
#>      6 30     10     5.67        0.3330        0.1890 3.3100  0.0689
#>      7 31      8     7.59        0.2580        0.2450 0.0218  0.8830
#>      8 30     14     8.72        0.4670        0.2910 3.2000  0.0735
#>      9 31     10    13.50        0.3230        0.4350 0.8960  0.3440
#>     10 30     15    18.90        0.5000        0.6320 0.8230  0.3640
#>  mean_survival mean_cumhaz
#>          0.950      0.0441
#>          0.931      0.0546
#>          0.906      0.0816
#>          0.865      0.1160
#>          0.818      0.1480
#>          0.707      0.1890
#>          0.655      0.2450
#>          0.538      0.2910
#>          0.441      0.4350
#>          0.221      0.6320
#> 
#> Overall: chi-sq = 12.7 on 9 df, p = 0.176

hzr_calibrate(avc$age, avc$dead, groups = 10, link = "logit")
#> Variable calibration (logit link, 10 groups)
#> 
#>  group  n events    mean     min     max  prob link_value
#>      1 30     11   3.519   1.051   5.388 0.367     -0.547
#>      2 31     11   8.665   5.421  11.532 0.355     -0.598
#>      3 30     13  15.194  11.631  18.497 0.433     -0.268
#>      4 31     11  23.077  18.990  27.828 0.355     -0.598
#>      5 30      7  43.544  28.124  57.167 0.233     -1.190
#>      6 31      3  72.066  59.730  86.408 0.097     -2.234
#>      7 30      2 101.154  86.507 117.522 0.067     -2.639
#>      8 31      3 162.739 121.169 203.733 0.097     -2.234
#>      9 30      4 247.051 205.343 297.140 0.133     -1.872
#>     10 31      3 530.623 324.573 790.981 0.097     -2.234

set.seed(1)
hzr_bootstrap(fit, n_boot = 20)  # small for vignette build
#> Bootstrap inference for hazard model
#> Mode: fixed refit 
#> Replicates: 20 successful, 0 failed
#> 
#>  parameter  n pct    mean     sd     min    max ci_lower ci_upper
#>         mu 20 100  0.0000 0.0000  0.0000 0.0000   0.0000   0.0000
#>         nu 20 100  0.2443 0.0176  0.2176 0.2867   0.2196   0.2799
#>        age 20 100 -0.0028 0.0016 -0.0070 0.0002  -0.0061  -0.0003
#>     status 20 100  0.7053 0.2266  0.3261 1.0926   0.3525   1.0779
#>        mal 20 100  0.4601 0.3280 -0.2047 0.9677  -0.1301   0.9321
#>     com_iv 20 100  0.7384 0.3397  0.0757 1.3481   0.1426   1.3177

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)
#> Competing risks cumulative incidence
#> Event types: 3  | Time points: 4 
#> Final survival: 0.9941 
#> Final incid_1 : 0.0059 
#> Final incid_2 : 0 
#> Final incid_3 : 0 
#> 
#>     time n_risk n_event_1 n_event_2 n_event_3 n_censor   surv incid_1 incid_2
#>  0.00068   1533         1         0         0        0 0.9993  0.0007       0
#>  0.00137   1532         4         0         0        0 0.9967  0.0033       0
#>  0.00171   1528         1         0         0        0 0.9961  0.0039       0
#>  0.00205   1527         3         0         0        0 0.9941  0.0059       0
#>  incid_3 se_surv   se_1 se_2 se_3
#>        0  0.0007 0.0007    0    0
#>        0  0.0015 0.0015    0    0
#>        0  0.0016 0.0016    0    0
#>        0  0.0020 0.0020    0    0

---

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.

%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.

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)
#> Multiphase hazard model (2 phases)
#>   observations: 310 
#>   predictors:   0 
#>   dist:         multiphase 
#>   phase 1:      early - cdf (early risk)
#>   phase 2:      constant - constant (flat rate)
#>   engine:       native-r-m2 
#>   converged:    TRUE 
#>   gradient:     relative 0.000123 (SAS/C requires <= 6.06e-06; not met, nlm code 2)
#>   log-lik:      -160.408 
#>   Note: Hessian ill-conditioned (rcond = 4.41e-12); standard errors may be unreliable.
#>   Not done in this run: none
#>   evaluations: fn=17, gr=1
#>   message:      continued with nlm() for 2 iterations (code 2) 
#> 
#> Coefficients (internal scale):
#> 
#>   Phase: early (cdf)
#>                   estimate    std_error    z_stat      p_value
#>   log_mu     -1.0484675704 7.527000e-01 -1.392942 1.636373e-01
#>   log_t_half -1.6580626653 3.631463e-01 -4.565826 4.975321e-06
#>   nu          1.4374160000 1.770559e-01  8.118432 4.722471e-16
#>   m           1.0000000000           NA        NA           NA
#>   age        -0.0320577400 9.381015e-03 -3.417300 6.324560e-04
#>   com_iv      1.3366750000 4.173461e-01  3.202798 1.360995e-03
#>   mal         0.6872028000 2.671643e-01  2.572210 1.010514e-02
#>   opmos      -0.0196337700 4.367582e-03 -4.495341 6.945864e-06
#>   op_age      0.0002086688 5.904656e-05  3.533971 4.093654e-04
#>   status      0.5169533000 1.571565e-01  3.289417 1.003951e-03
#> 
#>   Phase: constant (constant)
#>              estimate std_error    z_stat      p_value
#>   log_mu   -14.638385 3.0330322 -4.826321 1.390786e-06
#>   inc_surg   1.375285 0.6222371  2.210226 2.708945e-02
#>   orifice    3.117650 0.8999963  3.464070 5.320686e-04
#>   status     1.054988 0.5030359  2.097242 3.597215e-02

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:

# 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.

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):

%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;
);

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:

# 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

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