Learn R Programming

Multivariate Event Times (mets)

Implementation of various statistical models for multivariate event history data doi:10.1007/s10985-013-9244-x. Including multivariate cumulative incidence models doi:10.1002/sim.6016, and bivariate random effects probit models (Liability models) doi:10.1016/j.csda.2015.01.014. Modern methods for survival analysis, including regression modelling (Cox, Fine-Gray, Ghosh-Lin, Binomial regression) with fast computation of influence functions. Restricted mean survival time regression and years lost for competing risks. Average treatment effects and G-computation. All functions can be used with clusters and will work for large data.

Installation

install.packages("mets")

The development version may be installed directly from github (requires Rtools on windows and development tools (+Xcode) for Mac OS X):

remotes::install_github("kkholst/mets", dependencies="Suggests")

or to get development version

remotes::install_github("kkholst/mets",ref="develop")

Citation

To cite the mets package please use one of the following references

Thomas H. Scheike and Klaus K. Holst (2022). A Practical Guide to Family Studies with Lifetime Data. Annual Review of Statistics and Its Application 9, pp. 47-69. doi: http://dx.doi.org/10.1146/annurev-statistics-040120-024253

Thomas H. Scheike and Klaus K. Holst and Jacob B. Hjelmborg (2013). Estimating heritability for cause specific mortality based on twin studies. Lifetime Data Analysis. http://dx.doi.org/10.1007/s10985-013-9244-x

Klaus K. Holst and Thomas H. Scheike Jacob B. Hjelmborg (2015). The Liability Threshold Model for Censored Twin Data. Computational Statistics and Data Analysis. http://dx.doi.org/10.1016/j.csda.2015.01.014

BibTeX:

@Article{,
  title = {A Practical Guide to Family Studies with Lifetime Data},
  author = {Thomas H. Scheike and Klaus K. Holst},
  year = {2014},
  volume = {9},
  pages = {47-69},
  journal = {Annual Review of Statistics and Its Application},
  doi = {10.1146/annurev-statistics-040120-024253},
}

@Article{,
  title={Estimating heritability for cause specific mortality based on twin studies},
  author={Scheike, Thomas H. and Holst, Klaus K. and Hjelmborg, Jacob B.},
  year={2013},
  issn={1380-7870},
  journal={Lifetime Data Analysis},
  doi={10.1007/s10985-013-9244-x},
  url={http://dx.doi.org/10.1007/s10985-013-9244-x},
  publisher={Springer US},
  keywords={Cause specific hazards; Competing risks; Delayed entry;
        Left truncation; Heritability; Survival analysis},
  pages={1-24},
  language={English}
}

@Article{,
  title={The Liability Threshold Model for Censored Twin Data},
  author={Holst, Klaus K. and Scheike, Thomas H. and Hjelmborg, Jacob B.},
  year={2015},
  doi={10.1016/j.csda.2015.01.014},
  url={http://dx.doi.org/10.1016/j.csda.2015.01.014},
  journal={Computational Statistics and Data Analysis}
}

Examples: Twins Polygenic modelling

First considering standard twin modelling (ACE, AE, ADE, and more models)

# simulated data with pairs of observations in twins on long #data format
set.seed(1)
d <- twinsim(1000, b1=c(1,-1), b2=c(), acde=c(1,1,0,1))
# Polygenic model with Additive genetic effects, and shared and invidual environmental effects (ACE)
ace <- twinlm(y ~ 1, data=d, DZ="DZ", zyg="zyg", id="id")
ace
#>        Estimate Std. Error Z value  Pr(>|z|)
#> y     -0.019439   0.041817 -0.4649     0.642
#> sd(A)  0.902004   0.203739  4.4273 9.544e-06
#> sd(C)  1.137025   0.132852  8.5586 < 2.2e-16
#> sd(E)  1.728992   0.037408 46.2194 < 2.2e-16
#> 
#> MZ-pairs DZ-pairs 
#>     1000     1000 
#> 
#> Variance decomposition:
#>   Estimate 2.5%    97.5%  
#> A 0.15966  0.01867 0.30065
#> C 0.25370  0.13920 0.36820
#> E 0.58664  0.53677 0.63650
#> 
#> 
#>                          Estimate 2.5%    97.5%  
#> Broad-sense heritability 0.15966  0.01867 0.30065
#> 
#>                        Estimate 2.5%    97.5%  
#> Correlation within MZ: 0.41336  0.36229 0.46196
#> Correlation within DZ: 0.33353  0.27933 0.38561
#> 
#> 'log Lik.' -8779.953 (df=4)
#> AIC: 17567.91 
#> BIC: 17590.31
# An AE-model could be fitted as
ae <- twinlm(y ~ 1, data=d, DZ="DZ", zyg="zyg", id="id", type="ae")
# AIC
AIC(ae)-AIC(ace)
#> [1] 15.20656
# To adjust for the covariates we simply alter the formula statement
ace2 <- twinlm(y ~ x1+x2, data=d, DZ="DZ", zyg="zyg", id="id", type="ace")
 ## Summary/GOF
summary(ace2)
#>        Estimate Std. Error  Z value Pr(>|z|)
#> y     -0.026049   0.034844  -0.7476   0.4547
#> sd(A)  1.066060   0.072890  14.6256   <2e-16
#> sd(C)  0.980740   0.073569  13.3309   <2e-16
#> sd(E)  0.979980   0.021887  44.7736   <2e-16
#> y~x1   1.006963   0.021900  45.9807   <2e-16
#> y~x2  -0.993802   0.021962 -45.2512   <2e-16
#> 
#> MZ-pairs DZ-pairs 
#>     1000     1000 
#> 
#> Variance decomposition:
#>   Estimate 2.5%    97.5%  
#> A 0.37156  0.27300 0.47012
#> C 0.31446  0.22643 0.40250
#> E 0.31398  0.28381 0.34414
#> 
#> 
#>                          Estimate 2.5%    97.5%  
#> Broad-sense heritability 0.37156  0.27300 0.47012
#> 
#>                        Estimate 2.5%    97.5%  
#> Correlation within MZ: 0.68602  0.65467 0.71502
#> Correlation within DZ: 0.50024  0.45538 0.54257
#> 
#> 'log Lik.' -7449.697 (df=6)
#> AIC: 14911.39 
#> BIC: 14945

Examples: Twins Polygenic modelling time-to-events Data

In the context of time-to-events data we consider the “Liabilty Threshold model” with IPCW adjustment for censoring.

First we fit the bivariate probit model (same marginals in MZ and DZ twins but different correlation parameter). Here we evaluate the risk of getting cancer before the last double cancer event (95 years)

data(prt)
prt0 <-  force.same.cens(prt, cause="status", cens.code=0, time="time", id="id")
prt0$country <- relevel(prt0$country, ref="Sweden")
prt_wide <- fast.reshape(prt0, id="id", num="num", varying=c("time","status","cancer"))
prt_time <- subset(prt_wide,  cancer1 & cancer2, select=c(time1, time2, zyg))
tau <- 95
tt <- seq(70, tau, length.out=5) ## Time points to evaluate model in

b0 <- bptwin.time(cancer ~ 1, data=prt0, id="id", zyg="zyg", DZ="DZ", type="cor",
              cens.formula=Surv(time,status==0)~zyg, breaks=tau)
summary(b0)
#> 
#>                Estimate   Std.Err        Z   p-value    
#> (Intercept)   -1.348188  0.026276 -51.3086 < 2.2e-16 ***
#> atanh(rho) MZ  0.735992  0.087838   8.3789 < 2.2e-16 ***
#> atanh(rho) DZ  0.353023  0.068234   5.1737 2.295e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>  Total MZ/DZ Complete pairs MZ/DZ
#>  1994/3618   997/1809            
#> 
#>                            Estimate 2.5%    97.5%  
#> Tetrachoric correlation MZ 0.62672  0.51081 0.72024
#> Tetrachoric correlation DZ 0.33905  0.21584 0.45164
#> 
#> MZ:
#>                      Estimate 2.5%    97.5%  
#> Concordance          0.03504  0.02779 0.04409
#> Casewise Concordance 0.39458  0.31876 0.47584
#> Marginal             0.08880  0.08086 0.09743
#> Rel.Recur.Risk       4.44351  3.50521 5.38182
#> log(OR)              2.34131  1.87105 2.81157
#> DZ:
#>                      Estimate 2.5%    97.5%  
#> Concordance          0.01952  0.01449 0.02625
#> Casewise Concordance 0.21983  0.16667 0.28415
#> Marginal             0.08880  0.08086 0.09743
#> Rel.Recur.Risk       2.47556  1.81096 3.14016
#> log(OR)              1.23088  0.81020 1.65156
#> 
#>                          Estimate 2.5%    97.5%  
#> Broad-sense heritability 0.57533  0.25790 0.89276
#> 
#> 
#> Event of interest before time 95

Liability threshold model with ACE random effects structure

b1 <- bptwin.time(cancer ~ 1, data=prt0, id="id", zyg="zyg", DZ="DZ", type="ace",
           cens.formula=Surv(time,status==0)~zyg, breaks=tau)
summary(b1)
#> 
#>             Estimate  Std.Err        Z p-value    
#> (Intercept) -2.20664  0.16463 -13.4033  <2e-16 ***
#> log(var(A))  0.43260  0.39149   1.1050  0.2691    
#> log(var(C)) -1.98289  2.52342  -0.7858  0.4320    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>  Total MZ/DZ Complete pairs MZ/DZ
#>  1994/3618   997/1809            
#> 
#>                    Estimate 2.5%     97.5%   
#> A                   0.57533  0.25790  0.89276
#> C                   0.05139 -0.20836  0.31114
#> E                   0.37328  0.26874  0.47782
#> MZ Tetrachoric Cor  0.62672  0.51081  0.72024
#> DZ Tetrachoric Cor  0.33905  0.21584  0.45164
#> 
#> MZ:
#>                      Estimate 2.5%    97.5%  
#> Concordance          0.03504  0.02779 0.04409
#> Casewise Concordance 0.39458  0.31876 0.47584
#> Marginal             0.08880  0.08086 0.09743
#> Rel.Recur.Risk       4.44351  3.50520 5.38182
#> log(OR)              2.34131  1.87104 2.81157
#> DZ:
#>                      Estimate 2.5%    97.5%  
#> Concordance          0.01952  0.01449 0.02625
#> Casewise Concordance 0.21983  0.16667 0.28415
#> Marginal             0.08880  0.08086 0.09743
#> Rel.Recur.Risk       2.47556  1.81095 3.14017
#> log(OR)              1.23088  0.81020 1.65156
#> 
#>                          Estimate 2.5%    97.5%  
#> Broad-sense heritability 0.57533  0.25790 0.89276
#> 
#> 
#> Event of interest before time 95

Examples: Twins Concordance for time-to-events Data


data(prt) ## Prostate data example (sim)

## Bivariate competing risk, concordance estimates
p33 <- bicomprisk(Event(time,status)~strata(zyg)+id(id),
                  data=prt, cause=c(2,2), return.data=1, prodlim=TRUE)
#> Strata 'DZ'
#> Strata 'MZ'

p33dz <- p33$model$"DZ"$comp.risk
p33mz <- p33$model$"MZ"$comp.risk

## Probability weights based on Aalen's additive model (same censoring within pair)
prtw <- ipw(Surv(time,status==0)~country+zyg, data=prt,
            obs.only=TRUE, same.cens=TRUE, 
            cluster="id", weight.name="w")

## Marginal model (wrongly ignoring censorings)
bpmz <- biprobit(cancer~1 + cluster(id), 
                 data=subset(prt,zyg=="MZ"), eqmarg=TRUE)

## Extended liability model
bpmzIPW <- biprobit(cancer~1 + cluster(id),
                    data=subset(prtw,zyg=="MZ"),
                    weights="w")
smz <- summary(bpmzIPW)

## Concordance
plot(p33mz,ylim=c(0,0.1),axes=FALSE, automar=FALSE,atrisk=FALSE,background=TRUE,background.fg="white")
axis(2); axis(1)

abline(h=smz$prob["Concordance",],lwd=c(2,1,1),col="darkblue")
## Wrong estimates:
abline(h=summary(bpmz)$prob["Concordance",],lwd=c(2,1,1),col="lightgray",lty=2)

Examples: Cox model, RMST

We can fit the Cox model and compute many useful summaries, such as restricted mean survival and stanardized treatment effects (G-estimation). First estimating the standardized survival

 data(bmt)
 bmt$time <- bmt$time+runif(408)*0.001
 bmt$event <- (bmt$cause!=0)*1
 dfactor(bmt) <- tcell.f~tcell

 ss <- phreg(Surv(time,event)~tcell.f+platelet+age,bmt) 
 summary(survivalG(ss,bmt,50))
#> G-estimator :
#>       Estimate Std.Err   2.5%  97.5%    P-value
#> risk0   0.6539 0.02708 0.6008 0.7070 9.119e-129
#> risk1   0.5641 0.05973 0.4470 0.6811  3.600e-21
#> 
#> Average Treatment effect: difference (G-estimator) :
#>     Estimate Std.Err    2.5%   97.5% P-value
#> ps0 -0.08982 0.06293 -0.2132 0.03352  0.1535
#> 
#> Average Treatment effect: ratio (G-estimator) :
#> log-ratio: 
#>       Estimate  Std.Err       2.5%      97.5%   P-value
#> ps0 -0.1477619 0.109562 -0.3624994 0.06697567 0.1774462
#> ratio: 
#>  Estimate      2.5%     97.5% 
#> 0.8626365 0.6959347 1.0692695 
#> 
#> Average Treatment effect:  survival-difference (G-estimator) :
#>       Estimate    Std.Err        2.5%     97.5%   P-value
#> ps0 0.08981829 0.06292811 -0.03351854 0.2131551 0.1534889
#> 
#> Average Treatment effect: 1-G (survival)-ratio (G-estimator) :
#> log-ratio: 
#>     Estimate   Std.Err        2.5%     97.5%   P-value
#> ps0 0.230711 0.1504459 -0.06415759 0.5255796 0.1251491
#> ratio: 
#>  Estimate      2.5%     97.5% 
#> 1.2594952 0.9378572 1.6914390

 sst <- survivalGtime(ss,bmt,n=50)
 plot(sst,type=c("survival","risk","survival.ratio")[1])

Based on the phreg we can also compute the restricted mean survival time and years lost (via Kaplan-Meier estimates). The function does it for all times at once and can be plotted as restricted mean survival or years lost at the different time horizons

 out1 <- phreg(Surv(time,cause!=0)~strata(tcell,platelet),data=bmt)
 
 rm1 <- resmean_phreg(out1, times=c(50))
 summary(rm1)
#>                     strata times    rmean se.rmean    lower    upper years.lost
#> tcell=0, platelet=0      0    50 20.48245 1.411055 17.89542 23.44348   29.51755
#> tcell=0, platelet=1      1    50 28.33071 2.196175 24.33733 32.97934   21.66929
#> tcell=1, platelet=0      2    50 22.74596 4.053717 16.04005 32.25544   27.25404
#> tcell=1, platelet=1      3    50 26.11565 4.230688 19.01112 35.87517   23.88435
 par(mfrow=c(1, 2))
 plot(rm1,se=1)
 plot(rm1,years.lost=TRUE,se=1)

For competing risks the years lost can be decomposed into different causes and is based on the integrated Aalen-Johansen estimators for the different strata

 ## years.lost decomposed into causes
 drm1 <- cif_yearslost(Event(time,cause)~strata(tcell,platelet),data=bmt,times=50)
 par(mfrow=c(1,2)); plot(drm1,cause=1,se=1); title(main="Cause 1"); plot(drm1,cause=2,se=1); title(main="Cause 2")
 summary(drm1)
#> $estimate
#> $estimate$intF_1
#>                     strata times   intF_1 se.intF_1 lower_intF_1 upper_intF_1
#> tcell=0, platelet=0      0    50 21.36784  1.476647    18.661106     24.46717
#> tcell=0, platelet=1      1    50 12.97924  2.047516     9.527297     17.68191
#> tcell=1, platelet=0      2    50 12.64543  4.089981     6.708487     23.83649
#> tcell=1, platelet=1      3    50 11.80934  3.673701     6.418453     21.72807
#> 
#> $estimate$intF_2
#>                     strata times    intF_2 se.intF_2 lower_intF_2 upper_intF_2
#> tcell=0, platelet=0      0    50  8.149711  1.094520     6.263606     10.60376
#> tcell=0, platelet=1      1    50  8.690047  1.712441     5.905902     12.78669
#> tcell=1, platelet=0      2    50 14.608610  3.730259     8.856404     24.09685
#> tcell=1, platelet=1      3    50 12.075008  3.890207     6.421784     22.70487
#> 
#> 
#> $total.years.lost
#> [1] 29.51755 21.66929 27.25404 23.88435

Computations are again done for all time horizons at once as illustrated in the plot.

Examples: Cox model IPTW

We can fit the Cox model with inverse probabilty of treatment weights based on logistic regression. The treatment weights can be time-dependent and then mutiplicative weights are applied (see details and vignette).

data(bmt)
bmt$time <- bmt$time+runif(408)*0.001
bmt$id <- seq_len(nrow(bmt))
bmt$event <- (bmt$cause!=0)*1
dfactor(bmt) <- tcell.f~tcell

fit <- phreg_IPTW(Surv(time,event)~tcell.f+cluster(id),data=bmt,treat.model=tcell.f~platelet+age) 
summary(fit)
#> 
#>    n events
#>  408    248
#> 
#>  408 clusters
#> coefficients:
#>           Estimate      S.E.   dU^-1/2 P-value
#> tcell.f1 -0.108497  0.199556  0.089653  0.5867
#> 
#> exp(coefficients):
#>          Estimate    2.5%  97.5%
#> tcell.f1  0.89718 0.60676 1.3266
head(IC(fit))
#>    tcell.f1
#> 1 -1.639241
#> 2 -1.669074
#> 3 -1.749761
#> 4 -1.745988
#> 5 -1.625416
#> 6 -1.793372

Examples: Competing risks regression, Binomial Regression

We can fit the logistic regression model at a specific time-point with IPCW adjustment

data(bmt); bmt$time <- bmt$time+runif(408)*0.001
# logistic regresion with IPCW binomial regression 
out <- binreg(Event(time,cause)~tcell+platelet,bmt,time=50)
summary(out)
#>    n events
#>  408    160
#> 
#>  408 clusters
#> coeffients:
#>              Estimate   Std.Err      2.5%     97.5% P-value
#> (Intercept) -0.180371  0.126757 -0.428811  0.068068  0.1547
#> tcell       -0.418682  0.345438 -1.095729  0.258364  0.2255
#> platelet    -0.436959  0.240977 -0.909266  0.035349  0.0698
#> 
#> exp(coeffients):
#>             Estimate    2.5%  97.5%
#> (Intercept)  0.83496 0.65128 1.0704
#> tcell        0.65791 0.33430 1.2948
#> platelet     0.64600 0.40282 1.0360
head(IC(out))
#>        [,1]     [,2]    [,3]
#> 1 -2.834084 1.633524 2.52025
#> 2 -2.834084 1.633524 2.52025
#> 3 -2.834084 1.633524 2.52025
#> 4 -2.834084 1.633524 2.52025
#> 5 -2.834084 1.633524 2.52025
#> 6 -2.834084 1.633524 2.52025
 predict(out,data.frame(tcell=c(0,1),platelet=c(1,1)),se=TRUE)
#>        pred         se     lower     upper
#> 1 0.3503890 0.04848653 0.2553554 0.4454226
#> 2 0.2619201 0.06969710 0.1253138 0.3985265

Examples: Brier-score for Binomial Regression

 bmt$id <- seq_len(nrow(bmt))

 ## --- competing risks CIF/Survival (default) ---
 fit <- brier_binreg(
   Event(time, cause) ~ tcell + platelet + age + cluster(id),
   data = bmt, time = 50,
   rhs = list(small = ~age, full = ~tcell + platelet + age)
 )
 summary(fit)                        ## log scale
#> 
#> Cross-validated Brier score summary  [log scale]
#> Call: brier_binreg(formula = Event(time, cause) ~ tcell + platelet + 
#>     age + cluster(id), data = bmt, time = 50, rhs = list(small = ~age, 
#>     full = ~tcell + platelet + age))
#> Evaluation time(s): 50   Folds: 5   CI level: 95%
#> Estimator: binregStrata  |  IID: block-diagonal delta-method correction
#> 
#> Models:
#>   null (rhs0)        : ~+1
#>   small              : age
#>   full               : tcell + platelet + age
#> 
#> CV Brier scores by model:
#> 
#>   time = 50
#>          model   Brier     SE   lower   upper
#>  null.bbrierCV -1.4169 0.0195 -1.4551 -1.3787
#>          small -1.4477 0.0274 -1.5013 -1.3941
#>           full -1.4508 0.0308 -1.5112 -1.3905
#> -------------------------------------------------------- 
#> 
#> Model comparisons (delta CV Brier):
#> 
#>   time = 50
#>  time model     reference delta.Brier     SE   lower  upper      P
#>    50 small null.bbrierCV     -0.0308 0.0207 -0.0714 0.0097 0.1363
#>    50  full null.bbrierCV     -0.0339 0.0253 -0.0836 0.0157 0.1806
#>    50  full         small     -0.0031 0.0149 -0.0322 0.0260 0.8345
#> -------------------------------------------------------- 
#> 
#> Note: use summary(fit, transform = exp) for Brier scores on the original scale.
 summary(fit, transform = exp)       ## Brier scale
#> 
#> Cross-validated Brier score summary  [transformed scale]
#> Call: brier_binreg(formula = Event(time, cause) ~ tcell + platelet + 
#>     age + cluster(id), data = bmt, time = 50, rhs = list(small = ~age, 
#>     full = ~tcell + platelet + age))
#> Evaluation time(s): 50   Folds: 5   CI level: 95%
#> Estimator: binregStrata  |  IID: block-diagonal delta-method correction
#> 
#> Models:
#>   null (rhs0)        : ~+1
#>   small              : age
#>   full               : tcell + platelet + age
#> 
#> CV Brier scores by model:
#> 
#>   time = 50
#>          model  Brier     SE  lower  upper
#>  null.bbrierCV 0.2425 0.0047 0.2332 0.2517
#>          small 0.2351 0.0064 0.2225 0.2477
#>           full 0.2344 0.0072 0.2202 0.2485
#> -------------------------------------------------------- 
#> 
#> Model comparisons (delta CV Brier):
#> 
#>   time = 50
#>  time model     reference delta.Brier     SE   lower  upper      P
#>    50 small null.bbrierCV     -0.0074 0.0049 -0.0169 0.0022 0.1313
#>    50  full null.bbrierCV     -0.0081 0.0060 -0.0198 0.0036 0.1742
#>    50  full         small     -0.0007 0.0035 -0.0076 0.0061 0.8343
#> --------------------------------------------------------

Consequently, we can also compute Brier-score for RMST and RMTL regression for survival and competing risks data using the outcome model specification for the binreg function (see examples).

Examples: Competing risks regression, Fine-Gray/Logistic link

We can fit the Fine-Gray model and the logit-link competing risks model (using IPCW adjustment). Starting with the logit-link model

data(bmt)
bmt$time <- bmt$time+runif(nrow(bmt))*0.01
bmt$id <- 1:nrow(bmt)
## logistic link  OR interpretation
 or=cifreg(Event(time,cause)~strata(tcell)+platelet+age,data=bmt,cause=1)
summary(or)
#> 
#>    n events
#>  408    161
#> 
#>  408 clusters
#> coefficients:
#>           Estimate      S.E.   dU^-1/2 P-value
#> platelet -0.454687  0.235394  0.187994  0.0534
#> age       0.390212  0.097676  0.083637  0.0001
#> 
#> exp(coefficients):
#>          Estimate    2.5%  97.5%
#> platelet  0.63465 0.40010 1.0067
#> age       1.47729 1.21990 1.7890
par(mfrow=c(1,2))
 ## to see baseline 
plot(or)

 # predictions 
nd <- data.frame(tcell=c(1,0),platelet=0,age=0)
pll <- predict(or,nd)
plot(pll)

Similarly, the Fine-Gray model can be estimated using IPCW adjustment

 ## Fine-Gray model
 fg=cifreg(Event(time,cause)~strata(tcell)+platelet+age,data=bmt,cause=1,propodds=NULL)
 summary(fg)
#> 
#>    n events
#>  408    161
#> 
#>  408 clusters
#> coefficients:
#>           Estimate      S.E.   dU^-1/2 P-value
#> platelet -0.424370  0.180815  0.187822  0.0189
#> age       0.342156  0.079879  0.086293  0.0000
#> 
#> exp(coefficients):
#>          Estimate    2.5%  97.5%
#> platelet  0.65418 0.45897 0.9324
#> age       1.40798 1.20394 1.6466
## baselines 
plot(fg)
nd <- data.frame(tcell=c(1,0),platelet=0,age=0)
pfg <- predict(fg,nd,se=1)
plot(pfg,se=1)

## influence functions of regression coefficients
head(iid(fg))
#>      platelet          age
#> 1 0.004952983  0.000124585
#> 2 0.005348041 -0.002234595
#> 3 0.006068740 -0.008724149
#> 4 0.006042657 -0.008421442
#> 5 0.004731565  0.001184006
#> 6 0.006330810 -0.012173352

and we can get standard errors for predictions based on the influence functions of the baseline and the regression coefiicients (these are used in the predict function)

baseid <- iidBaseline(fg,time=40)
FGprediid(baseid,nd)
#>           pred     se-log     lower     upper
#> [1,] 0.2786898 0.23979272 0.1741843 0.4458954
#> [2,] 0.4505650 0.07265567 0.3907624 0.5195197

further G-estimation can be done

 dfactor(bmt) <- tcell.f~tcell
 fg1 <- cifreg(Event(time,cause)~tcell.f+platelet+age,bmt,cause=1,propodds=NULL)
 summary(survivalG(fg1,bmt,50))
#> G-estimator :
#>       Estimate Std.Err   2.5%  97.5%   P-value
#> risk0   0.4332 0.02749 0.3793 0.4870 6.332e-56
#> risk1   0.2727 0.05864 0.1577 0.3876 3.318e-06
#> 
#> Average Treatment effect: difference (G-estimator) :
#>     Estimate Std.Err   2.5%    97.5% P-value
#> ps0  -0.1605 0.06354 -0.285 -0.03595 0.01155
#> 
#> Average Treatment effect: ratio (G-estimator) :
#> log-ratio: 
#>       Estimate   Std.Err       2.5%       97.5%    P-value
#> ps0 -0.4627923 0.2212235 -0.8963823 -0.02920229 0.03644142
#> ratio: 
#>  Estimate      2.5%     97.5% 
#> 0.6295234 0.4080432 0.9712200

Examples: Marginal mean for recurrent events

We can estimate the expected number of events non-parametrically and get standard errors for this estimator

data(hfactioncpx12)
dtable(hfactioncpx12,~status)
#> 
#> status
#>    0    1    2 
#>  617 1391  124

gl1 <- recurrentMarginal(Event(entry,time,status)~strata(treatment)+cluster(id),hfactioncpx12,cause=1,death.code=2)
summary(gl1,times=1:5)
#> [[1]]
#>       new.time      mean         se   CI-2.5% CI-97.5% strata
#> 325          1 0.8737156 0.06783343 0.7503858 1.017315      0
#> 555          2 1.5718563 0.09572955 1.3949953 1.771140      0
#> 682          3 2.1184963 0.11385721 1.9066915 2.353829      0
#> 748          4 2.6815219 0.15451005 2.3951619 3.002118      0
#> 748.1        5 2.6815219 0.15451005 2.3951619 3.002118      0
#> 
#> [[2]]
#>       new.time      mean         se   CI-2.5%  CI-97.5% strata
#> 284          1 0.7815557 0.06908585 0.6572305 0.9293989      1
#> 499          2 1.4534055 0.10315606 1.2646561 1.6703258      1
#> 601          3 1.9240624 0.12165771 1.6998008 2.1779119      1
#> 645          4 2.3134997 0.14963892 2.0380418 2.6261880      1
#> 645.1        5 2.3134997 0.14963892 2.0380418 2.6261880      1
plot(gl1,se=1)

Examples: Ghosh-Lin for recurrent events

We can fit the Ghosh-Lin model for the expected number of events observed before dying (using IPCW adjustment and get predictions)

data(hfactioncpx12)
dtable(hfactioncpx12,~status)
#> 
#> status
#>    0    1    2 
#>  617 1391  124

gl1 <- recreg(Event(entry,time,status)~treatment+cluster(id),hfactioncpx12,cause=1,death.code=2)
summary(gl1)
#> 
#>     n events
#>  2132   1391
#> 
#>  741 clusters
#> coefficients:
#>             Estimate      S.E.   dU^-1/2 P-value
#> treatment1 -0.110404  0.078656  0.053776  0.1604
#> 
#> exp(coefficients):
#>            Estimate    2.5%  97.5%
#> treatment1  0.89547 0.76754 1.0447

## influence functions of regression coefficients
head(iid(gl1))
#>      treatment1
#> 1 -1.266428e-04
#> 2 -6.112340e-04
#> 3  2.885192e-03
#> 4  1.308207e-03
#> 5  5.404664e-05
#> 6  2.229380e-03

and we can get standard errors for predictions based on the influence functions of the baseline and the regression coefiicients

 nd=data.frame(treatment=levels(hfactioncpx12$treatment),id=1)
 pfg <- predict(gl1,nd,se=1)
 summary(pfg,times=1:5)
#> Predictions of type 'cumhaz'
#>   Showing subjects: 1, 2
#>   Showing times:    1, 2, 3, 4, 5
#> 
#> -- Subject 1 --
#>  time cumhaz     se  lower  upper
#>     1 0.8573 0.0572 0.7522 0.9771
#>     2 1.5923 0.0882 1.4285 1.7748
#>     3 2.1212 0.1096 1.9169 2.3473
#>     4 2.6354 0.1430 2.3696 2.9311
#>     5 2.6354 0.1430 2.3696 2.9311
#> 
#> -- Subject 2 --
#>  time cumhaz     se  lower  upper
#>     1 0.7677 0.0576 0.6627 0.8894
#>     2 1.4258 0.0950 1.2513 1.6246
#>     3 1.8995 0.1185 1.6809 2.1464
#>     4 2.3600 0.1484 2.0863 2.6695
#>     5 2.3600 0.1484 2.0863 2.6695
 plot(pfg,se=1)

The influence functions of the baseline and regression coefficients at a specific time-point can be obtained

baseid <- iidBaseline(gl1,time=2)
dd <- data.frame(treatment=levels(hfactioncpx12$treatment),id=1)
GLprediid(baseid,dd)
#>          pred     se-log    lower    upper
#> [1,] 1.596065 0.05530215 1.432113 1.778786
#> [2,] 1.429231 0.06660096 1.254329 1.628521

and G-computation

 hfactioncpx12$age <- (50+rnorm(741)*4)[hfactioncpx12$id]

 GLout <- recreg(Event(entry,time,status)~treatment+age,data=hfactioncpx12,cause=1,death.code=2)
 summary(GLout)
#> 
#>     n events
#>  2132   1391
#> 
#>  2132 clusters
#> coefficients:
#>              Estimate       S.E.    dU^-1/2 P-value
#> treatment1 -0.1037358  0.0641060  0.0538440  0.1056
#> age         0.0158460  0.0076257  0.0065148  0.0377
#> 
#> exp(coefficients):
#>            Estimate    2.5%  97.5%
#> treatment1  0.90146 0.79503 1.0222
#> age         1.01597 1.00090 1.0313
 summary(survivalG(GLout,hfactioncpx12,time=4))
#> G-estimator :
#>       Estimate Std.Err  2.5% 97.5%    P-value
#> risk0    2.636  0.1202 2.400 2.872 1.652e-106
#> risk1    2.376  0.1177 2.145 2.607  1.284e-90
#> 
#> Average Treatment effect: difference (G-estimator) :
#>    Estimate Std.Err    2.5%   97.5% P-value
#> p2  -0.2597  0.1602 -0.5736 0.05417  0.1049
#> 
#> Average Treatment effect: ratio (G-estimator) :
#> log-ratio: 
#>      Estimate    Std.Err       2.5%      97.5%   P-value
#> p2 -0.1037358 0.06410601 -0.2293812 0.02190971 0.1056215
#> ratio: 
#>  Estimate      2.5%     97.5% 
#> 0.9014635 0.7950254 1.0221515

Examples: Fixed time modelling for recurrent events

We can fit a log-link regression model at 2 years for the expected number of events observed before dying (using IPCW adjustment)

data(hfactioncpx12)

e2 <- recregIPCW(Event(entry,time,status)~treatment+cluster(id),hfactioncpx12,cause=1,death.code=2,time=2)
summary(e2)
#>    n events
#>  741   1052
#> 
#>  741 clusters
#> coeffients:
#>              Estimate   Std.Err      2.5%     97.5% P-value
#> (Intercept)  0.452430  0.060814  0.333236  0.571624  0.0000
#> treatment1  -0.078322  0.093560 -0.261696  0.105052  0.4025
#> 
#> exp(coeffients):
#>             Estimate    2.5%  97.5%
#> (Intercept)  1.57213 1.39548 1.7711
#> treatment1   0.92467 0.76974 1.1108
head(iid(e2))
#>            [,1]          [,2]
#> 1  1.959479e-04 -2.266440e-04
#> 2  2.237613e-03 -2.227140e-03
#> 3 -9.349773e-06  1.293789e-03
#> 4 -9.653029e-04  9.653029e-04
#> 5 -1.203962e-04  6.744236e-05
#> 6 -2.861359e-03  2.871831e-03

Examples: Regression for RMST/Restricted mean survival for survival and competing risks using IPCW

RMST can be computed using the Kaplan-Meier (via phreg) and the for competing risks via the cumulative incidence functions, but we can also get these estimates via IPCW adjustment and then we can do regression

 ### same as Kaplan-Meier for full censoring model 
 bmt$int <- with(bmt,strata(tcell,platelet))
 out <- resmeanIPCW(Event(time,cause!=0)~-1+int,bmt,time=30,
                         cens.model=~strata(platelet,tcell),model="lin")
 estimate(out)
#>                        Estimate Std.Err  2.5% 97.5%   P-value
#> inttcell=0, platelet=0    13.61  0.8314 11.98 15.24 3.450e-60
#> inttcell=0, platelet=1    18.90  1.2694 16.42 21.39 3.757e-50
#> inttcell=1, platelet=0    16.19  2.4057 11.48 20.91 1.678e-11
#> inttcell=1, platelet=1    17.77  2.4533 12.96 22.58 4.402e-13
 head(iid(out))
#>          [,1] [,2] [,3] [,4]
#> 1 -0.05341304    0    0    0
#> 2 -0.05341074    0    0    0
#> 3 -0.05344765    0    0    0
#> 4 -0.05341970    0    0    0
#> 5 -0.05341166    0    0    0
#> 6 -0.05342025    0    0    0
 ## same as 
 out1 <- phreg(Surv(time,cause!=0)~strata(tcell,platelet),data=bmt)
 rm1 <- resmean_phreg(out1,times=30)
 summary(rm1)
#>                     strata times    rmean  se.rmean    lower    upper
#> tcell=0, platelet=0      0    30 13.60589 0.8314022 12.07017 15.33701
#> tcell=0, platelet=1      1    30 18.90341 1.2691157 16.57270 21.56191
#> tcell=1, platelet=0      2    30 16.19423 2.4002608 12.11150 21.65324
#> tcell=1, platelet=1      3    30 17.76800 2.4418267 13.57249 23.26043
#>                     years.lost
#> tcell=0, platelet=0   16.39411
#> tcell=0, platelet=1   11.09659
#> tcell=1, platelet=0   13.80577
#> tcell=1, platelet=1   12.23200
 
 ## competing risks years-lost for cause 1  
 out1 <- resmeanIPCW(Event(time,cause)~-1+int,bmt,time=30,cause=1,
                       cens.model=~strata(platelet,tcell),model="lin")
 estimate(out1)
#>                        Estimate Std.Err   2.5%  97.5%   P-value
#> inttcell=0, platelet=0   12.103  0.8507 10.436 13.770 6.174e-46
#> inttcell=0, platelet=1    6.883  1.1739  4.582  9.184 4.538e-09
#> inttcell=1, platelet=0    7.259  2.3528  2.648 11.871 2.033e-03
#> inttcell=1, platelet=1    5.779  2.0922  1.679  9.880 5.737e-03
 ## same as 
 drm1 <- cif_yearslost(Event(time,cause)~strata(tcell,platelet),data=bmt,times=30)
 summary(drm1)
#> $estimate
#> $estimate$intF_1
#>                     strata times    intF_1 se.intF_1 lower_intF_1 upper_intF_1
#> tcell=0, platelet=0      0    30 12.102992 0.8506682    10.545456     13.89057
#> tcell=0, platelet=1      1    30  6.882771 1.1738748     4.927069      9.61475
#> tcell=1, platelet=0      2    30  7.259348 2.3528345     3.846040     13.70192
#> tcell=1, platelet=1      3    30  5.779466 2.0921601     2.842849     11.74956
#> 
#> $estimate$intF_2
#>                     strata times   intF_2 se.intF_2 lower_intF_2 upper_intF_2
#> tcell=0, platelet=0      0    30 4.291114 0.6160289     3.238711     5.685489
#> tcell=0, platelet=1      1    30 4.213818 0.9055750     2.765332     6.421022
#> tcell=1, platelet=0      2    30 6.546420 1.9700298     3.629525    11.807501
#> tcell=1, platelet=1      3    30 6.452530 2.0812174     3.429107    12.141686
#> 
#> 
#> $total.years.lost
#> [1] 16.39411 11.09659 13.80577 12.23200

Examples: Average treatment effects (ATE) for survival or competing risks

We can compute ATE for survival or competing risks data for the probabilty of dying

 bmt$event <- bmt$cause!=0; dfactor(bmt) <- tcell~tcell
 brs <- binregATE(Event(time,cause)~tcell+platelet+age,bmt,time=50,cause=1,
      treat.model=tcell~platelet+age)
 summary(brs)
#>    n events
#>  408    160
#> 
#>  408 clusters
#> coeffients:
#>              Estimate   Std.Err      2.5%     97.5% P-value
#> (Intercept) -0.198956  0.130987 -0.455687  0.057774  0.1288
#> tcell1      -0.636967  0.356613 -1.335917  0.061983  0.0741
#> platelet    -0.344844  0.246016 -0.827026  0.137338  0.1610
#> age          0.437263  0.107269  0.227021  0.647506  0.0000
#> 
#> exp(coeffients):
#>             Estimate    2.5%  97.5%
#> (Intercept)  0.81959 0.63401 1.0595
#> tcell1       0.52889 0.26292 1.0639
#> platelet     0.70833 0.43735 1.1472
#> age          1.54846 1.25486 1.9108
#> 
#> Average Treatment effects (G-formula) :
#>             Estimate    Std.Err       2.5%      97.5% P-value
#> treat0     0.4287614  0.0275128  0.3748373  0.4826854  0.0000
#> treat1     0.2899952  0.0659099  0.1608141  0.4191763  0.0000
#> treat:1-0 -0.1387662  0.0717785 -0.2794495  0.0019171  0.0532
#> 
#> Average Treatment effects (double robust) :
#>            Estimate   Std.Err      2.5%     97.5% P-value
#> treat0     0.428174  0.027614  0.374050  0.482297  0.0000
#> treat1     0.250455  0.064791  0.123466  0.377444  0.0001
#> treat:1-0 -0.177719  0.070146 -0.315203 -0.040234  0.0113
 head(brs$riskDR.iid)
#>         iidriskDR     iidriskDR
#> [1,] -0.001158955 -3.553336e-05
#> [2,] -0.001201019  7.581253e-05
#> [3,] -0.001326442  3.358536e-04
#> [4,] -0.001320301  3.246473e-04
#> [5,] -0.001140702 -9.122365e-05
#> [6,] -0.001398214  4.593736e-04
 head(brs$riskG.iid)
#>          iidriskG      iidriskG
#> [1,] -0.001190668 -0.0001527295
#> [2,] -0.001242368  0.0001090930
#> [3,] -0.001355202  0.0006919465
#> [4,] -0.001350615  0.0006680255
#> [5,] -0.001164435 -0.0002837919
#> [6,] -0.001404048  0.0009475763

or the the restricted mean survival or years-lost to cause 1

 out <- resmeanATE(Event(time,event)~tcell+platelet,data=bmt,time=40,treat.model=tcell~platelet)
 summary(out)
#>    n events
#>  408    241
#> 
#>  408 clusters
#> coeffients:
#>              Estimate   Std.Err      2.5%     97.5% P-value
#> (Intercept)  2.852868  0.062474  2.730420  2.975315  0.0000
#> tcell1       0.021394  0.122877 -0.219441  0.262228  0.8618
#> platelet     0.303388  0.090736  0.125549  0.481228  0.0008
#> 
#> exp(coeffients):
#>             Estimate     2.5%   97.5%
#> (Intercept) 17.33743 15.33933 19.5958
#> tcell1       1.02162  0.80297  1.2998
#> platelet     1.35444  1.13377  1.6181
#> 
#> Average Treatment effects (G-formula) :
#>           Estimate  Std.Err     2.5%    97.5% P-value
#> treat0    19.26530  0.95915 17.38539 21.14520  0.0000
#> treat1    19.68189  2.22757 15.31593 24.04785  0.0000
#> treat:1-0  0.41660  2.41045 -4.30779  5.14098  0.8628
#> 
#> Average Treatment effects (double robust) :
#>           Estimate  Std.Err     2.5%    97.5% P-value
#> treat0    19.28438  0.95797 17.40679 21.16196  0.0000
#> treat1    20.34765  2.54076 15.36786 25.32744  0.0000
#> treat:1-0  1.06327  2.70951 -4.24728  6.37382  0.6947
 head(out$riskDR.iid)
#>        iidriskDR  iidriskDR
#> [1,] -0.05143273 0.00589505
#> [2,] -0.05143115 0.00589505
#> [3,] -0.05145648 0.00589505
#> [4,] -0.05143730 0.00589505
#> [5,] -0.05143178 0.00589505
#> [6,] -0.05143768 0.00589505
 head(out$riskG.iid)
#>         iidriskG    iidriskG
#> [1,] -0.05185982 -0.01866233
#> [2,] -0.05185822 -0.01866186
#> [3,] -0.05188377 -0.01866936
#> [4,] -0.05186442 -0.01866368
#> [5,] -0.05185886 -0.01866205
#> [6,] -0.05186481 -0.01866379

 out1 <- resmeanATE(Event(time,cause)~tcell+platelet,data=bmt,cause=1,time=40,
                    treat.model=tcell~platelet)
 summary(out1)
#>    n events
#>  408    157
#> 
#>  408 clusters
#> coeffients:
#>              Estimate   Std.Err      2.5%     97.5% P-value
#> (Intercept)  2.806166  0.069619  2.669715  2.942617  0.0000
#> tcell1      -0.374071  0.247673 -0.859500  0.111359  0.1310
#> platelet    -0.491763  0.164949 -0.815058 -0.168469  0.0029
#> 
#> exp(coeffients):
#>             Estimate     2.5%   97.5%
#> (Intercept) 16.54636 14.43586 18.9654
#> tcell1       0.68793  0.42337  1.1178
#> platelet     0.61155  0.44261  0.8450
#> 
#> Average Treatment effects (G-formula) :
#>           Estimate  Std.Err     2.5%    97.5% P-value
#> treat0    14.52990  0.95694 12.65432 16.40547  0.0000
#> treat1     9.99553  2.37784  5.33505 14.65601  0.0000
#> treat:1-0 -4.53437  2.57483 -9.58095  0.51221  0.0782
#> 
#> Average Treatment effects (double robust) :
#>            Estimate   Std.Err      2.5%     97.5% P-value
#> treat0     14.51179   0.95790  12.63433  16.38924  0.0000
#> treat1      9.36367   2.41670   4.62702  14.10031  0.0001
#> treat:1-0  -5.14812   2.59759 -10.23931  -0.05693  0.0475

Here event is 0/1 thus leading to restricted mean and cause taking the values 0,1,2 produces regression for the years lost due to cause 1.

Examples: While Alive estimands for recurrent events

We consider an RCT and aim to describe the treatment effect via while alive estimands

data(hfactioncpx12)

dtable(hfactioncpx12,~status)
#> 
#> status
#>    0    1    2 
#>  617 1391  124
dd <- WA_recurrent(Event(entry,time,status)~treatment+cluster(id),hfactioncpx12,time=2,death.code=2)
summary(dd)
#> While-Alive summaries:  
#> 
#> RMST,  E(min(D,t)) 
#>            Estimate Std.Err  2.5% 97.5% P-value
#> treatment0    1.859 0.02108 1.817 1.900       0
#> treatment1    1.924 0.01502 1.894 1.953       0
#>  
#>                           Estimate Std.Err    2.5%    97.5% P-value
#> [treatment0] - [treat.... -0.06517 0.02588 -0.1159 -0.01444  0.0118
#> mean events, E(N(min(D,t))): 
#>            Estimate Std.Err  2.5% 97.5%   P-value
#> treatment0    1.572 0.09573 1.384 1.759 1.375e-60
#> treatment1    1.453 0.10315 1.251 1.656 4.376e-45
#>  
#>                           Estimate Std.Err    2.5%  97.5% P-value
#> [treatment0] - [treat....   0.1185  0.1407 -0.1574 0.3943     0.4
#> _______________________________________________________ 
#> Ratio of means E(N(min(D,t)))/E(min(D,t)) 
#>    Estimate Std.Err   2.5%  97.5%   P-value
#> p1   0.8457 0.05264 0.7425 0.9488 4.411e-58
#> p2   0.7555 0.05433 0.6490 0.8619 5.963e-44
#>  
#>             Estimate Std.Err     2.5%  97.5% P-value
#> [p1] - [p2]  0.09022 0.07565 -0.05805 0.2385   0.233
#> _______________________________________________________ 
#> Mean of Events per time-unit E(N(min(D,t))/min(D,t)) 
#>        Estimate Std.Err   2.5%  97.5%   P-value
#> treat0   1.0725  0.1222 0.8331 1.3119 1.645e-18
#> treat1   0.7552  0.0643 0.6291 0.8812 7.508e-32
#>  
#>                     Estimate Std.Err    2.5%  97.5% P-value
#> [treat0] - [treat1]   0.3173  0.1381 0.04675 0.5879 0.02153

dd <- WA_recurrent(Event(entry,time,status)~treatment+cluster(id),hfactioncpx12,time=2,
           death.code=2,trans=.333)
summary(dd,type="log")
#> While-Alive summaries, log-scale:  
#> 
#> RMST,  E(min(D,t)) 
#>            Estimate  Std.Err   2.5%  97.5% P-value
#> treatment0   0.6199 0.011340 0.5977 0.6421       0
#> treatment1   0.6543 0.007807 0.6390 0.6696       0
#>  
#>                           Estimate Std.Err     2.5%     97.5% P-value
#> [treatment0] - [treat.... -0.03446 0.01377 -0.06145 -0.007478 0.01231
#> mean events, E(N(min(D,t))): 
#>            Estimate Std.Err   2.5%  97.5%   P-value
#> treatment0   0.4523 0.06090 0.3329 0.5716 1.119e-13
#> treatment1   0.3739 0.07097 0.2348 0.5130 1.376e-07
#>  
#>                           Estimate Std.Err    2.5%  97.5% P-value
#> [treatment0] - [treat....  0.07835 0.09352 -0.1049 0.2616  0.4022
#> _______________________________________________________ 
#> Ratio of means E(N(min(D,t)))/E(min(D,t)) 
#>            Estimate Std.Err    2.5%    97.5%   P-value
#> treatment0  -0.1676 0.06224 -0.2896 -0.04563 7.081e-03
#> treatment1  -0.2804 0.07192 -0.4214 -0.13947 9.651e-05
#>  
#>                           Estimate Std.Err     2.5%  97.5% P-value
#> [treatment0] - [treat....   0.1128 0.09511 -0.07361 0.2992  0.2356
#> _______________________________________________________ 
#> Mean of Events per time-unit E(N(min(D,t))/min(D,t)) 
#>        Estimate Std.Err    2.5%   97.5%   P-value
#> treat0  -0.3833 0.04939 -0.4801 -0.2865 8.487e-15
#> treat1  -0.5380 0.05666 -0.6491 -0.4270 2.191e-21
#>  
#>                     Estimate Std.Err     2.5%  97.5% P-value
#> [treat0] - [treat1]   0.1548 0.07517 0.007459 0.3021 0.03948

Copy Link

Version

Install

install.packages('mets')

Monthly Downloads

17,484

Version

1.3.12

License

Apache License (== 2.0)

Issues

Pull Requests

Stars

Forks

Maintainer

Klaus Holst

Last Published

July 30th, 2026

Functions in mets (1.3.12)

binregCasewise

Estimate Casewise Concordance Using Binomial Regression
binregATE

Average Treatment Effect for Censored Competing Risks Data using Binomial Regression
binregG

G-Estimator for Binomial Regression Model (Standardized Estimates)
WA_reg

While-Alive Regression for Recurrent Events
binregRatio

Percentage of Years Lost Due to a Cause Regression
WA_recurrent

While-Alive Estimands for Recurrent Events
aalenMets

Fast Additive Hazards Model with Robust Standard Errors
binreg

Binomial Regression for Censored Competing Risks Data
bicomprisk

Estimation of Concordance in Bivariate Competing Risks Data
binomial_twostage

Fits Clayton-Oakes or bivariate Plackett (OR) models for binary data using marginals that are on logistic form. If clusters contain more than two times, the algoritm uses a compososite likelihood based on all pairwise bivariate models.
binregTSR

Two-Stage Randomization for Survival or Competing Risks Data
blocksample

Block sampling
binregStrata

Stratified Binomial Regression with cross-validation option
casewise

Estimate Casewise Concordance from prodlim Objects
biprobit

Bivariate Probit model
binreg_IPTW

IPTW logistic regression, Inverse Probabibilty of Treatment Weighted binreg
brier_binreg

Cross-validated Brier score for competing risks, RMST and RMTL regression using stratified leave-fold-out fits (binregStrata)
bptwin

Liability model for twin data
bmt

The Bone Marrow Transplant Data
calgb8923

CALGB 8923, twostage randomization SMART design
casewise_bin

Casewise Concordance from Concordant/Discordant Counts
cifregFG

Fine-Gray Cumulative Incidence Function Regression
coarse_clust

Coarsen Cluster Identifiers
cluster_index

Finds subjects related to same cluster
concordanceCor

Concordance Computes concordance and casewise concordance
rr_cif

Non-parametric Cumulative Incidence Functions
cif

Cumulative Incidence with Robust Standard Errors
cifreg

Cumulative Incidence Function (CIF) Regression
cif_yearslost

Restricted Mean Time Lost for Competing Risks
cor_cif

Cross-odds-ratio, OR or RR risk regression for competing risks
dcor

summary, tables, and correlations for data frames
count_history

Compute cumulative event counts as time-dependent covariates
dermalridgesMZ

Dermal ridges data (monozygotic twins)
divide_conquer

Split a data set and run function
daggregate

aggregating for for data frames
cumoddsreg

Cumulative Odds Regression for Discrete Time Data
dcut

Cutting, sorting, rm (removing), rename for data frames
dby

Calculate summary statistics grouped by
diabetes

The Diabetic Retinopathy Data
dermalridges

Dermal ridges data (families)
drelevel

relev levels for data frames
dtransform

Transform that allows condition
dsort

Sort data frame
dspline

Simple linear spline
drop.specials

Remove Special Terms from a Formula
estimate.binregCV

Joint estimate across CV models for a binregCV object
dlag

Lag operator
dprint

list, head, print, tail
dreg

Regression for data frames with dutility call
dtable

tables for data frames
familyclusterWithProbands_index

Finds all pairs within a cluster (famly) with the proband (case/control)
fast.pattern

Fast pattern
fast.cluster

Fast Cluster Index Conversion
fast.reshape

Fast reshape
event_split

event_split (SurvSplit).
familycluster_index

Finds all pairs within a cluster (family)
eventpois

Extract survival estimates from lifetable analysis
fast.approx

Fast approximation
extendCums

Extend Cumulative Hazard Functions to Common Time Range
event_split2

Event split with two time-scales, time and gaptime
faster.reshape

Fast Reshape from Long to Wide Format
haplo

haplo fun data
gofM_phreg

Goodness-of-Fit for Cox Covariates (Model Matrix)
grouptable

Create Group Contingency Table from Clustered Data
gof.phreg

Goodness-of-Fit for Cox PH Regression (Proportionality)
gofZ_phreg

Goodness-of-Fit for Cox Covariates (Linearity)
haplo_surv_discrete

Discrete Time-to-Event Haplotype Analysis
force.same.cens

Force Same Censoring Within Clusters
folds

Generate Random Fold Indices for Cross-Validation
hfactioncpx12

hfaction, subset of block randomized study HF-ACtion from WA package
lifetable.matrix

Life table
iidBaseline

Influence Functions or IID Decomposition of Baseline
jumptimes

Extract Event (Jump) Times
ilap

Inverse Laplace Transform Helper
ipw2

Inverse Probability of Censoring Weights
lifecourse

Life-course plot
interval_logitsurv_discrete

Discrete Time-to-Event Analysis with Interval Censoring
ipw

Inverse Probability of Censoring Weights
km

Kaplan-Meier with Robust Standard Errors
logitSurv

Proportional Odds Survival Model
migr

Migraine data
mlogit

Multinomial Regression Based on phreg
medweight

Computes mediation weights
mediatorSurv

Mediation analysis in survival context
multcif

Multivariate Cumulative Incidence Function example data set
mets-package

mets: Analysis of Multivariate Event Times
melanoma

The Melanoma Survival Data
mets-simulation

Simulation Helper Functions
mets.options

Set global options for mets
mena

Menarche data set
plot.phreg

Plotting the baselines of stratified Cox
robust.basehaz.phreg

Robust Baseline Hazard Standard Errors
plot_twin

Scatter plot function
phreg_rct

Lu-Tsiatis More Efficient Log-Rank for Randomized Studies with Baseline Covariates
plack_cif

plack Computes concordance for or.cif based model, that is Plackett random effects model
np

np data set
phreg

Fast Cox Proportional Hazards Regression
phreg_IPTW

IPTW Cox Regression (Inverse Probability of Treatment Weighted)
pmvn

Multivariate normal distribution function
phreg_weibull

Weibull-Cox regression
prob_exceed_recurrent

Estimate the probability of exceeding k recurrent events by time t
print.casewise

prints Concordance test
prt

Prostate data set
predict.phreg

Predictions from Proportional Hazards Model
predict.mlogit

Predictions from Multinomial Regression
print.summary_binregCV_multi

Print method for summary_binregCV_multi objects
rchazl

Multiple Cause Piecewise Constant Hazard Simulation
ratioATE

Ratio of Average Treatment Effects
rchaz

Simulation of Piecewise Constant Hazard Model (Cox)
random_cif

Random effects model for competing risks data
recregIPCW

IPCW Estimator for Recurrent Events
recurrent_marginal

Marginal mean estimation for recurrent events with a terminal event
sim_ClaytonOakes

Simulate from the Clayton-Oakes frailty model
resmean_phreg

Restricted Mean for Stratified Kaplan-Meier or Cox Model
resmeanATE

Average Treatment Effect for Restricted Mean Time
rweibullcox

Simulate observations from a Weibull distribution
rcrisk

Simulation of Piecewise constant hazard models with two causes (Cox).
recreg

Recurrent Events Regression with Terminal Event
resmeanIPCW

Restricted IPCW Mean for Censored Survival Data
reexports

Objects exported from other packages
sim_multistateII

Illness-Death Competing Risks with Two Causes of Death
sim_recurrent

Simulate recurrent events with a single event type and a terminal event
sim_GLcox

Simulation of Two-Stage Recurrent Events Data
sim_ClaytonOakesWei

Simulate from the Clayton-Oakes frailty model
sim_recurrentII

Simulate recurrent events with two event types and a terminal event
sim_recurrentTS

Simulate recurrent events from a two-stage model with structured gamma frailties
sim_cif

Simulation of Output from Cumulative Incidence Regression Model
sim_phreg

Simulation of Output from Cox Model
sim_phregs

Simulation of Cause-Specific Cox Models
sim_multistate

Simulation of Illness-Death Model
survival-helpers

Survival Twostage Helpers
summaryGLM

Reporting OR (exp(coef)) from glm with binomial link and glm predictions
strata-numeric

Stratified Cumulative and Summary Operations
summary.cor

Summary for dependence models for competing risks
surv_boxarea

Bivariate Survival Data on Rectangular Regions
sim_recurrent_ts

Simulate recurrent events from a two-stage Cox or Ghosh-Lin model
survivalG

G-Estimator for Cox and Fine-Gray Models
summaryTimeobject

Summarize a Time-Varying Estimate with Confidence Bands
summary.binregCV

Summary method for binregCV objects
summary.binregCV_list

Summary method for binregCV_list objects
p11_binomial_twostage_RV

Concordance Probability from Twostage Model
test_marginalMean

Pepe-Mori Test for Marginal Mean Comparison
test_conc

Compare Two Concordance Estimates
tie_breaker

Break ties in event times for recurrent event data
tetrachoric

Estimate parameters from odds-ratio
test_casewise

Test for Independence Using Casewise Concordance
survival_twostage

Twostage Survival Model for Multivariate Survival Data
test_logrankRecurrent

Logrank-type test for comparing recurrent event marginal means between groups
ttpd

ttpd discrete survival data on interval form
twinbmi

BMI data set
twinstut

Stutter data set
twostageREC

Fitting of Two-Stage Recurrent Events Random Effects Model
twostageMLE

Twostage Survival Model Fitted by Pseudo MLE
twinlm

Classic twin model for quantitative traits
twinsim

Simulate twin data
LinSpline

Simple linear spline
Event

Event history object
IC.binreg

Influence curve components for binomial regression ATE
ClaytonOakes

Clayton-Oakes model with piece-wise constant hazards
IC.phreg

Influence Functions for phreg objects
Dbvn

Derivatives of the bivariate normal cumulative distribution function
CPH_HPN_CRBSI

Rates for HPN program for patients of Copenhagen Cohort
Grandom_cif

Additive Random effects model for competing risks data for polygenetic modelling
TRACE

The TRACE study group of myocardial infarction
ACTG175

ACTG175, block randomized study from speff2trial package