Handle non-standard factor contrasts

DESCRIPTION

@@ -70,4 +70,4 @@ BugReports: https://github.com/chjackson/flexsurv/issues
 VignetteBuilder: knitr
 LazyData: yes
 LinkingTo: Rcpp
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.0

NEWS.md

@@ -1,4 +1,4 @@
-# Version 2.3 (2023-??-??)
+# Version 2.3 (2024-??-??)
 
 * Analytic Hessian calculation for models where this is possible, that
   is, Weibull, Gompertz, exponential, and spline models in hazard and
@@ -10,6 +10,8 @@
   consistently with other models, as a weighted sum of individual
   log-likelihoods.
 
+* Non-default factor contrasts now handled.
+
 * `pmatrix.simfs` can now accept a vector of times `t` and has a
   `tidy` output option.
 

R/flexsurvreg.R

@@ -866,11 +866,11 @@ flexsurvreg <- function(formula, anc=NULL, data, weights, bhazard, rtrunc, subse
     temp[[1]] <- as.name("model.frame")
 
     f2 <- concat.formulae(formula,forms, data)
-    temp[["formula"]] <- f2
+    temp[["formula"]] <- terms(f2)
     if (missing(data)) temp[["data"]] <- environment(formula)
     m <- eval(temp, parent.frame())
 
-    m <- droplevels(m) # remove unused factor levels after subset applied
+    m <- droplevels_keepcontrasts(m) # remove unused factor levels after subset applied
     attr(m,"covnames") <- attr(f2, "covnames") # for "newdata" in summary
     attr(m,"covnames.orig") <- intersect(colnames(m), attr(f2, "covnames.orig")) # for finding factors in plot method
     Y <- check.flexsurv.response(model.extract(m, "response"))
@@ -1390,3 +1390,22 @@ deriv_supported <- function(Y){
   left_cens <- is.finite(Y[!event, "time2"])
   !any(left_cens)
 }
+
+
+droplevels_factor_keepcontrasts <- function (x, exclude = if (anyNA(levels(x))) NULL else NA, ...) {
+  ct <- attr(x, "contrasts")
+  x <- factor(x, exclude = exclude)
+  contrasts(x) <- ct
+  x
+}
+
+droplevels_keepcontrasts <- function (x, except = NULL, exclude, ...) 
+{
+    ix <- vapply(x, is.factor, NA)
+    if (!is.null(except)) 
+        ix[except] <- FALSE
+    x[ix] <- if (missing(exclude)) 
+        lapply(x[ix], droplevels_factor_keepcontrasts)
+    else lapply(x[ix], droplevels_factor_keepcontrasts, exclude = exclude)
+    x
+}

inst/NEWS

@@ -1,6 +1,6 @@
 -*- text -*-
 
-Version 2.3 (2023-??-??)
+Version 2.3 (2024-??-??)
 -------------------------
 
 * Analytic Hessian calculation for models where this is possible, that
@@ -13,6 +13,8 @@ Version 2.3 (2023-??-??)
   consistently with other models, as a weighted sum of individual
   log-likelihoods.
 
+* Non-default factor contrasts now handled.
+
 * `pmatrix.simfs` can now accept a vector of times `t` and has a
   `tidy` output option.
 

tests/testthat/test_contrasts.R

@@ -0,0 +1,22 @@
+## Thanks to @https://github.com/anddis
+## Fixed https://github.com/chjackson/flexsurv/issues/178
+
+test_that("Non-default factor contrasts", { 
+  veteran$celltype <- factor(veteran$celltype)
+  fit.noc <- flexsurvreg(Surv(time, status) ~ celltype,
+                         data = veteran, 
+                         dist = "exp")
+
+  veteran$celltype2 <- veteran$celltype
+  contrasts(veteran$celltype2) <- stats::contr.helmert(4)
+  fit.c <- flexsurvreg(Surv(time, status) ~ celltype2,
+                       data = veteran, 
+                       dist = "exp")
+
+  fit.poi <- glm(status ~ celltype2 + offset(log(time)),
+                 data = veteran,
+                 family = "poisson")
+
+  expect_equivalent(coef(fit.c), coef(fit.poi))
+  expect_true(coef(fit.c)[1] != coef(fit.noc)[1])
+})