Add method for coxphf models

[dchiu911]

Feature

Can we add a tidy.coxphf() method for coxphf::coxphf() Firth penalized likelihood Cox regression models? In the reprex below we get the error Error in co[, -2, drop = FALSE]: incorrect number of dimensions but the tidied information can be easily extracted:

library(survival)
library(coxphf)
library(broom)
library(tibble)

data(breast)

fit1 <- coxph(Surv(TIME, CENS) ~ `T` + N + G + CD, breast)
#> Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
#> Loglik converged before variable 3 ; coefficient may be infinite.
tidy(fit1)
#> # A tibble: 4 × 5
#>   term  estimate std.error statistic p.value
#>   <chr>    <dbl>     <dbl>     <dbl>   <dbl>
#> 1 T        1.28      0.502   2.55     0.0109
#> 2 N        0.946     0.425   2.23     0.0260
#> 3 G       18.4    4317.      0.00427  0.997 
#> 4 CD       0.400     0.443   0.902    0.367

fit2 <- coxphf(Surv(TIME, CENS) ~ `T` + N + G + CD, breast)
tidy(fit2)
#> coxphf(formula = Surv(TIME, CENS) ~ T + N + G + CD, data = breast)
#> 
#> Model fitted by Penalized ML
#> Confidence intervals and p-values by Profile Likelihood 
#> 
#>         coef  se(coef) exp(coef) lower 0.95  upper 0.95    Chisq           p
#> T  1.2244388 0.4916044  3.402256  1.3627461    9.472184 6.983773 0.008225204
#> N  0.9188882 0.4225734  2.506502  1.1204552    5.832863 5.004409 0.025282830
#> G  2.4244141 1.4735463 11.295610  1.4656675 1451.945938 6.090654 0.013589876
#> CD 0.3971181 0.4418554  1.487532  0.6268672    3.511784 0.822321 0.364502440
#> 
#> Likelihood ratio test=35.95142 on 4 df, p=2.961054e-07, n=100
#> Wald test = 23.20007 on 4 df, p = 0.000115489
#> 
#> Covariance-Matrix:
#>              T           N           G          CD
#> T   0.24167491 -0.00079123 -0.06806805 -0.07374037
#> N  -0.00079123  0.17856824 -0.04416386 -0.05117879
#> G  -0.06806805 -0.04416386  2.17133881 -0.02606956
#> CD -0.07374037 -0.05117879 -0.02606956  0.19523619
#> Error in co[, -2, drop = FALSE]: incorrect number of dimensions

tibble(
  term = names(coef(fit2)),
  estimate = coef(fit2),
  std.error = sqrt(diag(fit2$var)),
  statistic = qchisq(1 - fit2$prob, 1),
  p.value = fit2$prob
)
#> # A tibble: 4 × 5
#>   term  estimate std.error statistic p.value
#>   <chr>    <dbl>     <dbl>     <dbl>   <dbl>
#> 1 T        1.22      0.492     6.98  0.00823
#> 2 N        0.919     0.423     5.00  0.0253 
#> 3 G        2.42      1.47      6.09  0.0136 
#> 4 CD       0.397     0.442     0.822 0.365

^{Created on 2025-02-27 with reprex v2.1.1}