finegray

Version 1.3.8 | 2026-10-26

Quick Start · Installation · Worked Examples · Options · Stata-Tools · Report a bug

finegray fits Fine-Gray subdistribution hazard models for one event type in the presence of competing events, then turns the fit into cumulative incidence for a covariate profile, a group, or each observation. It handles delayed entry, stratified baselines, piecewise time-varying effects, and sampling weights, working on stset data without expanding it.

Quick Start

After installing, model local relapse (failtype == 1) in Stata's hypoxia data, with distant relapse as the competing event, and draw the cumulative incidence at the covariate means with a 95% band:

webuse hypoxia, clear
gen byte status = failtype
stset dftime, failure(dfcens==1) id(stnum)
finegray ifp tumsize pelnode, compete(status) cause(1)
finegray_cif, ci ylabel(0(0.1)0.5) yscale(range(0 0.5))

_{Drawn by demo/demo_finegray.do, which also sets Stata's sj scheme with a white graph region and an unboxed legend, and adds title(), ytitle(), and xtitle("Analysis time (years)").}

Features

• Fits the Fine and Gray (1999) model by a Mata scan of the stset records; no pseudo-observations or expanded data are created.
• finegray_cif reports cumulative incidence at chosen horizons or as a curve, for a covariate profile (at()) or one curve per group (over()), with analytic influence-function or bootstrap intervals; finegray_predict does the same per observation.
• The default variance is a fixed-weight sandwich; nuisance adds the censoring-weight (psi) term of Fine and Gray (1999) eq. 7-8, which is the variance stcrreg reports; cluster() gives the Zhou et al. (2012) cluster-robust form, and every fit records its choice in e(vce_meat).
• Supports delayed entry with stratified entry and censoring weights (Zhang, Zhang and Fine 2011), a stratified baseline subdistribution hazard (bstrata(), Zhou et al. 2011), piecewise-constant effects (tvc()), pweight and fweight, and mi estimate.
• Factor variables and interactions post e(b) and e(V) the way stcox does, so test, margins, contrast, and estimates table address them directly.
• Covariates that vary within subject, and option combinations with no published or derived variance, stop with r(198) and a message giving the reason.

Requirements

• Stata 16 or later.
• Data declared with stset. Without id() each record is one subject, as stset itself and stcrreg treat it; id() is needed only when a subject contributes more than one record.
• No external software or additional Stata package is required.

Installation

Install the released package from the Stata-Tools repository:

capture ado uninstall finegray
net install finegray, from("https://raw.githubusercontent.com/tpcopeland/Stata-Tools/main/finegray") replace

For a local checkout, replace the from() directory with the path to its finegray folder:

net install finegray, from("/path/to/Stata-Tools/finegray") replace

Commands

Command	Purpose
finegray	Fit a Fine-Gray subdistribution hazard model.
finegray_predict	Generate linear predictors, cumulative incidence, Schoenfeld residuals, or the baseline cumulative subdistribution hazard.
finegray_cif	Evaluate cumulative incidence at profiles or time points, optionally with confidence intervals or bootstrap confidence limits.
finegray_phtest	Inspect proportional-hazards behavior using raw Schoenfeld residual-time correlations.

Methods, formulas, grounding and the refusal rationales live in a separate help topic, help finegray_methods. It documents no command; each option's help entry links into it.

How It Works

finegray estimates regression coefficients for the subdistribution hazard of the integer event value selected by cause() in the compete() variable. Competing events remain represented in the risk-set construction, while inverse-probability-of-censoring weights account for right censoring; the implementation uses a native Mata scan rather than expanding the data into pseudo-observations.

With a proportional subdistribution-hazards model, the exponentiated coefficient is a subdistribution hazard ratio (SHR), and the fitted baseline subdistribution hazard is combined with the linear predictor to obtain a cumulative-incidence function (CIF). The default display is exponentiated coefficients; use noshr for log-SHR coefficients.

The SHR is a ratio of subdistribution hazards, not a relative risk: exp(b) is not the ratio of cumulative incidences at any horizon, although an SHR above 1 does mean a higher CIF at every horizon, and the two ratios are close only while the cumulative incidence is small. Its risk set keeps subjects after a competing event, so a covariate can raise the cause-specific hazard of the cause of interest and still carry an SHR below 1 when it raises the competing hazard more. Like any hazard ratio it is non-collapsible. Report the CIF from finegray_cif for absolute effects and consider a cause-specific stcox model alongside (Latouche et al. 2013; Austin, Lee and Fine 2016; Austin and Fine 2017). Tied cause-event times use the Breslow approximation, as stcrreg does; there is no ties() option.

The usual workflow is to declare survival-time data, fit one model, use finegray_predict for row-level quantities, use finegray_cif for a covariate profile or a curve, and use finegray_phtest to explore time-varying effects. xb can score compatible new data. Point cif and basecshazard predictions can also use compatible data while the active fit's cached event-time structure or a posted e(basehaz) is available; those quantities need _t or timevar(). finegray_cif, finegray_phtest, and the ci, schoenfeld, and bootstrap() options of finegray_predict require the original, unchanged stset data.

For delayed entry, declare enter(time ...) in stset and fit with the package's Weight 1 risk-set construction. When entry or censoring depends on a discrete factor, name it in both strata() and truncstrata(); that gives the published stratified weight (zzf1_stratified). Naming a grouping in only one option, or different groupings in the two, gives the factorized extension (zzf1_factorized), which has no published derivation and is experimental. Inspect e(lt_weight), e(lt_vce), and the weight diagnostics before interpreting the result.

Three different things are called strata here, and they are three separate options. bstrata(varname) stratifies the baseline subdistribution hazard: the stratified model of Zhou, Latouche, Rocha and Fine (2011), one unconstrained baseline per level with the coefficient vector shared, and the only one of the three that means what stcox, strata() means. strata() stratifies the censoring distribution, truncstrata() the entry distribution. The axes are independent and compose; bstrata(v) strata(v) is the paper's regularly stratified regime. bstrata() is right-censoring only and is refused with delayed entry, and it requires bstratum(#) when asking finegray_cif for a curve.

Performance

There are three ways to fit a Fine-Gray model in Stata: finegray, Stata's built-in stcrreg, and stcrprep (Lambert, SSC) followed by a weighted stcox. stcrreg forms the weighted risk sets internally at every iteration; stcrprep expands the data once with time-dependent censoring weights and hands the expanded rows to stcox; finegray scans the original rows in Mata and expands nothing.

Simulated competing-risks data, three covariates, median of three timed runs after one untimed warm-up, on Linux x86-64 under Stata 17 with c(processors) = 16:

N	finegray	stcrreg	stcrprep + stcox	stcox alone	vs stcrreg	vs stcrprep
109 (hypoxia)	0.009s	0.045s	0.020s	0.006s	5.0x	2.2x
500	0.031s	1.279s	0.130s	0.034s	41.3x	4.2x
1,000	0.049s	3.160s	0.353s	0.126s	64.5x	7.2x
2,000	0.087s	11.835s	1.635s	0.520s	136.0x	18.8x
5,000	0.214s	76.496s	14.848s	5.896s	357.5x	69.4x
10,000	0.358s	334.140s	76.550s	35.374s	933.4x	213.8x

The stcrprep + stcox column is the full sequence needed for one Fine-Gray fit. stcox alone is what a second model on the same expanded weights costs, since stcrprep is designed to compute the weights once and reuse them. Even so, it is slower than a complete finegray fit from N=1,000 upward (35.374s vs 0.358s at N=10,000). The expansion is why: it turns 5,000 subjects into 487,458 weighted rows and 10,000 into 1,790,932, and every operation afterwards pays for them.

A tvc() fit runs the same scan once per interval, so it costs roughly J times a proportional fit and keeps the same shape in N. Measured on simulated competing-risks data with three covariates, one time-varying, boundaries at the cause-event quartiles, set processors 1:

N	proportional	tvc(), J=2	tvc(), J=4
2,000	0.057s	0.094s (1.65x)	0.199s (3.49x)
10,000	0.249s	0.435s (1.75x)	0.827s (3.32x)
50,000	1.554s	2.629s (1.69x)	4.970s (3.20x)

The ratio is flat across a 25-fold change in N, because the piecewise fit is still linear in the original rows and still expands nothing.

Absolute seconds are machine-dependent; the ratios are the portable quantity. All three commands fit the same model; the benchmark scripts print the maximum relative coefficient difference across the three fits alongside the timings. Methodology, the on-demand stcrprep install, and how to reproduce the table: demo/README.md.

Which standard error am I getting?

One table, and one machine-readable counterpart per row: e(vce_meat) for the coefficients, r(se_method) for a CIF interval.

Coefficients (e(b), e(V)):

You type	You get	e(vce_meat)
(default)	fixed-weight sandwich	fixed_weight
nuisance	eta+psi, Fine & Gray (1999) eq. 7–8	nuisance_adjusted
norobust	model-based inverse information	not_applicable
cluster()	cluster-robust, Zhou et al. (2012)	fixed_weight
bstrata()	per-stratum sandwich, summed	fixed_weight
bstrata() nuisance	Zhou et al. (2011) §4.1 Σ_r	nuisance_adjusted
tvc()	sandwich from each subject's joint score across intervals (cluster-summed under cluster())	fixed_weight
tvc() nuisance	piecewise eta+psi	nuisance_adjusted
delayed entry	fixed-weight sandwich	fixed_weight
delayed entry + nuisance	Zhang, Zhang & Fine (2011) Appendix B terms added (pooled weight)	nuisance_adjusted
delayed entry + nuisance + strata()/truncstrata()	refused, r(198)	—
[pweight=]	sandwich with meat Σ (w_i s_i)², N/(N−1); Wogu et al. (2021) eq. 3 score	fixed_weight
[pweight=] cluster()	the same, summed within cluster	fixed_weight
[pweight=] norobust	refused, r(198)	—
[fweight=]	replication semantics: meat Σ w_i s_i s_i', N = Σ w; norobust allowed	fixed_weight

stcrreg's variance is the nuisance row, not the default: its scores are eta+psi ([ST] stcrreg, Methods and formulas), and finegray, nuisance reproduces its e(V) to numerical precision on untied right-censored data. The default fixed-weight e(V) differs from stcrreg's by the psi terms.

How far apart the default and nuisance variances are, on average and on a single dataset, is reported from simulation in help finegray_methods. The fixed-weight variance omits the psi term, so it approximates Fine and Gray's variance. Use nuisance when you want the eta+psi variance and the design allows it.

CIF intervals (finegray_cif, finegray_predict):

You type	You get	r(se_method)
ci	analytic influence function	analytic
ci bootstrap(#)	resampled SD over replications	bootstrap

Analytic CIF intervals are available after supported converged fits, including tvc(), bstrata(), and weighted fits. Bootstrap CIF intervals are available after the same fits except fweights; expand frequency-weighted data to individual records and refit first. Both methods require the original estimation data and are unavailable after an mi fit. Bootstrap resampling is by subject, or by whole cluster after cluster(), and requires at least 25 requested and successful refits. Analytic intervals treat the weights as fixed in all cases (they do not propagate the uncertainty in Ĝ even after a nuisance fit), so only bootstrap() propagates weight re-estimation.

What is still refused, and why. Three cells on the delayed-entry branch, and the design-weight cells, all r(198), each for its own reason:

Cell	Reason
nuisance + delayed entry + weight strata	The pooled-weight term (Zhang, Zhang & Fine 2011, Appendix B) is implemented; for the stratified weight their own Appendix E ships the first part only, "treating the weight function known", which is the default sandwich.
bstrata() + delayed entry	This package does not implement a baseline-stratified delayed-entry extension. Its baseline-stratification grounding covers right censoring; Kim et al. (2020), Section 7, discusses future work for case-cohort left-truncated data when independent truncation may fail.
tvc() + delayed entry	No source for time-varying subdistribution coefficients under left truncation at all, and the delayed-entry branch is already this package's own extension.
weights + nuisance	The estimated-Ĝ term is not derived under design weights in this release (Wogu et al. 2021 §4 write it with a weighted at-risk count that differs from their §3 estimator).
weights + strata() / truncstrata()	Not yet checked against another implementation. For strata() such a check is possible: a strata() term in survival::finegray's formula does stratify the censoring Kaplan-Meier, and the same call takes weights. For truncstrata() no source weights the delayed-entry side at all.
weights + bstrata() / tvc()	Linear in the same per-subject terms, but each combination needs its own external check before it is allowed.
weights + delayed entry	The delayed-entry branch is this package's own extension of Zhang, Zhang & Fine (2011); no source weights it.
finegray_phtest after a weighted fit	The correlation summary has no weighted form in the corpus; predict, schoenfeld still returns the residuals.

Choosing a Workflow

Goal	Command or pattern	Main considerations
Fit a model	finegray x1 x2, compete(status) cause(1)	Declare stset first; id() only for multiple-record data.
Score observations	finegray_predict newvar, xb	Compatible new data can be scored for the linear predictor.
Estimate row-level CIFs	finegray_predict newvar, cif timevar(t)	Use timevar() to evaluate all rows at a common horizon.
Compare a covariate profile	finegray_cif, at(...) attime(...)	at() supplies a profile; default profile values are estimation-sample means, weighted by the fit's weights when the fit was weighted.
Draw a CIF curve	finegray_cif, timepoints(...)	Use ci, nograph, saving(), and twoway options as needed.
Explore proportional hazards	finegray_predict stub, schoenfeld; then finegray_phtest	The residual variables are optional for inspection; the diagnostic reports correlations, not an omnibus chi-squared test.
Account for delayed entry	stset ..., enter(time entry) plus finegray, strata(...) truncstrata(...)	Name the discrete factors that entry or censoring depends on in both options; check positivity and the posted weight diagnostics.

Worked Examples

The examples below are intended to be run separately from a clean Stata session. Most use Stata's webuse hypoxia example dataset; examples 3a and 12 use webuse hiv_si and webuse pneumonia, the datasets of [ST] stcrreg examples 4 and 5.

1. Fit a basic model and inspect the estimates

webuse hypoxia, clear
gen byte status = failtype
stset dftime, failure(dfcens == 1) id(stnum)
finegray ifp tumsize pelnode, compete(status) cause(1)
ereturn list

2. Generate linear predictors and common-horizon CIFs

webuse hypoxia, clear
gen byte status = failtype
stset dftime, failure(dfcens == 1) id(stnum)
finegray ifp tumsize pelnode, compete(status) cause(1)
finegray_predict double xb_hat, xb
gen double horizon5 = 5
finegray_predict double cif5, cif timevar(horizon5) ci level(95)
summarize xb_hat cif5 cif5_lci cif5_uci

3. Estimate a profile-specific CIF table

webuse hypoxia, clear
gen byte status = failtype
stset dftime, failure(dfcens == 1) id(stnum)
finegray ifp tumsize pelnode, compete(status) cause(1)
finegray_cif, at(pelnode=1 ifp=20 tumsize=5) attime(1 3 5 8) ci
return list

3a. Compare cumulative incidence between two exposure groups

hiv_si is the Amsterdam Cohort dataset of [ST] stcrreg example 4: appearance of the SI HIV phenotype as the event of interest, an AIDS diagnosis as the competing event, and one binary covariate.

webuse hiv_si, clear
gen byte any_event = status > 0
stset time, failure(any_event == 1) id(patnr)
finegray ccr5, compete(status) cause(2)
finegray_cif, over(ccr5) attime(2 5 10) ci
finegray_cif, over(ccr5) ci

over(varname) evaluates one curve per level of a model variable in a single call: a table per level with attime(), otherwise one graph with a legend entry per level and each band shaded in its own curve's color. Every curve is the same computation as the standalone at(ccr5=#) call, r(table) carries a sixth column over naming each row's level, and r(at) has one row per curve. The same option draws every fitted baseline stratum after a bstrata() fit (over(pelnode) after bstrata(pelnode)), and a factor variable in an interaction with the other covariate fixed by at():

webuse hypoxia, clear
gen byte status = failtype
stset dftime, failure(dfcens == 1) id(stnum)
finegray i.pelnode c.ifp i.pelnode#c.ifp tumsize, compete(status) cause(1)
finegray_cif, over(pelnode) at(ifp=20) ci

4. Save a confidence-banded CIF curve on a custom grid

webuse hypoxia, clear
gen byte status = failtype
stset dftime, failure(dfcens == 1) id(stnum)
finegray ifp tumsize pelnode, compete(status) cause(1)
tempfile cifout
finegray_cif, timepoints(1 2 3 4 5 6 7 8) ci nograph saving("`cifout'", replace)
use "`cifout'", clear
list time cif se lci uci, noobs

5. Use factor variables and inspect proportional hazards

webuse hypoxia, clear
gen byte status = failtype
stset dftime, failure(dfcens == 1) id(stnum)
finegray i.pelnode c.ifp##c.tumsize, compete(status) cause(1)
finegray_predict double schoenfeld, schoenfeld
finegray_phtest, time(log) detail

6. Fit with censoring strata and clustered robust inference

webuse hypoxia, clear
gen byte status = failtype
gen int site = ceil(_n / 10)
stset dftime, failure(dfcens == 1) id(stnum)
finegray ifp tumsize pelnode, compete(status) cause(1) strata(pelnode) cluster(site)
display "clusters = " e(N_clust)

6a. Fit with a stratified baseline subdistribution hazard

webuse hypoxia, clear
gen byte status = failtype
stset dftime, failure(dfcens == 1) id(stnum)
finegray ifp tumsize, compete(status) cause(1) bstrata(pelnode)
display "baseline strata = " e(k_bstrata)
finegray_cif, attime(1 5) bstratum(0) ci nograph
finegray_cif, attime(1 5) bstratum(1) ci nograph
finegray_predict double h0, basecshazard
tabstat h0, by(pelnode) stat(min max)

6b. Fit a piecewise-constant time-varying effect

webuse hypoxia, clear
gen byte status = failtype
stset dftime, failure(dfcens == 1) id(stnum)
finegray ifp tumsize pelnode, compete(status) cause(1) tvc(pelnode) tsplit(1)
test [tvc1]pelnode = [tvc2]pelnode
finegray_predict double xb_late, xb attime(2)
finegray_cif, at(ifp=20 tumsize=5 pelnode=0) attime(1 3 5) nograph

The coefficient table gains one equation per interval (main for the effects held proportional, then tvc1, tvc2, …), with the interval bounds and the cause-event count behind each one printed underneath. test [tvc1]x = [tvc2]x is the Wald test of whether the effect is in fact constant.

stcrreg fits the same two-interval model through a threshold time interaction, and the two parameterizations map onto each other: stcrreg's main coefficient is finegray's [tvc1], and main + tvc is [tvc2].

stset dftime, failure(status == 1) id(stnum)
stcrreg ifp tumsize pelnode, compete(status == 2) tvc(pelnode) texp(_t > 1) noshr

6c. Fit with sampling weights

webuse hypoxia, clear
gen byte status = failtype
gen double sw = cond(pelnode == 1, 2, 1)
stset dftime, failure(dfcens == 1) id(stnum)
finegray ifp tumsize [pweight = sw], compete(status) cause(1)
finegray_cif, attime(1 5) ci nograph
display "`e(wtype)' `e(wexp)'  sum of weights = " e(sum_w)

[pweight=] multiplies every subject's risk-set contribution, event term and Breslow increment by its weight, while estimating the censoring distribution without weights from the analysis sample, as in survival::finegray. Population interpretation requires sampling to preserve the relevant censoring distribution and independent censoring. This is not a general case-cohort or outcome-dependent sampling estimator: score weights do not correct a censoring KM distorted by selection. The variance is the fixed-weight sandwich with meat Σ (w_i s_i)², omitting censoring-estimation uncertainty, so norobust is refused; [fweight=] is replication throughout and allows it. Weights must be constant within id() and are for right-censored data only. A weight declared in stset is not inherited: unlike stcox, finegray takes its weight from the command line, and fitting on weighted stset data with no weight typed is r(198), not an unweighted fit. _n and _N are refused inside a weight expression, because they name no variable and so cannot enter e(datasignature); anything else the expression reads that is not a variable (a scalar, an e() value, a subscript such as w[1]) is caught the other way, because every post-estimation command checks the total of the rebuilt weight column against e(sum_w) and refuses with r(459) when it has moved.

7. Use bootstrap confidence intervals for a CIF profile

webuse hypoxia, clear
gen byte status = failtype
gen int site = ceil(_n / 10)
stset dftime, failure(dfcens == 1) id(stnum)
finegray ifp tumsize pelnode, compete(status) cause(1) cluster(site)
finegray_cif, attime(1 5 8) ci bootstrap(25) seed(24680)
return list

The minimum is 25 replications; use substantially more replications for substantive inference.

8. Fit a delayed-entry model and inspect the weight path

Left truncation needs a cohort in which entry is not a function of the outcome time. In the cohort below, entry depends on z1, censoring does not, and z1 is a model covariate. Naming z1 in both strata() and truncstrata() gives the published stratified weight (zzf1_stratified), which remains valid when censoring does not depend on z1.

clear
set seed 20260713
set obs 24000
gen byte z1 = runiform() < 0.5
gen double z2 = rnormal()
gen byte g4 = ceil(runiform() * 4)
gen double ez = exp(0.5*z1 - 0.5*z2)
gen double p1 = 1 - (1 - 0.5)^ez
gen byte cause = cond(runiform() < p1, 1, 2)
gen double v = runiform()
gen double event_time = -ln(1 - (1 - (1 - v*p1)^(1/ez))/0.5) if cause == 1
replace event_time = rexponential(1/(0.5*exp(0.5*z1 + 0.5*z2))) if cause == 2
gen double censor_time = min(rexponential(1/0.15), 6)
gen double entry_time = rexponential(1/cond(z1 == 1, 1.6, 0.5))
gen double time = min(event_time, censor_time)
gen byte status = cond(event_time <= censor_time, cause, 0)
drop if !(entry_time < time)
keep in 1/4000
gen long id = _n
gen byte any_event = status > 0
stset time, failure(any_event == 1) id(id) enter(time entry_time)

finegray z1 z2, compete(status) cause(1) strata(z1) truncstrata(z1)
display "weight method = " "`e(lt_weight)'"
display "smallest weight probability = " e(min_weight_prob)
display "largest entry weight = " e(max_lt_weight)

Pooling an entry distribution that is not in fact common changes the estimates. Dropping strata(z1) truncstrata(z1) on these data leaves the posted weight as zzf1_geskus and attenuates the log-SHR on z1 from 0.472 to 0.368, about 22%:

finegray z1 z2, compete(status) cause(1) noshr
display "weight method = " "`e(lt_weight)'"

strata() and truncstrata() are separate options. Naming a variable in only one of them, or different variables in the two, gives the factorized extension (zzf1_factorized). It has no published derivation and is experimental: finegray says so in a note at fit time, and it is best kept for sensitivity analyses.

finegray z1 z2, compete(status) cause(1) truncstrata(z1)
finegray z1 z2, compete(status) cause(1) strata(g4) truncstrata(z1)
finegray_cif, attime(1 3 5) ci nograph

9. Request the baseline hazard and baseline cumulative subdistribution hazard

webuse hypoxia, clear
gen byte status = failtype
stset dftime, failure(dfcens == 1) id(stnum)
finegray ifp tumsize pelnode, compete(status) cause(1) basehaz
matrix list e(basehaz)
finegray_predict double baseline_subhaz, basecshazard
summarize baseline_subhaz

10. Pool a fit over multiply imputed data

webuse hypoxia, clear
gen byte status = failtype
replace ifp = . in 1/12
stset dftime, failure(status==1) id(stnum)
sts generate na1 = na
stset dftime, failure(status==2) id(stnum)
sts generate na2 = na
stset, clear
mi set wide
mi register imputed ifp
mi register regular tumsize pelnode status dftime dfcens stnum na1 na2
mi impute regress ifp = tumsize pelnode i.status na1 na2, add(10) rseed(20260825)
mi stset dftime, failure(dfcens == 1) id(stnum)
mi estimate, cmdok eform("SHR"): finegray ifp tumsize pelnode, compete(status) cause(1)

The imputation model includes the outcome: an indicator for each event type (i.status) and each cause's Nelson-Aalen cumulative hazard at the subject's exit time (na1, na2); omitting them biases the imputed covariate's SHR towards the null.

cmdok is required because mi estimate's supported-command list is internal to Stata. eform("SHR") labels the pooled column on the scale the coefficients are reported on elsewhere; without it mi estimate prints unlabeled log-SHRs. Post-estimation on an mi fit is refused with r(301); use mi extract #, clear and refit to predict.

11. Bootstrap the coefficients

finegray_cif, bootstrap() resamples the CIF. Coefficient inference that also propagates weight-estimation uncertainty has to resample the whole estimation sequence, which means re-running stset on the resampled subject identifiers, which the bootstrap prefix cannot do on its own. Wrap the sequence:

webuse hypoxia, clear
gen byte status = failtype
stset dftime, failure(dfcens == 1) id(stnum)
program define fgboot, eclass
    version 16.0
    capture drop _st _d _t _t0
    quietly stset dftime, failure(dfcens == 1) id(newid)
    finegray ifp tumsize pelnode, compete(status) cause(1) noshr nolog
end
bootstrap _b, reps(200) seed(13579) cluster(stnum) idcluster(newid): fgboot

12. What happens with an internal time-varying covariate

pneumonia is the multiple-record dataset of [ST] stcrreg example 5: ICU patients, death in the ICU as the cause of interest, discharge as the competing event, and an exposure that switches from 0 to 1 mid-stay.

webuse pneumonia, clear
gen byte status = cond(died == 1, 1, cond(discharged == 1, 2, 0))
bysort id (ndays): gen byte outcome = status[_N]
gen byte any_event = status > 0
stset ndays, id(id) failure(any_event == 1)
finegray age pneumonia, compete(outcome) cause(1)

That fit stops with r(198): pneumonia varies within id(). stcrreg handles these data by carrying each subject's last exposure value into the risk sets after a competing event; finegray refuses instead, because once a subject has failed from the competing risk its covariate path has no meaning while the subdistribution risk set still retains it. The exposure at admission is subject-constant and is accepted. That is a baseline-exposure model, a different estimand from the time-updated coefficient stcrreg reports:

bysort id (ndays): gen byte pneu0 = pneumonia[1]
finegray age pneu0, compete(outcome) cause(1)

Demo

From the root of a Stata-Tools checkout, run the demonstration:

stata-mp -b do finegray/demo/demo_finegray.do

The script adds the local package directory to the session-only adopath, exercises fitting, prediction, diagnostics, CIF tables, grouped CIF curves on hiv_si, the internal time-varying covariate refusal on pneumonia, delayed entry with entry strata, string identifiers, CIF and coefficient bootstrap inference, a stratified baseline subdistribution hazard, a piecewise-constant time-varying effect, mi estimate, cmdok eform("SHR"):, graph export, and regenerates both figures below in Stata's sj scheme. It also writes demo/finegray_results.xlsx with the sibling tabtools package (Stata 17+): a Table 1 from table1_tc, a side-by-side table of four model specifications from collect: and regtab, and formatted CIF and proportional-hazards-diagnostic tables from puttab. Agreement with stcrreg on two datasets, the tvc()/texp() parameterization mapping, and the split-record reduction are recomputed and checked with assert.

The repository also includes demo/benchmark_finegray.do and demo/benchmark_large.do for optional timing comparisons.

The main demo regenerates this cumulative-incidence curve with a pointwise 95% confidence band:

It also regenerates the bstrata() comparison below: one unconstrained baseline subdistribution hazard per stratum, one shared coefficient vector, with the two curves evaluated at the estimation-sample covariate means, drawn in one call by finegray_cif, over(pelnode).

Command Reference

finegray

finegray varlist [if] [in] [pweight fweight], compete(varname) cause(integer) [censvalue(integer) noshr level(#) strata(varlist) truncstrata(varlist) bstrata(varname) tvc(varlist) tsplit(numlist) cluster(varname) norobust noadjust nolog basehaz nuisance iterate(integer) tolerance(#)]

finegray [, noshr level(#) coeflegend]

The second form replays the last finegray results without refitting; noshr and level() are honored on replay, and coeflegend redisplays the table with each coefficient's _b[] name in place of its statistics (replay only; not combinable with level() or noshr).

varlist accepts numeric factor-variable notation, including continuous terms, indicators, interactions, and expanded interactions. compete() identifies the event-type variable and cause() selects the integer event value whose subdistribution hazard is modeled; censvalue() defaults to 0. The stset declaration supplies analysis time, failure coding, subject identifier, and optional delayed-entry time.

The default variance is a finite-sample-adjusted fixed-weight sandwich variance; with delayed entry, both censoring and entry weights are treated as fixed. The default display is the SHR scale. norobust requests the model-based inverse-information variance; noadjust suppresses the finite-sample sandwich adjustment and is available only with the robust variance. cluster() changes the sandwich clustering and strata() stratifies the right-censoring model.

truncstrata() is for delayed-entry risk-set weighting and requires delayed entry. bstrata() fits a stratified baseline subdistribution hazard with shared coefficients (Zhou et al. 2011); it is right-censoring only and is refused with delayed entry (r(198)). nuisance adds the weight-estimation (psi) term to the sandwich variance; it composes freely with bstrata() and tvc() (Zhou et al. 2011 sec. 4.1 and the piecewise eta+psi form respectively, both posting e(vce_meat) = nuisance_adjusted), and under delayed entry it adds the Zhang, Zhang & Fine (2011) Appendix B terms for the pooled weight, and is refused only when combined with strata() or truncstrata() there. nuisance requires the sandwich and is refused with norobust. basehaz posts the optional baseline matrix needed for direct inspection. nolog, iterate(), and tolerance() control optimization output and convergence.

[pweight=exp] multiplies every subject's risk-set contribution, event term and Breslow increment by its weight; computation rescales the pweights over rows contributing to cause-event risk sets, so results are invariant to a common rescaling, while e(sum_w), e(ll), e(ll_0) and the iteration log stay on the scale of the supplied weights, and a smallest-to-largest contributing-weight ratio below the normal double range is refused with r(430). The censoring Kaplan-Meier is unweighted on the analysis sample (as in survival::finegray); population interpretation requires that sampling preserve the required censoring distribution and independent censoring. This is not a general case-cohort or outcome-dependent sampling estimator. The variance is the fixed-weight sandwich with meat Σ (w_i s_i)² (cluster-summed under cluster()), omitting censoring-estimation uncertainty; norobust is refused. [fweight=exp] has replication semantics throughout (an fweighted fit equals the fit of the expanded data) and allows norobust. Weights are for right-censored data and must be constant within id(); they do not compose with nuisance, strata(), truncstrata(), bstrata(), tvc() or delayed entry, and finegray_phtest is refused after a weighted fit. A weight carried by stset is not inherited and is refused with r(198); _n and _N are refused inside the weight expression. e(wtype), e(wexp) and e(sum_w) record the weight; post-estimation re-evaluates e(wexp), whose variables are in e(datasignaturevars) and therefore covered by e(datasignature), so a changed weight is refused with r(459) rather than answered under different weights.

finegray_predict

finegray_predict [type] newvar [if] [in], [xb cif schoenfeld basecshazard timevar(varname) ci level(#) bootstrap(#) seed(#) attime(#)]

An optional storage type such as double may precede newvar. The default is xb. cif generates cumulative incidence at each observation's _t, or at the numeric values in timevar(). ci adds newvar_lci and newvar_uci; confidence intervals require cif, and bootstrap intervals additionally require bootstrap() with at least 25 replications. seed() requires bootstrap().

schoenfeld generates a stub for the raw Schoenfeld residuals, followed by numbered variables for additional model terms. basecshazard generates the fitted baseline cumulative subdistribution hazard and cannot be combined with ci or bootstrap(). The prediction types are mutually exclusive.

After a tvc() fit, xb is evaluated at each row's own _t unless attime(#) names a single time; cif accumulates the baseline interval by interval; basecshazard is unchanged. schoenfeld is refused. Both the analytic ci and ci bootstrap(#) are available (ci bootstrap(#) is refused after an fweight fit); analytic intervals treat the weights as fixed.

finegray_cif

finegray_cif [, at(string) over(varname) attime(numlist) timepoints(numlist) bstratum(#) ci level(#) saving(filename[, replace]) bootstrap(#) seed(#) nograph plotopts(line_options) plot#opts(line_options) ciopts(area_options) ci#opts(area_options) twoway_options]

Use at() to define a covariate profile; unspecified covariates default to estimation-sample means: the weighted means, taken under the fit's own [pweight=]/[fweight=] column, when the fit was weighted, so that an fweight fit reports the same default curve as the fit of the expanded data. A factor indicator's default is its weighted sample proportion. over(varname) draws one curve per level of a model variable, or per baseline stratum after bstrata(), in one call; each curve is bit-identical to its standalone at()/bstratum() call, and the results are stacked with an over column. Use attime() for a table at requested horizons or timepoints() for a curve grid; the two options are mutually exclusive. With neither option, the command uses the distinct baseline event-time grid, thinned when necessary.

bstratum(#) names the baseline stratum and is required after a bstrata() fit with more than one level: once the baselines are free, a covariate profile alone no longer identifies a curve. ci requests pointwise confidence limits, nograph suppresses the graph, and saving() writes time, cif, se, lci, and uci to a Stata dataset. bootstrap() refits the model for at least 25 subject- or cluster-level replications; seed() requires bootstrap(). plotopts() and ciopts() style every CIF line and every confidence band, and plot#opts() and ci#opts() style curve # only, overriding the built-in lwidth(medthick) line and color(%30) lwidth(none) band (the marginsplot convention). Remaining options are passed to the underlying twoway graph. The default graph has no note (add one with note()) and labels the y axis with a leading zero, with decimals chosen from the tick step; a format() in your own ylabel() replaces that. In attime() mode, graph options are ignored with a note.

finegray_phtest

finegray_phtest [, time(function) detail]

The default time(rank) scale uses event-time ranks; time(log) uses log event time and time(identity) uses event time. The command reports term-level correlations between raw Schoenfeld residuals and the selected time scale, along with event counts. detail displays the first 20 event-level residual contributions. This is an exploratory diagnostic and does not post an omnibus chi-squared statistic or p-value.

Options

Estimation options

Option	Default and use
compete(varname)	Required event-type variable; it must agree with the stset failure coding.
cause(#)	Required integer event code to model; it must differ from censvalue().
censvalue(#)	0; value treated as right censoring in compete().
strata(varlist)	None; stratifies the right-censoring model.
truncstrata(varlist)	None; delayed-entry truncation strata for Weight 1 estimation.
bstrata(varname)	None; stratifies the baseline subdistribution hazard with a shared coefficient vector. Composes with nuisance and tvc(). Right censoring only; not available with delayed entry.
tvc(varlist)	None; covariates whose coefficient is piecewise constant in analysis time. Requires tsplit(). Composes with bstrata() and nuisance. Right censoring only; not available with delayed entry.
tsplit(numlist)	None; the J−1 interior interval boundaries for tvc(), strictly ascending and positive. Requires tvc(). Intervals are (lower, upper], so an event exactly on a boundary falls in the earlier interval.
cluster(varname)	None; cluster variable for the sandwich variance.
norobust	Off; use model-based rather than sandwich variance.
noadjust	Off; suppress the finite-sample sandwich adjustment.
noshr	Off; display log-SHR coefficients instead of exponentiated coefficients.
level(#)	c(level); confidence level on Stata's cilevel rule (10 to 99.99 inclusive, at most two decimals) for displayed intervals and postestimation.
nolog	Off; suppress the iteration log.
basehaz	Off; post the baseline cumulative subdistribution hazard in e(basehaz).
nuisance	Off; add the weight-estimation (psi) term to the sandwich variance. Composes with bstrata(), tvc() and delayed entry (pooled weight); refused with norobust, and with strata()/truncstrata() under delayed entry.
iterate(#)	200; maximum number of optimization iterations.
tolerance(#)	1e-8; convergence tolerance.

norobust, noadjust, and nuisance have compatibility restrictions documented in the command help. In particular, noadjust is meaningful only with the sandwich variance, and the model-based variance cannot be combined with cluster().

Prediction options

Option	Use
xb	Generate the linear predictor; this is the default prediction type. After a tvc() fit it is a function of time, evaluated at each row's own _t.
attime(#)	Evaluate xb at the single time # instead of each row's _t. Requires xb and a tvc() fit.
cif	Generate cumulative incidence at _t or at timevar().
schoenfeld	Generate raw Schoenfeld residual variables on the original estimation data.
basecshazard	Generate the baseline cumulative subdistribution hazard.
timevar(varname)	Numeric evaluation time for cif or basecshazard.
ci	Add lower- and upper-confidence-limit variables for cif.
level(#)	c(level); confidence level on Stata's cilevel rule (10 to 99.99, at most two decimals) for ci.
bootstrap(#)	At least 25 requested and successful refits for bootstrap CIF confidence limits; requires ci and is refused after an fweight fit.
seed(#)	Reproducible bootstrap seed; requires bootstrap().

CIF and diagnostic options

Option	Use
at(string)	Profile values for finegray_cif; unspecified covariates use estimation-sample means (weighted by the fit's weights after a weighted fit).
over(varname)	One curve per level of a model variable (at most 20 distinct values), or per fitted baseline stratum when varname is the bstrata() variable; may not be combined with at() on the same variable or with bstratum().
attime(numlist)	Fixed horizons for a returned CIF table.
timepoints(numlist)	Evaluation grid for a returned CIF curve.
bstratum(#)	Baseline stratum to evaluate; required after a bstrata() fit with more than one level, because a covariate profile alone no longer identifies a curve.
ci	Add pointwise confidence limits to the finegray_cif table and curve.
level(#)	c(level); confidence level on Stata's cilevel rule (10 to 99.99, at most two decimals) for finegray_cif; requires ci.
bootstrap(#)	At least 25 subject- or cluster-level refits for a bootstrap finegray_cif band; requires ci, and is refused after an fweight fit.
seed(#)	Reproducible bootstrap seed for finegray_cif; requires bootstrap().
saving(filename[, replace])	Save CIF curve data with time, cif, se, lci, and uci (plus over with over()).
nograph	Suppress the CIF graph.
plotopts(), plot#opts()	Line options for every CIF line, or for curve # only; override the built-in lwidth(medthick).
ciopts(), ci#opts()	Area options for every confidence band, or for band # only; override the built-in color(%30) lwidth(none); require ci.
twoway_options	Graph options passed to the CIF twoway graph.
time(function)	Time scale for finegray_phtest: rank, log, or identity; default rank.
detail	Show the first 20 rows of the raw Schoenfeld residual matrix for finegray_phtest.

For finegray_cif, attime() and timepoints() cannot be combined. ci uses pointwise influence-function limits; bootstrap() uses refitted subject- or cluster-level samples and can report requested, successful, and failed replications in r().

Stored Results

After finegray

Standard estimation results include e(b), e(V), e(sample), e(N), e(depvar), and e(properties). e(depvar) is _t, the stset analysis time, as it is after stcrreg; the event-type variable is e(compete). For a factor-variable fit, e(b) and e(V) are named with the terms you typed (1.pelnode, 1.pelnode#c.ifp) and carry each base level (0b.pelnode) as a zero coefficient with a zero row and column, as stcox posts them, so test, testparm, estimates table, margins, contrast, and pwcompare address them directly; the internal design columns, one per estimated coefficient, are listed in e(designvars). The command also posts these scalars (with e(N_clust) when clustering is used):

e(N_fail), e(N_compete), e(N_cens), e(sum_w) (only with weights), e(ll), e(ll_0), e(chi2), e(p), e(df_m), e(rank), e(N_clust), e(converged), e(N_delayed), e(N_G_trunc), e(k_bstrata), e(n_intervals), e(k_tvc), e(level), e(cause), e(censvalue), e(iterate), e(tolerance), e(N_weight_strata), e(min_weight_prob), e(max_lt_weight), e(N_prob_warn), e(N_weight_warn), e(N_lt_prehole), and e(wsig_n) (only with weights), and e(rowsig_n), the number of estimation observations behind e(rowsig).

Posted macros are e(cmd), e(cmdline), e(refitcmd), e(predict), e(compete), e(compete_values), e(designvars), e(entryvar), e(idvar), e(mi_data), e(postest), e(fvvarlist), e(fvsemantic), e(strata), e(truncstrata), e(bstrata), e(bstrata_noevent), e(bstrata_noevent_x), e(bstrata_level) and e(bstrata_level_x) (single-level bstrata() fits only), e(tvc), e(tsplit), e(tvc_covariates), e(tvc_pos), e(tsplit_nfail), e(clustvar), e(monotone) (only when a term's likelihood looks monotone), e(wtype) and e(wexp) (only with weights), e(lt_weight), e(lt_norm) (stratified delayed-entry fits only), e(lt_vce), e(bh_seq), e(bh_key), e(weight_warn_strata), e(wsig) (only with weights), e(wsigkeyvars) (only with weights; the estimation variables plus the stset id() and entry variables), e(vce), e(vcetype) (Robust, unset under norobust), e(vce_meat), e(vce_adjust), e(title), e(marginsok), e(datasignature), e(datasignaturevars), e(rowsigvars) (the signature variables plus the stset id() and entry variables, kept out of e(datasignaturevars) so compress or recast of an id is accepted), e(rowsig) (an order-invariant checksum of the complete original estimation records; finite and noncryptographic, not a proof of dataset equality), and, only when the fit has factor terms, e(fvsupport_vars) with numbered e(fvsupport1), e(fvsupport2), … holding every raw factor level observed at fit time.

The optional matrix e(basehaz) is posted only with basehaz: its columns are time and cumhazard for an unstratified or single-level baseline, and bstratum, time, cumhazard for multiple baseline strata, with one row per distinct cause-event time within each stratum. With basehaz and tvc(), e(basehaz_tvc) also preserves the cumulative mass within each interval; its columns are bstratum, interval, time, and cumhazard. e(N_clust) is posted when clustering is used; delayed-entry diagnostics are populated when delayed-entry weighting applies.

For TVC CIFs, a saved fit without interval masses needs the original estimation data for rebuilding, or a refit with basehaz. Estimates saved by finegray 1.3.7 or earlier must be refit for this purpose.

After finegray_cif

The returned matrix r(table) contains the CIF table or curve (with over(), one block per curve and a sixth column over), and r(at) contains the evaluated profile (one row per curve with over()). Scalars are r(level) and r(cause), plus r(bstratum) and macro r(bstrata) after a bstrata() fit, and macros r(over) and r(levels) plus the matrix r(levels_mat) with over(). r(levels) and the row names of r(at) are display spellings; r(levels_mat) and the over column of r(table) carry the same levels as machine doubles and are the forms to compare or feed back into at()/bstratum(). The macro r(vce_adjust) is always none: the finite-sample factor e(vce_adjust) reports applies to the coefficient variance only, never to the cumulative-incidence variance. The macro r(se_method) names how the standard errors in r(table) were computed: analytic for the influence-function method, bootstrap when bootstrap() supplied them. With bootstrap(), the command additionally posts r(bootstrap_requested), r(bootstrap_success), and r(bootstrap_failed); r(profile_vars) identifies the profile variables used.

After finegray_phtest

The command posts r(N_fail), r(time), r(residual_scale), and matrix r(phtest). The matrix contains term-level correlations and event counts; no omnibus chi-squared result is posted.

After finegray_predict

finegray_predict clears r() and does not post estimation results. It creates the requested variable, plus _lci and _uci variables for cif, ci, or numbered Schoenfeld-residual variables when the requested stub represents multiple model terms.

Assumptions and Limits

• stset must define the analysis time and failure indicator; id() is required only when a subject contributes more than one record, and without it each record is one subject and e(idvar) is empty. finegray rejects an empty analysis sample and inconsistent competing-event coding.
• Multiple records per subject are supported when intervals are contiguous and covariates, censoring strata, truncation strata, and cluster membership are constant within subject. The estimation sample is reduced to one subject-level record for the fitted model, while interval records remain available for relevant postestimation quantities.
• Time-varying covariates and by prefixes are not supported. A covariate that varies within id() is refused with r(198) rather than reduced to one of its values; worked example 12 shows the refusal on webuse pneumonia and the subject-constant baseline-exposure fit that is accepted instead. Factor variables and interactions are supported, but a scoring dataset must contain compatible values and factor levels.
• Time-varying effects are supported through tvc() with tsplit(): the coefficient on a fixed covariate may be piecewise constant in analysis time. This is a different thing from a time-varying covariate: the covariate value stays known for every subject at every time, including after a competing event, so the subdistribution risk set and the CIF remain well defined. Every interval must carry at least one cause-of-interest event or the fit is refused (r(459)), and the per-interval event counts are printed with the results because an interval resting on a handful of events can reach a monotone likelihood and still converge. tvc() composes with bstrata() and nuisance but is refused with delayed entry (r(198)); after such a fit finegray_predict, schoenfeld and finegray_phtest are refused, while CIF confidence intervals are available both analytically (ci, fixed-weight) and by bootstrap().
• The default variance treats estimated censoring weights as fixed and, with delayed entry, also treats entry weights as fixed; it uses a sandwich estimator. norobust requests model-based information-matrix variance and is not a general replacement when weights are estimated. Worked example 11 shows the wrapper that bootstraps the coefficients over the whole estimation sequence, stset included.
• nuisance adds the weight-estimation (psi) term to the sandwich variance: Fine & Gray (1999) eq. 7–8 without delayed entry, the Zhang, Zhang & Fine (2011) Appendix B terms under delayed entry with the pooled weight; it is refused with strata()/truncstrata() under delayed entry. The default delayed-entry variance is fixed-weight; e(lt_vce) says which was computed.
• Delayed entry uses the package's Weight 1 construction. It checks censoring/truncation positivity and observed grouping support; weights can exceed 1. strata() defines censoring groups and truncstrata() defines entry groups; observed combinations form joint weight strata, with at most 100 cells and at least 20 estimation-sample subjects per cell. The same grouping in both gives the published stratified weight; anything else gives the experimental factorized extension. Continuous covariate-dependent entry is not supported.
• Postestimation requires a converged fit. finegray_cif, finegray_phtest, and the ci, schoenfeld, and bootstrap() options of finegray_predict require the original, unchanged stset data. Point xb predictions work on compatible new data; point cif and basecshazard predictions can also do so while the active fit has its cached or posted baseline. Package-created factor-design columns may be dropped and rebuilt on demand, but retain _fg_entry for postestimation after a multiple-record fit; modifying a present _fg_* design column in place is rejected.
• finegray_phtest is a residual-correlation diagnostic rather than a formal omnibus test. Interpret it alongside the scientific model and the observed event-time support.
• Tied cause-event times are handled by the Breslow approximation, matching stcrreg and cmprsk::crr; there is no ties() option. With many ties Breslow is a less accurate approximation than Efron ([ST] stcox). In Cox regression it tends to bias coefficients toward zero, and the bias grows as ties become heavier (Hertz-Picciotto and Rockhill 1997). Heavily tied day- or year-scale times are where this matters; a finer time scale is the remedy.
• Restricting follow-up. The stset options exit(), origin(), and scale() are honoured through the _t, _t0, and _d they produce. exit(time #) administratively censors subjects whose event falls after the limit, so their _d becomes 0 while compete() still holds the event code, and finegray stops with compete() and stset failure indicator do not match. Recode the competing-event variable to 0 wherever _d == 0 before fitting. The same message appears when stset was declared on the cause of interest alone, as stcrreg requires: finegray expects stset on any event, with compete() naming the type.
• With many thin bstrata() strata (tens of subjects and a handful of cause events each) the sandwich standard error runs below the Monte Carlo standard deviation: 5 to 7 percent below, with 95 percent coverage 0.936 to 0.938, for ten strata of 30 to 50 subjects in simulation. Prefer larger strata or bootstrap the coefficients as in worked example 11.
• Point cif and basecshazard predictions from finegray_predict score whatever data are in memory while the fit's cached baseline is available, without checking that they are the estimation data; finegray_cif refuses changed data with r(459). Re-run finegray if the data have changed and the estimation sample is what you mean. The commands that recompute from the estimation data (finegray_cif, finegray_predict, cif ci, Schoenfeld residuals, finegray_phtest) check every original estimation record in full: a change, including a permutation of one variable's values across observations, is r(459), and estimates saved by finegray 1.3.7 or earlier (no e(rowsig)) must be refit for every form of finegray_cif (including a point CIF), finegray_predict, cif ci, finegray_predict, schoenfeld and finegray_phtest, all r(301). On such estimates finegray_predict, xb still works, and cif and basecshazard work only when the saved estimates hold e(basehaz) (fitted with basehaz); otherwise they are r(301) too. Prediction accepts every factor level observed at fit time, including levels not typed as terms (for example 2.grp); estimates without e(fvsupport_vars) keep the narrower typed-level rule and need a refit. attime() and timepoints() ranges that numlist would round are r(198); list those times one by one.

References

• Austin PC and Fine JP (2017). Practical recommendations for reporting Fine-Gray model analyses for competing risk data. Statistics in Medicine, 36(27), 4391–4400. doi:10.1002/sim.7501.
• Austin PC, Lee DS, and Fine JP (2016). Introduction to the analysis of survival data in the presence of competing risks. Circulation, 133(6), 601–609. doi:10.1161/CIRCULATIONAHA.115.017719.
• Latouche A, Allignol A, Beyersmann J, Labopin M, and Fine JP (2013). A competing risks analysis should report results on all cause-specific hazards and cumulative incidence functions. Journal of Clinical Epidemiology, 66(6), 648–653. doi:10.1016/j.jclinepi.2012.09.017.
• Hertz-Picciotto I and Rockhill B (1997). Validity and efficiency of approximation methods for tied survival times in Cox regression. Biometrics, 53(3), 1151–1156. doi:10.2307/2533573.
• Fine JP and Gray RJ (1999). A proportional hazards model for the subdistribution of a competing risk. Journal of the American Statistical Association, 94(446), 496–509. doi:10.1080/01621459.1999.10474144.
• Zhang X, Zhang M-J, and Fine J (2011). A proportional hazards regression model for the subdistribution with right-censored and left-truncated competing risks data. Statistics in Medicine, 30(16), 1933–1951. doi:10.1002/sim.4264.
• Geskus RB (2011). Cause-specific cumulative incidence estimation and the Fine and Gray model under both left truncation and right censoring. Biometrics, 67(1), 39–49. doi:10.1111/j.1541-0420.2010.01420.x.
• Bellach A, Kosorok MR, Gilbert PB, and Fine JP (2020). General regression model for the subdistribution of a competing risk under left-truncation and right-censoring. Biometrika, 107(4), 949–964. doi:10.1093/biomet/asaa034.
• Kawaguchi ES, Shen JI, Suchard MA, and Li G (2021). Scalable algorithms for large competing risks data. Journal of Computational and Graphical Statistics, 30(3), 685–693. doi:10.1080/10618600.2020.1841650.
• Li J, Scheike TH, and Zhang MJ (2015). Checking Fine and Gray subdistribution hazards model with cumulative sums of residuals. Lifetime Data Analysis, 21(2), 197–217. doi:10.1007/s10985-014-9313-9.
• Zhou B, Latouche A, Rocha V, and Fine J (2011). Competing risks regression for stratified data. Biometrics, 67(2), 661–670. doi:10.1111/j.1541-0420.2010.01493.x.
• Zhou B, Fine J, Latouche A, and Labopin M (2012). Competing risks regression for clustered data. Biostatistics, 13(3), 371–383. doi:10.1093/biostatistics/kxr032.
• Wogu AF, Zhao S, Nichols HB, and Cai J (2021). Proportional subdistribution hazards model for competing risks in case-cohort studies. American Journal of Applied Mathematics, 9(5), 165–185. doi:10.11648/j.ajam.20210905.12.

Testing

The test files and how to run them are described in qa/README.md.

Version History

• 1.3.8 (2026-10-19) Pending SSC revision: pweights rescaled to mean one internally, so results no longer depend on the weight scale; pooled delayed-entry fits outside the positivity domain refused; postestimation checks a checksum of the estimation data (estimates from 1.3.7 or earlier must be refit for CIF CIs, Schoenfeld residuals, and finegray_phtest); literal attime()/timepoints() horizons kept exact alongside ranges; prediction accepts all fit-time factor levels; weighted finegray_cif works after mata clear or estimates use; positive delayed-entry probabilities are retained; inactive extreme weights no longer distort risk arithmetic or convergence; postestimation uses the same weight variable it checks; TVC interval masses are preserved separately; tiny CIFs, finite hazard products and extreme cloglog endpoints use stable arithmetic. Methods help cites Austin, Latouche and Fine (2020) and Poguntke et al. (2018) for the time-varying-covariate limitation. Help files: stored results list e(sample) under Functions; finegray.sthlp gains a References section; the methods help reference list adds Woodroofe (1985) and corrects the Kim et al. (2020) entry; source comments cite the published equation numbers of Zhang, Zhang and Fine (2011). Refits refused for a tvc() variable not in the model, a tsplit() interval without cause events, too few clusters, or an undersized delayed-entry weight stratum now leave the previous fit usable; negative-pweight and noninteger-fweight refusals are judged on the estimation sample only, and the weight expression is now evaluated once, so the draw that is checked is the draw that is used (a weight expression calling a random-number function can therefore give different results for a given seed than 1.3.8 as first released); factor variables with names of 27 or more characters (26 with ##) are no longer refused; a monotone likelihood is flagged with a warning and e(monotone); finegray_cif support notes print stored event times exactly and flag horizons just below a float-stored event time; the bstratum() refusal lists levels in a form it accepts back; predict, cif after estimates use without basehaz names estimates esample: as the repair.
• 1.3.7 (2026-09-29): e(vcetype) posted; censoring-floor note counts only weights that read a floored G(t-); factor terms with a base level (i.grp#c.x, i.a#i.b) now estimated as stcox does; margins on interaction-only fits fixed; labels containing $ or backticks carried verbatim; finegray_cif gains plotopts()/ciopts() styling, exact event-time grids, and cleaner over() legends; tsplit() boundaries kept exact with ranges.
• 1.3.6 (2026-09-14): tsplit() boundaries kept at full precision (were rounded to 9 significant digits); methods help no longer equates the default variance with stcrreg's (that is nuisance).

Earlier versions

• 1.3.5 (2026-09-13): Delayed-entry weights fixed for event/censoring ties and observation gaps; e(refitcmd) stores the fit-time factor expansion; bootstrap and string-literal weight fixes.
• 1.3.4 (2026-09-13): Stratified delayed-entry normalizer estimated on the identifiable region; e(N_lt_prehole) added.
• 1.3.3 (2026-09-12): Stratified delayed-entry weight gains the per-stratum normalizer; nuisance psi term under strata() corrected; finegray_cif no longer extends a flat tail past the timepoints() grid.
• 1.3.2 (2026-09-08): Current SSC Release. stset without id() accepted.
• 1.3.1 (2026-09-08): pweight scope corrected; postestimation fixes for restored estimates and multiple-record fits.
• 1.3.0 (2026-09-04): mi estimate support, bstrata(), tvc()/tsplit(), finegray_cif, over(), pweight/fweight, margins support. Breaking: e(covariates) renamed e(designvars).
• 1.2.0 (2026-08-16): Delayed-entry weights, nuisance-adjusted sandwich, baseline-hazard output, expanded CIF and diagnostic workflows.
• 1.1.0 (2026-07-10): CIF curves, multiple-record support, stratified censoring, postestimation confidence intervals.
• 1.0.0 (2026-04-06): First SSC Release.

Author

Timothy P Copeland, Karolinska Institutet

Questions and bug reports: GitHub issues. Please include the first line of which finegray and your Stata version.

To cite: Copeland TP. finegray: Fast Fine-Gray competing-risks regression. Stata-Tools, https://github.com/tpcopeland/Stata-Tools/tree/main/finegray. Please give the version you used (which finegray prints it).

License

MIT (see LICENSE)