Introduction

Many phenotypes in biobanks are naturally time-to-event outcomes: age at diagnosis, time from enrollment to a first event, or time to death. Analyzing such a phenotype as a binary case/control indicator throws away the follow-up time and the censoring pattern, and it loses power when the event is rare. SAIGEgds therefore supports time-to-event outcomes directly, via trait.type="survival" in seqFitNullGLMM_SPA() together with the new event.time argument.

The implementation follows the GATE method [2]. The model is a Cox proportional-hazards frailty model:

\[\lambda_i(t) \;=\; \lambda_0(t)\,\exp(X_i\beta + b_i), \qquad b \sim N(0,\ \tau_G \cdot \mathrm{GRM})\]

where $\lambda_0(t)$ is an unspecified baseline hazard, $X_i$ the fixed-effect covariates and $b$ a random effect (“frailty”) with covariance proportional to the genetic relationship matrix (GRM), which accounts for sample relatedness and population structure exactly as in the binary and quantitative traits.

Two ingredients make this scalable:

Cox fitted as a Poisson model. With the Breslow estimator of the baseline cumulative hazard $\Lambda_0$, the Cox score equations are identical to those of a Poisson model with mean $\mu_i = \Lambda_0(t_i)\exp(\eta_i)$ and “count” $y_i = $ event status. SAIGEgds exploits this equivalence and reuses the same IRLS + AI-REML machinery (with PCG and the same sparse genotype storage) that is used for the other trait types, so no risk-set expansion of the data is needed.
A Poisson saddlepoint approximation (SPA). For rare variants combined with a low event rate, the score statistic is far from normal and the usual normal approximation gives anticonservative p-values. The single-variant score test therefore switches to a saddlepoint approximation of the Poisson score whenever the normal p-value is below the threshold (0.0455 by default), equivalent to abs(z) > 2, where z is the standardized score statistic under the null hypothesis.

This vignette shows the workflow: fitting the null model, single-variant tests and set-based tests. It assumes familiarity with the main vignette (“SAIGEgds Tutorial (single variant tests)”), which covers GDS files, LD pruning and the GRM in more detail.

Input data requirements

For trait.type="survival":

The response in the formula is the 0/1 event status (1 = event, 0 = censored), not the time. event.time names a separate column of data holding the follow-up time (time to event or to censoring). It must be numeric, non-negative and non-missing.
A GRM is optional, supplied either as genotypes (gdsfile) or as a user-defined dense/sparse GRM (grm.mat). If neither is given, the frailty term is dropped and a standard Cox proportional-hazards model is fitted, assuming independent samples; the Poisson saddlepoint approximation is still applied in the association tests.
At least one covariate besides the intercept is needed. The intercept is not identifiable in a Cox model – it is absorbed by the baseline hazard – so it is removed from the design matrix. For the same reason, X.transform and use.offset are forced to FALSE for survival, whatever the user passes.
Subjects censored before the first observed event time are dropped before fitting: they contribute no risk-set information ($\Lambda_0 = \mu = 0$). The number of removed subjects is reported when verbose=TRUE.
Ties in the event times are handled by the Breslow approximation.

Example

library(SeqArray)
library(SAIGEgds)

# the example GDS file: 1,000 samples and 10,000 SNPs
(fn <- system.file("extdata", "grm1k_10k_snp.gds", package="SAIGEgds"))

## [1] "/tmp/Rtmp97fK15/Rinstfad3647542c/SAIGEgds/extdata/grm1k_10k_snp.gds"

gdsfile <- seqOpen(fn)

The package ships a small phenotype table; here a time-to-event outcome is simulated from an exponential model with an independent exponential censoring time, so that the outcome depends on the covariates x1 and x2 but not on any SNP (i.e. the global null hypothesis holds).

phenofn <- system.file("extdata", "pheno.txt.gz", package="SAIGEgds")
pheno <- read.table(phenofn, header=TRUE, as.is=TRUE)

set.seed(1000)
n <- nrow(pheno)
ftime <- rexp(n, rate=exp(0.3*pheno$x1 - 0.2*pheno$x2) * 0.05)  # time to event
ctime <- rexp(n, rate=0.03)                                     # censoring time

pheno$status <- as.integer(ftime <= ctime)  # 0/1 event status, the response
pheno$atime  <- pmin(ftime, ctime)          # observed follow-up time

head(pheno)

##   sample.id y     yy      x1 x2 status     atime
## 1        s1 0 4.5542  1.5118  1      0  4.982329
## 2        s2 0 3.7941  0.3898  1      1 11.250576
## 3        s3 0 5.0411 -0.6212  1      0 11.814053
## 4        s4 0 5.6394 -2.2147  1      0 22.077174
## 5        s5 0 4.2134  1.1249  1      1  8.322178
## 6        s6 0 4.6145 -0.0449  1      1  4.136543

table(pheno$status)  # 1: event, 0: censored

## 
##   0   1 
## 407 593

Fitting the null model

glmm <- seqFitNullGLMM_SPA(status ~ x1 + x2, pheno, gdsfile,
    trait.type="survival", event.time="atime", sample.col="sample.id")

## SAIGE association analysis:
## 2026-09-24 19:32:42
##     survival: removing 1 subject censored before the first event time
## Filtering variants:
## [..................................................]  0%, ETC: --- (1/1)    [==================================================] 100%, used 0s (1/1)    [==================================================] 100%, complete, 0s    
## MAC category for estimating variance ratio:
##     MAC[20, Inf):    30+ randomly from 9,649 variants
## Fit the null model: status ~ x1 + x2 + var(GRM)
##     # of samples: 999
##     # of variants in GRM: 9,649
##     MAF threshold for GRM: >= 0.01
##     using 1 thread
## Use dense genetic relationship matrix in the file:
##     /tmp/Rtmp97fK15/Rinstfad3647542c/SAIGEgds/extdata/grm1k_10k_snp.gds
## Loading SNP genotypes from the GDS file:
## [..................................................]  0%, ETC: ---    [==================================================] 100%, used 0s    
##     using 2.6M (stored in a sparse form)
## Survival outcome (event status): status
##     # of events: 593 (59.36%), # censored: 406
## Initial fixed-effect coefficients:
##     (Intercept)        x1         x2
##       0.5391031 0.2740208 -0.2859998
## Initial variance component estimates, tau (Sigma_E, Sigma_G):
##     tau: (1, 0.1)
##     fixed coeff: (0.2962767, -0.2238389)
## Iteration 1:
##     tau: (1, 0.0999238)
##     fixed coeff: (0.2962722, -0.2238348)
## Iteration 2:
##     tau: (1, 0.08603075)
##     fixed coeff: (0.2954319, -0.2231507)
## Iteration 3:
##     tau: (1, 0.07898628)
##     fixed coeff: (0.2949885, -0.2227349)
## Iteration 4:
##     tau: (1, 0.07531609)
##     fixed coeff: (0.2947523, -0.2225135)
## Final tau: (1, 0.07338431)
##     fixed coeff: (0.2946266, -0.2223867)
## 2026-09-24 19:32:46
## Calculate the average ratio of variances:
##     1, maf: 0.10911, mac: 218,   ratio: 0.8835 (var1: 0.486, var2: 0.55)
##     2, maf: 0.02653, mac: 53,    ratio: 0.8986 (var1: 0.469, var2: 0.522)
##     3, maf: 0.04955, mac: 99,    ratio: 0.8923 (var1: 0.48, var2: 0.538)
##     4, maf: 0.02352, mac: 47,    ratio: 0.8923 (var1: 0.585, var2: 0.656)
##     5, maf: 0.11411, mac: 228,   ratio: 0.8988 (var1: 0.438, var2: 0.487)
##     .........................
##     ratio avg: 0.8917671, sd: 0.01297065, CV: 0.0004848296
## 2026-09-24 19:32:47
## Done.

The log reports the number of events and censored subjects, the AI-REML iterations for the frailty variance, and the variance ratio used to approximate the score variance for the association tests.

glmm$coefficients   # log hazard ratios for x1 and x2 (no intercept)

##         x1         x2 
##  0.2946266 -0.2223867

glmm$tau            # (Sigma_E, Sigma_G); Sigma_G is the frailty variance

##    Sigma_E    Sigma_G 
## 1.00000000 0.07338431

glmm$converged

## [1] TRUE

Note what the components of a survival null model mean:

Components of a survival null model.
Component	Meaning for `trait.type="survival"`
`coefficients`	fixed effects on the log hazard scale; no intercept
`tau`	`Sigma_E` is fixed at 1 (Poisson scale), `Sigma_G` is the frailty variance on the log-hazard scale
`fitted.values`	$\mu_i = \hat\Lambda_0(t_i)\exp(\hat\eta_i)$, the expected number of events
`residuals`	$y_i - \mu_i$, the martingale residuals (summing to zero)
`var.ratio`	variance ratio(s) for the score test, computed with the Cox risk-set information correction

summary(glmm$fitted.values)   # expected numbers of events

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0006554 0.1717397 0.4377340 0.5935936 0.8744162 3.4909630

sum(glmm$residuals)           # martingale residuals sum to ~0

## [1] -3.556322e-13

head(glmm$var.ratio)

##     id        maf mac      var1      var2     ratio
## 1 2292 0.10910911 218 0.4858737 0.5499333 0.8835138
## 2 5689 0.02652653  53 0.4688258 0.5217384 0.8985841
## 3 6088 0.04954955  99 0.4797865 0.5377030 0.8922890
## 4 6067 0.02352352  47 0.5854085 0.6560404 0.8923361
## 5 6558 0.11411411 228 0.4380347 0.4873452 0.8988181
## 6 1612 0.02952953  59 0.5508638 0.6292359 0.8754489

Sanity check against `coxph()`

Without a GRM (neither gdsfile nor grm.mat), seqFitNullGLMM_SPA() fits a standard Cox proportional-hazards model without random effects, and its coefficients reproduce coxph() with Breslow ties. With the GRM the fixed effects are close but not identical: the frailty model estimates conditional (subject-specific) log hazard ratios, which are typically slightly larger in magnitude than the marginal estimates when the frailty variance is positive.

# the same model without a GRM (no random effects)
glmm0 <- seqFitNullGLMM_SPA(status ~ x1 + x2, pheno, trait.type="survival",
    event.time="atime", verbose=FALSE)

cf <- rbind(SAIGEgds.GRM=glmm$coefficients, SAIGEgds.noGRM=glmm0$coefficients)
if (requireNamespace("survival", quietly=TRUE))
{
    cx <- survival::coxph(survival::Surv(atime, status) ~ x1 + x2, pheno,
        ties="breslow")
    cf <- rbind(coxph=coef(cx), cf)
}
cf

##                       x1         x2
## coxph          0.2888592 -0.2167932
## SAIGEgds.GRM   0.2946266 -0.2223867
## SAIGEgds.noGRM 0.2888584 -0.2168195

Single-variant association tests

seqAssocGLMM_SPA() is called exactly as for a binary trait; it detects the trait type from the model object and switches to the Poisson score test with the saddlepoint approximation.

assoc <- seqAssocGLMM_SPA(gdsfile, glmm, mac=10)

## SAIGE association analysis:
## 2026-09-24 19:32:47
##     trait type: survival
##     # of samples: 999
##     # of variants: 10,000
##     MAF threshold: no
##     MAC threshold: >= 10
##     missing proportion threshold: <= 0.05
##     variance ratio for approximation: 0.8917671
##     genetic model: additive
##     # of processes: 1
## [..................................................]  0%, ETC: --- (1/1)    [==================================================] 100%, used 0s (1/1)    [==================================================] 100%, complete, 0s    
## # of variants after filtering by MAF, MAC and missing thresholds: 9,976
## P-value:
##      [0,5e-10]  (5e-10,5e-08]  (5e-08,5e-06] (5e-06,0.0005]     (0.0005,1] 
##              0              0              0              5           9971 
## 2026-09-24 19:32:47
## Done.

head(assoc)

##   id chr pos rs.id ref alt     AF.alt mac num        beta         SE      pval
## 1  1   1   1   rs1   1   2 0.03053053  61 999 -0.03274893 0.16779026 0.8452538
## 2  2   1   2   rs2   1   2 0.03803804  76 999 -0.18799559 0.14459258 0.1935412
## 3  3   1   3   rs3   1   2 0.02152152  43 999  0.23081107 0.23566363 0.3273780
## 4  4   1   4   rs4   1   2 0.38938939 778 999  0.07259130 0.06347674 0.2527941
## 5  5   1   5   rs5   1   2 0.03903904  78 999  0.11642187 0.16023682 0.4674948
## 6  6   1   6   rs6   1   2 0.05255255 105 999  0.08627289 0.14191033 0.5432276
##   method    p.norm converged
## 1 Normal 0.8452538      TRUE
## 2 Normal 0.1935412      TRUE
## 3 Normal 0.3273780      TRUE
## 4 Normal 0.2527941      TRUE
## 5 Normal 0.4674948      TRUE
## 6 Normal 0.5432276      TRUE

The output columns specific to this test:

beta, SE: the estimated log hazard ratio per copy of the alternative allele and its standard error;
pval: the reported p-value, from the Poisson SPA when it was applied;
method: Normal if the normal approximation was used, SPA if the saddlepoint approximation was applied (the ER method is used for binary traits only);
p.norm: the normal-approximation p-value, always kept for reference;
converged: whether the saddlepoint equation converged; if FALSE, treat the p-value with caution.

table(assoc$method)

## 
## Normal    SPA     ER 
##   9479    497      0

Kaplan-Meier curve for the top variant

The following code selects the variant with the smallest finite p-value and plots the unadjusted Kaplan-Meier estimate by alternative-allele dosage. The sample order is matched to the fitted null model before joining the genotypes to the phenotype data.

top <- assoc[which.min(replace(assoc$pval, !is.finite(assoc$pval), Inf)), ]
top[c("id", "chr", "pos", "ref", "alt", "AF.alt", "pval")]

##      id chr pos ref alt    AF.alt         pval
## 227 227   1 227   1   2 0.1826827 3.449984e-05

seqFilterPush(gdsfile)
seqSetFilter(gdsfile, sample.id=glmm$sample.id, variant.id=top$id,
    verbose=FALSE)
sample.id <- seqGetData(gdsfile, "sample.id")
dosage <- drop(seqGetData(gdsfile, "$dosage_alt2"))
seqFilterPop(gdsfile)

ii <- match(sample.id, pheno$sample.id)
plot.data <- data.frame(
    time = pheno$atime[ii],
    status = pheno$status[ii],
    genotype = factor(dosage))
plot.data <- plot.data[is.finite(dosage), ]
levels(plot.data$genotype) <- paste0("ALT dosage = ", levels(plot.data$genotype))

km <- survival::survfit(survival::Surv(time, status) ~ genotype, data=plot.data)
cols <- seq_along(km$strata)
plot(km, col=cols, lwd=2, mark.time=TRUE,
    xlab="Follow-up time", ylab="Survival probability")
legend("bottomleft", sub("^genotype=", "", names(km$strata)),
    col=cols, lwd=2, bty="n")

Because the variant was selected using the same association results, this curve is descriptive and should not be interpreted as an independent significance test. It is also unadjusted for the covariates and relatedness accounted for by the model.

Since no SNP affects the simulated outcome, the p-values should be uniform; the genomic inflation factor is close to 1:

p <- assoc$pval[is.finite(assoc$pval)]
median(qchisq(p, 1, lower.tail=FALSE)) / qchisq(0.5, 1)   # lambda_GC

## [1] 0.9809162

# QQ plot
plot(-log10(ppoints(length(p))), -log10(sort(p)), pch=20, cex=0.6,
    xlab=expression(Expected~~-log[10](italic(p))),
    ylab=expression(Observed~~-log[10](italic(p))))
abline(0, 1, col="red")

Set-based (aggregate) tests

Burden and ACAT-V tests support survival outcomes: a burden test is a single score test on the collapsed genotype, and ACAT-V is a Cauchy combination of single-variant p-values, so both inherit the Cox-via-Poisson score test.

units <- seqUnitSlidingWindows(gdsfile, win.size=500, win.shift=250)

## Chromosome 1, # of units: 40
## Chromosome 2, # of units: 3
## # of units in total: 43

burden <- seqAssocGLMM_Burden(gdsfile, glmm, units, verbose=FALSE)
head(burden)

##   chr start end maxMAF numvar macmin macmed macmax summac weight        beta
## 1   1     0 499   0.01     14      7     15     19    210  (1,1) -0.03321554
## 2   1     0 499   0.01     14      7     15     19    210 (1,25) -0.03445993
## 3   1     0 499   0.01      2    NaN    NaN    NaN    NaN Cauchy         NaN
## 4   1   250 749   0.01     13      7     14     19    184  (1,1) -0.15157347
## 5   1   250 749   0.01     13      7     14     19    184 (1,25) -0.15100487
## 6   1   250 749   0.01      2    NaN    NaN    NaN    NaN Cauchy         NaN
##           SE      pval method    p.norm converged
## 1 0.09362730 0.7227668 Normal 0.7227668      TRUE
## 2 0.09350925 0.7124863 Normal 0.7124863      TRUE
## 3        NaN 0.7176942   <NA>       NaN      TRUE
## 4 0.10139717 0.1349538 Normal 0.1349538      TRUE
## 5 0.10112838 0.1353849 Normal 0.1353849      TRUE
## 6        NaN 0.1351690   <NA>       NaN      TRUE

acatv <- seqAssocGLMM_ACAT_V(gdsfile, glmm, units, verbose=FALSE)
head(acatv)

##   chr start end maxMAF numvar macmin macmed macmax summac weight n_single
## 1   1     0 499   0.01     14      7     15     19    210  (1,1)       13
## 2   1     0 499   0.01     14      7     15     19    210 (1,25)       13
## 3   1     0 499   0.01      2    NaN    NaN    NaN    NaN Cauchy       NA
## 4   1   250 749   0.01     13      7     14     19    184  (1,1)       10
## 5   1   250 749   0.01     13      7     14     19    184 (1,25)       10
## 6   1   250 749   0.01      2    NaN    NaN    NaN    NaN Cauchy       NA
##   n_collapse      pval
## 1          1 0.1015476
## 2          1 0.1114369
## 3         NA 0.1062713
## 4          3 0.7332731
## 5          3 0.7336823
## 6         NA 0.7334778

SKAT and ACAT-O are not available for survival outcomes: the GATE method defines the score/saddlepoint test for single-variant and burden-style statistics only, and the SKAT variance-component test has no Poisson/Cox analogue. ACAT-O combines burden + ACAT-V + SKAT and is therefore unavailable as well. Both functions stop with an informative error instead of returning a questionable p-value.

seqClose(gdsfile)

Session Information

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] SAIGEgds_2.13.3  Rcpp_1.1.2       SeqArray_1.53.3  gdsfmt_1.49.8   
## [5] Matrix_1.7-6     BiocStyle_2.41.0
## 
## loaded via a namespace (and not attached):
##  [1] jsonlite_2.0.0       compiler_4.6.1       BiocManager_1.30.27 
##  [4] crayon_1.5.3         GenomicRanges_1.65.4 Biostrings_2.81.9   
##  [7] parallel_4.6.1       survey_4.5           jquerylib_0.1.4     
## [10] splines_4.6.1        IRanges_2.47.5       Seqinfo_1.3.2       
## [13] yaml_2.3.12          fastmap_1.2.0        lattice_0.23-1      
## [16] R6_2.6.1             XVector_0.53.0       generics_0.1.4      
## [19] knitr_1.52           BiocGenerics_0.59.12 maketools_1.3.2     
## [22] DBI_1.3.0            bslib_0.12.0         rlang_1.3.0         
## [25] cachem_1.1.0         xfun_0.61            sass_0.4.10         
## [28] sys_3.4.3            RcppParallel_6.2.1   cli_3.6.6           
## [31] SKAT_2.2.5           digest_0.6.39        grid_4.6.1          
## [34] SPAtest_3.1.2        CompQuadForm_1.4.4   lifecycle_1.0.5     
## [37] S4Vectors_0.51.10    RSpectra_0.16-2      evaluate_1.0.5      
## [40] mitools_2.7          buildtools_1.0.0     survival_3.8-12     
## [43] stats4_4.6.1         rmarkdown_2.32       tools_4.6.1         
## [46] htmltools_0.5.9

References

Zheng X, Davis J.Wade. SAIGEgds – an efficient statistical tool for large-scale PheWAS with mixed models. Bioinformatics (2021). Mar;37(5):728-730. DOI: 10.1093/bioinformatics/btaa731.
Dey R, Zhou W, Kiiskinen T, Havulinna A, Elliott A, Karjalainen J, Kurki M, Qin A, FinnGen, Lee S, Palotie A, Neale B, Daly M, Lin X. Efficient and accurate frailty model approach for genome-wide survival association analysis in large-scale biobanks. Nat Commun (2022). 13:5437. DOI: 10.1038/s41467-022-32885-x (the GATE method)
Breslow NE. Covariance analysis of censored survival data. Biometrics (1974). 30(1):89-99.

Time-to-event (survival) analysis in SAIGEgds

Introduction

Input data requirements

Example

Fitting the null model

Sanity check against `coxph()`

Single-variant association tests

Kaplan-Meier curve for the top variant

Set-based (aggregate) tests

Session Information

References

See also

Time-to-event (survival) analysis in SAIGEgds

Introduction

Input data requirements

Example

Fitting the null model

Sanity check against coxph()

Single-variant association tests

Kaplan-Meier curve for the top variant

Set-based (aggregate) tests

Session Information

References

See also

Sanity check against `coxph()`