Introduction
xhaz is an R package for fitting excess hazard regression models, with or without the assumption of proportional hazards. In the presence of hierarchical data, inter-cluster heterogeneity can affect excess hazard estimates. To address this issue, the mexhaz R package implements multilevel excess hazard regression models. If the user wants to check the validity of the non-comparability bias, an easy way is to use the mexhazLT function of the xhaz R package. The user can also specify whether the framework is consistent with classical excess hazard modeling, i.e., assuming that the expected mortality of the individuals studied is appropriate, using the pophaz argument is equivalent to classical. The user can also consider a different framework by specifying that pophaz equals rescaled: the expected mortality available in the life table is not accurate, i.e., there is a non-comparability source of bias with respect to the expected mortality of the study population Goungounga et al. (2023).
(Here’s the abstract from Goungounga et al. (2023) paper: In the presence of competing causes of event occurrence (e.g., death), the interest might not only be in the overall survival but also in the so-called net survival, that is, the hypothetical survival that would be observed if the disease under study were the only possible cause of death. Net survival estimation is commonly based on the excess hazard approach in which the hazard rate of individuals is assumed to be the sum of a disease-specific and expected hazard rate, supposed to be correctly approximated by the mortality rates obtained from general population life tables…
A related software package can be found at https://CRAN.R-project.org/package=xhaz.
Installation
The most recent version of xhaz can be installed
directly from the cran repository using
xhaz depends on the stats,
survival, optimParallel,
numDeriv, statmod, gtools and
splines packages which can be installed directly from
CRAN.
It also utilizes survexp.fr, the R package containing
the French life table. For example, to install survexp.fr
follow the instructions available at the RStudio page on R
and survexp.fr.
First, install the R package via github.
Then, when these other packages are installed, please load the xhaz R package.
Fitting an classical excess hazard regression model with a cubic B-splines for the baseline hazard as in mexhaz
For this illustration, we used a simulated dataset from the original
paper by Goungounga et al. (2023). This dataset consists of 4,978
patients with information on their treatment arm and clinical center ID,
as this information may affect the excess mortality rate. The US life
table can be used to estimate the model parameters. mexhazLT() function
will be used with argument (`pophaz = classic).
library(survexp.fr)
#> Warning: le package 'survexp.fr' a été compilé avec la version R 4.5.2
data("breast", package = "xhaz")
head(breast)
#> temps statut age agecr date SEX armt hosp dept
#> 3 5.050550 1 30.95628 -2.213589 1985-12-23 2 0 3 09
#> 8 1.415929 1 35.20680 -1.788538 1985-04-30 2 1 8 81
#> 11 1.904082 1 31.80505 -2.128713 1985-10-15 2 0 11 06
#> 16 15.000000 0 42.52067 -1.057150 1985-09-16 2 1 16 10
#> 17 2.504911 0 32.01450 -2.107767 1985-01-29 2 0 17 94
#> 18 5.319706 0 31.77941 -2.131276 1985-06-04 2 1 18 55
dim(breast)
#> [1] 4978 9
# armt: numbers and percentages
breast$armt2 <- factor(breast$armt, levels = c("0","1"))
table(breast$armt2)
#>
#> 0 1
#> 2471 2507
# The life table to be used is survexp.us. Note that SEX is coded 2 instead of female in survexp.us.
breast$sexe <- "female"
fit.haz <- exphaz(
formula = Surv(temps, statut) ~ 1,
data = breast, ratetable = survexp.us,
only_ehazard = FALSE,
rmap = list(age = 'age', sex = 'sexe',
year = 'date'))
breast$expected <- fit.haz$ehazard
breast$expectedCum <- fit.haz$ehazardInt
qknots <- quantile(breast[breast$statut==1,]$temps, probs=c(1:2/3))
mod.bs <- mexhazLT(
formula = Surv(temps, statut) ~ agecr + armt2,
data = breast, degree = 3,
knots = qknots, expected = "expected",
expectedCum = "expectedCum",
base = "exp.bs", pophaz = "classic")
#> iteration = 0
#> Step:
#> [1] 0 0 0 0 0 0 0 0
#> Parameter:
#> [1] -1 -1 -1 -1 -1 -1 0 0
#> Function Value
#> [1] 11032.51
#> Gradient:
#> [1] -1281.60395 -490.49700 -810.60363 -112.90836 72.86473 103.05848
#> [7] -1994.09940 -415.61972
#>
#> iteration = 39
#> Parameter:
#> [1] -1.5083019 0.8413747 0.5236963 -0.2330524 -0.3193138 -1.1142808 0.3759031
#> [8] -0.3239171
#> Function Value
#> [1] 9784.837
#> Gradient:
#> [1] -0.0003931511 0.0013787940 -0.0010113581 -0.0001018634 -0.0007694325
#> [6] 0.0007721411 -0.0000509317 0.0001600711
#>
#> Gradient relatif proche de zéro.
#> L'itération courante est probablement la solution.
#>
#>
#> Data
#> Name N.Obs.Tot N.Obs N.Events N.Clust
#> data 4978 4978 4363 1
#>
#> Details
#> Iter Eval Base Nb.Leg Nb.Aghq Optim Method Code LogLik Total.Time
#> 39 289 exp.bs 20 10 nlm --- 1 -9784.837 2.03
mod.bs
#> Call:
#> mexhazLT(formula = Surv(temps, statut) ~ agecr + armt2, data = breast,
#> expected = "expected", expectedCum = "expectedCum", pophaz = "classic",
#> base = "exp.bs", degree = 3, knots = qknots)
#>
#> Coefficients:
#> Estimate StdErr z.value p.value
#> Intercept -1.508302 0.088649 -17.0143 < 2.2e-16 ***
#> BS3.1 0.841375 0.131145 6.4156 1.403e-10 ***
#> BS3.2 0.523696 0.090616 5.7793 7.501e-09 ***
#> BS3.3 -0.233052 0.188543 -1.2361 0.2164
#> BS3.4 -0.319314 0.255004 -1.2522 0.2105
#> BS3.5 -1.114281 0.248572 -4.4827 7.369e-06 ***
#> agecr 0.375903 0.013603 27.6334 < 2.2e-16 ***
#> armt21 -0.323917 0.031372 -10.3249 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> log-likelihood: -9784.837 (for 8 degree(s) of freedom)
#>
#> number of observations: 4978, number of events: 4363Fitting a rescaled excess hazard regression model with a cubic B-splines for the baseline hazard as in mexhaz
The new parameter to be added to mexhazLT() function is
pophaz="rescaled": it allows to rescale the life table
available for the estimation of the excess hazard parameters. This model
concerns that proposed by Goungounga et al. (2016).
mod.bs2 <- mexhazLT(
formula = Surv(temps, statut) ~ agecr + armt2,
data = breast, degree = 3,
knots = qknots, expected = "expected",
expectedCum = "expectedCum",
base = "exp.bs", pophaz = "rescaled")
#> iteration = 0
#> Step:
#> [1] 0 0 0 0 0 0 0 0 0
#> Parameter:
#> [1] -1 -1 -1 -1 -1 -1 0 0 0
#> Function Value
#> [1] 11032.51
#> Gradient:
#> [1] -1281.60395 -490.49700 -810.60363 -112.90836 72.86473 103.05848
#> [7] -1994.09940 -415.61972 -140.32355
#>
#> iteration = 51
#> Parameter:
#> [1] -1.5243849 0.8517578 0.5329562 -0.2609964 -0.3386736 -1.2554154 0.3655821
#> [8] -0.3321785 0.5658880
#> Function Value
#> [1] 9784.532
#> Gradient:
#> [1] 1.523795e-03 4.783942e-04 6.621121e-04 3.292371e-04 5.093170e-05
#> [6] 5.071200e-05 -2.364686e-05 8.985808e-04 1.746230e-04
#>
#> Gradient relatif proche de zéro.
#> L'itération courante est probablement la solution.
#>
#>
#> Data
#> Name N.Obs.Tot N.Obs N.Events N.Clust
#> data 4978 4978 4363 1
#>
#> Details
#> Iter Eval Base Nb.Leg Nb.Aghq Optim Method Code LogLik Total.Time
#> 51 392 exp.bs 20 10 nlm --- 1 -9784.532 2.75
mod.bs2
#> Call:
#> mexhazLT(formula = Surv(temps, statut) ~ agecr + armt2, data = breast,
#> expected = "expected", expectedCum = "expectedCum", pophaz = "rescaled",
#> base = "exp.bs", degree = 3, knots = qknots)
#>
#> Coefficients:
#> Estimate StdErr z.value p.value
#> Intercept -1.524385 0.092895 -16.4097 < 2.2e-16 ***
#> BS3.1 0.851758 0.134411 6.3370 2.344e-10 ***
#> BS3.2 0.532956 0.093281 5.7134 1.107e-08 ***
#> BS3.3 -0.260996 0.198195 -1.3169 0.187883
#> BS3.4 -0.338674 0.271588 -1.2470 0.212393
#> BS3.5 -1.255415 0.327769 -3.8302 0.000128 ***
#> agecr 0.365582 0.019092 19.1480 < 2.2e-16 ***
#> armt21 -0.332178 0.033808 -9.8253 < 2.2e-16 ***
#> log(Alpha) 0.565888 0.552889 1.0235 0.306066
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> log-likelihood: -9784.5317 (for 9 degree(s) of freedom)
#>
#> number of observations: 4978, number of events: 4363Fitting a random-effects excess hazard regression model with a cubic B-splines for the baseline hazard as in mexhaz
As in mexhaz, the new parameter to be added to the
mexhazLT() function is random="hosp": it
allows to specify the variable indicating the cluster levels. However,
pophaz argument is set to be equal to
"classic". This excess hazard regression model is the one
proposed by Charvat et al (2016). In the output, loglik
corresponds to the total log-likelihood including the sum of the
cumulative expected hazards, which is often removed in classical excess
hazard regression models because it is considered a nuisance
parameter.
mod.bs3 <- mexhazLT(
formula = Surv(temps, statut) ~ agecr + armt2,
data = breast, degree = 3,
knots = qknots, expected = "expected",
expectedCum = "expectedCum",
base = "exp.bs", pophaz = "classic", random = "hosp")
#> iteration = 0
#> Step:
#> [1] 0 0 0 0 0 0 0 0 0
#> Parameter:
#> [1] -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.0 0.0 0.1
#> Function Value
#> [1] 10129.93
#> Gradient:
#> [1] -39.69633 -148.95592 -395.16138 18.72774 118.58665 119.00834
#> [7] -2896.28313 329.60901 39.47875
#>
#> iteration = 53
#> Parameter:
#> [1] -1.90783127 1.00354112 1.11855695 0.90280950 0.65681302 0.05890742
#> [7] 0.49124655 -0.45036907 -0.20709299
#> Function Value
#> [1] 8999.233
#> Gradient:
#> [1] 0.0026219409 0.0022965273 0.0006163271 0.0006457412 -0.0012478267
#> [6] 0.0002801244 -0.0007330527 -0.0002019078 0.0019390427
#>
#> Gradient relatif proche de zéro.
#> L'itération courante est probablement la solution.
#>
#> Computation of the covariance matrix of the shrinkage estimators
#>
#> Data
#> Name N.Obs.Tot N.Obs N.Events N.Clust
#> data 4978 4978 4363 100
#>
#> Details
#> Iter Eval Base Nb.Leg Nb.Aghq Optim Method Code LogLik Total.Time
#> 53 422 exp.bs 20 10 nlm --- 1 -8999.233 16.95
mod.bs3
#> Call:
#> mexhazLT(formula = Surv(temps, statut) ~ agecr + armt2, data = breast,
#> expected = "expected", expectedCum = "expectedCum", pophaz = "classic",
#> base = "exp.bs", degree = 3, knots = qknots, random = "hosp")
#>
#> Coefficients:
#> Estimate StdErr z.value p.value
#> Intercept -1.907831 0.122002 -15.6376 < 2.2e-16 ***
#> BS3.1 1.003541 0.131641 7.6233 2.465e-14 ***
#> BS3.2 1.118557 0.093232 11.9976 < 2.2e-16 ***
#> BS3.3 0.902810 0.192646 4.6864 2.781e-06 ***
#> BS3.4 0.656813 0.257079 2.5549 0.010622 *
#> BS3.5 0.058907 0.251562 0.2342 0.814855
#> agecr 0.491247 0.014198 34.5990 < 2.2e-16 ***
#> armt21 -0.450369 0.032262 -13.9596 < 2.2e-16 ***
#> hosp [log(sd)] -0.207093 0.074504 -2.7796 0.005442 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> log-likelihood: -8999.2334 (for 9 degree(s) of freedom)
#>
#> number of observations: 4978, number of events: 4363Fitting a rescaled random-effects excess hazard regression model with a cubic B-splines for the baseline hazard as in mexhaz
As in mexhaz, the new parameter to be added to the
mexhazLT() function is random="hosp": it
allows to specify the variable indicating the cluster levels. However,
pophaz is set to be equal to “rescaled”
(`pophaz = "rescaled"). This excess hazard regression model
is the one proposed by Goungounga et al (2023).
mod.bs4 <- mexhazLT(
formula = Surv(temps, statut) ~ agecr + armt2,
data = breast, degree = 3,
knots = qknots, expected = "expected",
expectedCum = "expectedCum",
base = "exp.bs", pophaz = "rescaled", random = "hosp")
#> iteration = 0
#> Step:
#> [1] 0 0 0 0 0 0 0 0 0 0
#> Parameter:
#> [1] -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.0 0.0 0.0 0.1
#> Function Value
#> [1] 10129.93
#> Gradient:
#> [1] -39.69633 -148.95592 -395.16138 18.72774 118.58665 119.00834
#> [7] -2896.28313 329.60901 -82.81752 39.47875
#>
#> iteration = 63
#> Parameter:
#> [1] -1.9623310 1.0285195 1.1490601 0.8954798 0.6356856 -0.1502895
#> [7] 0.4755500 -0.4705399 0.9120266 -0.1643994
#> Function Value
#> [1] 8997.261
#> Gradient:
#> [1] -2.002219e-04 8.489046e-05 2.263724e-04 6.366463e-05 -4.129106e-04
#> [6] 2.237357e-04 -1.246008e-03 -6.893970e-04 -4.947651e-04 1.564331e-04
#>
#> Gradient relatif proche de zéro.
#> L'itération courante est probablement la solution.
#>
#> Computation of the covariance matrix of the shrinkage estimators
#>
#> Data
#> Name N.Obs.Tot N.Obs N.Events N.Clust
#> data 4978 4978 4363 100
#>
#> Details
#> Iter Eval Base Nb.Leg Nb.Aghq Optim Method Code LogLik Total.Time
#> 63 532 exp.bs 20 10 nlm --- 1 -8997.261 21.1
mod.bs4
#> Call:
#> mexhazLT(formula = Surv(temps, statut) ~ agecr + armt2, data = breast,
#> expected = "expected", expectedCum = "expectedCum", pophaz = "rescaled",
#> base = "exp.bs", degree = 3, knots = qknots, random = "hosp")
#>
#> Coefficients:
#> Estimate StdErr z.value p.value
#> Intercept -1.962331 0.130425 -15.0457 < 2.2e-16 ***
#> BS3.1 1.028519 0.137173 7.4980 6.484e-14 ***
#> BS3.2 1.149060 0.098175 11.7042 < 2.2e-16 ***
#> BS3.3 0.895480 0.205051 4.3671 1.259e-05 ***
#> BS3.4 0.635686 0.285647 2.2254 0.026053 *
#> BS3.5 -0.150290 0.310309 -0.4843 0.628157
#> agecr 0.475550 0.016621 28.6121 < 2.2e-16 ***
#> armt21 -0.470540 0.035247 -13.3499 < 2.2e-16 ***
#> log(Alpha) 0.912027 0.304652 2.9937 0.002756 **
#> hosp [log(sd)] -0.164399 0.078232 -2.1014 0.035604 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> log-likelihood: -8997.2614 (for 10 degree(s) of freedom)
#>
#> number of observations: 4978, number of events: 4363We can compare the output of these two models using AIC or BIC criteria.
compared_models <- list(mod.bs,mod.bs2, mod.bs3, mod.bs4)
names(compared_models) <- c("mod.bs","mod.bs2", "mod.bs3", "mod.bs4")
sapply(compared_models, function(i) {
AIC(i)
})
#> mod.bs mod.bs2 mod.bs3 mod.bs4
#> 19585.67 19587.06 18016.47 18014.52A statistical comparison between two nested models can be performed with a likelihood ratio test calculated by function anova method implemented in xhaz.
For example, suppose we want to test whether we can drop the rescaling parameter between the different excess hazard regression models.
As in survival package, we compare the models using anova(), i.e.,
anova(mod.bs,mod.bs2)
#> Assumption: Model 1 nested within Model 2
#>
#> Likelihood ratio test
#> Model 1:
#> Surv(temps, statut) ~ agecr + armt2
#> Model 2:
#> Surv(temps, statut) ~ agecr + armt2
#> Model.df loglik df Chisq Pr(>Chisq)
#> [1,] 8.0 -9784.8 NA NA NA
#> [2,] 9.0 -9784.5 1.0 0.305 0.4346
anova(mod.bs3,mod.bs4)
#> Assumption: Model 1 nested within Model 2
#>
#> Likelihood ratio test
#> Model 1:
#> Surv(temps, statut) ~ agecr + armt2
#> Model 2:
#> Surv(temps, statut) ~ agecr + armt2
#> Model.df loglik df Chisq Pr(>Chisq)
#> [1,] 9.0 -8999.2 NA NA NA
#> [2,] 10.0 -8997.3 1.0 1.972 0.04704 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Note that it is the user’s responsibility to provide properly nested hazard models for the LRT to be valid. The results suggest that by correcting for between-cluster heterogeneity and non-comparability bias with a rescaled random effects excess hazard regression model, the user will be able to provide more accurate estimates of net survival.
Plot of net survival and excess hazard for different models
One could be interested in the prediction of net survival and excess hazard for the individual with the same characteristics as individual 1 in the breast dataset (age 30.95, armt equals 0) as performed in mexhaz.
# hazard and survival in the two arms
predict_mod_amr0 <- predict(mod.bs,
time.pts=seq(0.1,10,by=0.1),
data.val=data.frame(agecr = 0,
armt2 = "0"))
predict_mod_amr1 <- predict(mod.bs,
time.pts=seq(0.1,10,by=0.1),
data.val=data.frame(agecr = 0,
armt2 = "1"))
predict_mod2_arm0 <- predict(mod.bs2,
time.pts=seq(0.1,10,by=0.1),
data.val=data.frame(agecr = 0,
armt2 ="0"))
predict_mod2_arm1 <- predict(mod.bs2,
time.pts=seq(0.1,10,by=0.1),
data.val=data.frame(agecr = 0,
armt2 ="1"))
predict_mod3_arm0 <- predict(mod.bs3,
time.pts=seq(0.1,10,by=0.1),
data.val=data.frame(agecr = 0,
armt2 = "0"),
marginal = TRUE)
predict_mod3_arm1 <- predict(mod.bs3,
time.pts=seq(0.1,10,by=0.1),
data.val=data.frame(agecr = 0,
armt2 = "1"),
marginal = TRUE)
predict_mod4_arm0 <- predict(mod.bs4,
time.pts=seq(0.1,10,by=0.1),
data.val=data.frame(agecr = 0,
armt2 = "0"),
marginal = TRUE)
predict_mod4_arm1 <- predict(mod.bs4,
time.pts=seq(0.1,10,by=0.1),
data.val=data.frame(agecr = 0,
armt2 = "1"),
marginal = TRUE)The graphs can be obtained as follows:
old.par <- par(no.readonly = TRUE)
on.exit({ layout(1); par(old.par) })
## ----- 1) Make a 2-row layout: top row = legend, bottom row = 2 plots -----
layout(
matrix(c(1,1, 2,3), nrow = 2, byrow = TRUE), # top spans both columns
heights = c(1.2, 8) # adjust top strip height
)
## ----- 2) TOP STRIP: shared model legend (outside the plots) -----
par(mar = c(0, 0, 0, 0))
plot.new()
legend("center", bty = "n", horiz = TRUE, title = "Model",
legend = c("mod.bs", "mod.bs2", "mod.bs3", "mod.bs4"),
lty = 1, lwd = 2, col = c("black","blue","red","green"), cex = 1)
## ----- 3) LEFT PANEL -----
par(mar = c(4, 4, 2, 1))
plot(predict_mod_amr0$results$time.pts, predict_mod_amr0$results$hazard,
type = "l", lwd = 2, xlab = "Time (years)", ylab = "Excess hazard",
ylim = c(0, 0.5))
lines(predict_mod_amr1$results$time.pts, predict_mod_amr1$results$hazard, lwd = 2, lty = 2)
lines(predict_mod2_arm0$results$time.pts, predict_mod2_arm0$results$hazard, lwd = 2, col = "blue", lty = 1)
lines(predict_mod2_arm1$results$time.pts, predict_mod2_arm1$results$hazard, lwd = 2, col = "blue", lty = 2)
lines(predict_mod3_arm0$results$time.pts, predict_mod3_arm0$results$hazard, lwd = 2, col = "red", lty = 1)
lines(predict_mod3_arm1$results$time.pts, predict_mod3_arm1$results$hazard, lwd = 2, col = "red", lty = 2)
lines(predict_mod4_arm0$results$time.pts, predict_mod4_arm0$results$hazard, lwd = 2, col = "green", lty = 1)
lines(predict_mod4_arm1$results$time.pts, predict_mod4_arm1$results$hazard, lwd = 2, col = "green", lty = 2)
grid()
legend("topright", bty = "n", title = "Treatment arm",
legend = c("No immunotherapy","Immunotherapy"),
lty = c(1,2), lwd = 2, col = "black", cex = 0.9)
## ----- 4) RIGHT PANEL -----
par(mar = c(4, 4, 2, 1))
plot(predict_mod_amr0$results$time.pts, predict_mod_amr0$results$surv,
type = "l", lwd = 2, xlab = "Time (years)", ylab = "Net survival",
ylim = c(0, 1))
lines(predict_mod_amr1$results$time.pts, predict_mod_amr1$results$surv, lwd = 2, lty = 2)
lines(predict_mod2_arm0$results$time.pts, predict_mod2_arm0$results$surv, lwd = 2, col = "blue", lty = 1)
lines(predict_mod2_arm1$results$time.pts, predict_mod2_arm1$results$surv, lwd = 2, col = "blue", lty = 2)
lines(predict_mod3_arm0$results$time.pts, predict_mod3_arm0$results$surv, lwd = 2, col = "red", lty = 1)
lines(predict_mod3_arm1$results$time.pts, predict_mod3_arm1$results$surv, lwd = 2, col = "red", lty = 2)
lines(predict_mod4_arm0$results$time.pts, predict_mod4_arm0$results$surv, lwd = 2, col = "green", lty = 1)
lines(predict_mod4_arm1$results$time.pts, predict_mod4_arm1$results$surv, lwd = 2, col = "green", lty = 2)
grid()
legend("topright", bty = "n", title = "Treatment arm",
legend = c("No immunotherapy","Immunotherapy"),
lty = c(1,2), lwd = 2, col = "black", cex = 0.9)
License
GPL 3.0, for academic use.
Acknowledgments
We are grateful to the members of the CENSUR Survival Group for their helpful comments.
Session Info
sessionInfo()
#> R version 4.5.1 (2025-06-13 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26100)
#>
#> Matrix products: default
#> LAPACK version 3.12.1
#>
#> locale:
#> [1] LC_COLLATE=C LC_CTYPE=French_France.utf8
#> [3] LC_MONETARY=French_France.utf8 LC_NUMERIC=C
#> [5] LC_TIME=French_France.utf8
#>
#> time zone: Europe/Paris
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] survexp.fr_1.2 xhaz_2.1.0 mexhaz_2.6 survival_3.8-3
#>
#> loaded via a namespace (and not attached):
#> [1] vctrs_0.6.5 WriteXLS_6.8.0 cli_3.6.5
#> [4] knitr_1.50 rlang_1.1.6 xfun_0.52
#> [7] stringi_1.8.7 jsonlite_2.0.0 glue_1.8.0
#> [10] RcppParallel_5.1.11-1 statmod_1.5.0 htmltools_0.5.8.1
#> [13] sass_0.4.10 rmarkdown_2.29 grid_4.5.1
#> [16] evaluate_1.0.5 jquerylib_0.1.4 MASS_7.3-65
#> [19] fastmap_1.2.0 yaml_2.3.10 lifecycle_1.0.4
#> [22] numDeriv_2016.8-1.1 stringr_1.5.1 compiler_4.5.1
#> [25] Rcpp_1.1.0 optimParallel_1.0-2 lamW_2.2.5
#> [28] rstudioapi_0.17.1 lattice_0.22-7 digest_0.6.37
#> [31] R6_2.6.1 parallel_4.5.1 splines_4.5.1
#> [34] magrittr_2.0.3 bslib_0.9.0 Matrix_1.7-4
#> [37] tools_4.5.1 cachem_1.1.0