model_parameters struggles with time transform coxph.penal model

Fixes #286
easystats · Nov 16, 2024 · dfde2d3 · dfde2d3
1 parent 2636aff
commit dfde2d3
Show file tree

Hide file tree

Showing 7 changed files with 61 additions and 15 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: parameters
 Title: Processing of Model Parameters
-Version: 0.23.0.8
+Version: 0.23.0.9
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",

diff --git a/NEWS.md b/NEWS.md
@@ -5,6 +5,10 @@
 * The `robust` argument, which was deprecated for a long time, is now no longer
   supported. Please use `vcov` and `vcov_args` instead.
 
+## Changes
+
+* Added support for `coxph.panel` models.
+
 ## Bug fixes
 
 * Fixed bug when extracting 'pretty labels' for model parameters, which could

diff --git a/R/ci_generic.R b/R/ci_generic.R
@@ -3,12 +3,8 @@
                         ci = 0.95,
                         method = "wald",
                         dof = NULL,
-                        effects = c("fixed", "random", "all"),
-                        component = c(
-                          "all", "conditional", "zi", "zero_inflated",
-                          "dispersion", "precision", "scale",
-                          "smooth_terms", "full", "marginal"
-                        ),
+                        effects = "fixed",
+                        component = "all",
                         vcov = NULL,
                         vcov_args = NULL,
                         verbose = TRUE,
@@ -26,8 +22,14 @@
     )
   )
 
-  effects <- match.arg(effects)
-  component <- match.arg(component)
+  effects <- insight::validate_argument(effects, c("fixed", "random", "all"))
+  component <- insight::validate_argument(
+    component,
+    c(
+      "all", "conditional", "zi", "zero_inflated", "dispersion", "precision",
+      "scale", "smooth_terms", "full", "marginal"
+    )
+  )
 
   if (method == "ml1") { # nolint
     return(ci_ml1(model, ci = ci))
@@ -90,7 +92,7 @@
   }
 
   # check if all estimates are non-NA
-  params <- .check_rank_deficiency(params, verbose = FALSE)
+  params <- .check_rank_deficiency(model, params, verbose = FALSE)
   # for polr, we need to fix parameter names
   params$Parameter <- gsub("Intercept: ", "", params$Parameter, fixed = TRUE)
 
@@ -176,7 +178,9 @@
   if ("Effects" %in% names(params) && effects != "fixed") out$Effects <- params$Effects
   if ("Response" %in% names(params)) out$Response <- params$Response
 
-  if (anyNA(params$Estimate)) {
+  # for cox-panel models, we have non-linear parameters with NA coefficient,
+  # but test statistic and p-value - don't check for NA estimates in this case
+  if (anyNA(params$Estimate) && !inherits(model, "coxph.penal")) {
     out[stats::complete.cases(out), ]
   } else {
     out

diff --git a/R/extract_parameters.R b/R/extract_parameters.R
@@ -58,7 +58,7 @@
   statistic <- insight::get_statistic(model, component = component)
 
   # check if all estimates are non-NA
-  parameters <- .check_rank_deficiency(parameters)
+  parameters <- .check_rank_deficiency(model, parameters)
 
 
   # ==== check if we really have a component column
@@ -437,7 +437,7 @@
   statistic <- insight::get_statistic(model, component = "all")
 
   # check if all estimates are non-NA
-  parameters <- .check_rank_deficiency(parameters)
+  parameters <- .check_rank_deficiency(model, parameters)
 
   # sometimes, due to merge(), row-order messes up, so we save this here
   original_order <- parameters$.id <- seq_len(nrow(parameters))
@@ -1014,7 +1014,12 @@
 # tools -------------------------
 
 
-.check_rank_deficiency <- function(p, verbose = TRUE) {
+.check_rank_deficiency <- function(model, p, verbose = TRUE) {
+  # for cox-panel models, we have non-linear parameters with NA coefficient,
+  # but test statistic and p-value - don't check for NA estimates in this case
+  if (!is.null(model) && inherits(model, "coxph.penal")) {
+    return(p)
+  }
   if (anyNA(p$Estimate)) {
     if (isTRUE(verbose)) {
       insight::format_alert(

diff --git a/R/p_value_kenward.R b/R/p_value_kenward.R
@@ -75,7 +75,7 @@ p_value_kenward.lmerMod <- function(model, dof = NULL) {
   params <- insight::get_parameters(model, component = component)
 
   # check if all estimates are non-NA
-  params <- .check_rank_deficiency(params, verbose = FALSE)
+  params <- .check_rank_deficiency(model, params, verbose = FALSE)
 
   if (is.null(statistic)) {
     statistic <- insight::get_statistic(model, component = component)

diff --git a/tests/testthat/_snaps/coxph.md b/tests/testthat/_snaps/coxph.md
@@ -0,0 +1,16 @@
+# model_parameters coxph-panel
+
+    Code
+      print(model_parameters(mod))
+    Output
+      Parameter         | Coefficient |       SE |        95% CI | Chi2(1) |      p
+      -----------------------------------------------------------------------------
+      ph ecog [ok]      |        0.36 |     0.20 | [-0.03, 0.75] |    3.19 | 0.074 
+      ph ecog [limited] |        0.87 |     0.23 | [ 0.41, 1.33] |   13.87 | < .001
+      age, linear       |        0.01 | 9.36e-03 | [-0.01, 0.03] |    1.30 | 0.253 
+      age, nonlin       |             |          |               |    2.83 | 0.093 
+    Message
+      
+      Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
+        using a Wald z-distribution approximation.
+
diff --git a/tests/testthat/test-coxph.R b/tests/testthat/test-coxph.R
@@ -62,3 +62,20 @@ test_that("model_parameters", {
   unloadNamespace("survey")
   unloadNamespace("survival")
 })
+
+
+skip_if_not_installed("withr")
+
+withr::with_package(
+  "survival",
+  test_that("model_parameters coxph-panel", {
+    set.seed(123)
+    # a time transform model
+    mod <- survival::coxph(
+      survival::Surv(time, status) ~ ph.ecog + tt(age),
+      data = lung,
+      tt = function(x, t, ...) pspline(x + t / 365.25)
+    )
+    expect_snapshot(print(model_parameters(mod)))
+  })
+)