Add predictEventProb function for pamm and list objects

[jc-espinosa]

This update introduces two new S3 methods predictEventProb.pamm and predictEventProb.list to the predict.R script. These additions enable compatibility with the pec package when working in a competing risks setting using pammtools.

[fabian-s]

Review of PR #252 — predictEventProb for pammtools

Overview

The PR adds predictEventProb.pamm() and predictEventProb.list() for pec package compatibility in competing risks. Valuable goal, but the implementation has several bugs and design issues that need fixing before merge.

Bugs (must fix)

1. sp_all_cause - 1e-5 produces negative CIF values

The existing get_cif.default uses - 1e-20. The PR uses - 1e-5, which is 15 orders of magnitude larger. When overall survival drops below 1e-5 (right tail), CIF increments go negative — a CIF is a probability and must be non-decreasing. Fix: use - 1e-20 or .Machine$double.eps to match existing code.

2. object[["attr_ped"]]$risks is likely NULL

In predictEventProb.pamm(), all_causes is extracted from object[["attr_ped"]]$risks. But looking at get_cif.default (line ~776 of add-functions.R), this exact approach was abandoned in favor of levels(newdata[[cause_var]]) because attr_ped on the model object isn't reliable. If all_causes is NULL, the for loop never executes, no csh* columns are created, and the final pull() silently fails.

3. Hardcoded defaults in predictEventProb.list

predictEventProb.list <- function(
    object=pam_csh,        # ← debugging leftover
    newdata=test_fourD,    # ← debugging leftover
    times=times_eval,      # ← debugging leftover

These are development objects from the author's session. If defaults are ever evaluated, R throws a confusing "object not found" error. The signature must match the generic: function(object, newdata, times, cause, ...).

4. id_var undefined when newdata is already PED

id_var is assigned only inside the if (!is.ped(...)) block but used unconditionally later in arrange(), group_by(), and the final map(). Passing PED input will error with "object 'id_var' not found".

5. Bare arrange(id, times) — uses unquoted id

The code uses arrange(id, times) instead of arrange(.data[[id_var]], .data$times). This (a) hardcodes the column name id instead of using id_var, and (b) triggers an R CMD check NOTE for no visible binding.

Design Issues (should fix)

6. Dispatching on plain list is too broad

Any list passed to predictEventProb() would match — data frames, gam objects, etc. The principled fix is a custom class (e.g., pamm_list) with a constructor that validates all elements are pamm objects with compatible cut points. At minimum, add validation at the top of the method.

7. get() inside dplyr is non-idiomatic

lapply(cause_vars, function(var) cumsum(get(var) * intlen))

get() bypasses the dplyr data mask protocol. Replace with .data[[var]] or rewrite using rowSums(across(...)):

mutate(
  h_all = rowSums(across(all_of(cause_vars))),
  sp_all_cause = exp(-cumsum(h_all * intlen))
)

8. predict() without unpam()

The existing predictSurvProb.pamm uses predict(unpam(object), ...). The PR calls predict(object, ...) directly. It works today because no predict.pamm() exists, but it's inconsistent and fragile if one is added later.

9. @importFrom pammtools get_intervals — a package cannot import from itself. Remove.

10. ~80% code duplication

The PED preparation and CIF calculation are nearly identical between the two methods. Extract shared logic into internal helpers.

11. Missing newline at end of file — trivial but will trigger a NOTE.

Mathematical Note

The CIF formula uses S(t_j) (survival at end of interval) instead of S(t_{j-1}) (survival entering the interval) as in the standard Aalen-Johansen estimator. However, this matches the existing get_cif.default implementation, so it's internally consistent. For fine time grids the bias is small. Worth noting but not a blocker for this PR specifically.

Missing Tests

No tests are included. At minimum:

• CIF values in [0, 1], non-decreasing
• Sum of cause-specific CIFs + S(t) ≈ 1
• Output matrix dimensions: nrow(newdata) × length(times)
• .pamm and .list methods produce equivalent results
• Error on times > max(brks)
• Error on invalid/non-existent cause
• pec::pec() integration smoke test

Interaction with Other Open PRs

• PR Remove internal columns from make_newdata output #269 (Remove internal columns from make_newdata output): Removes intlen from make_newdata() output and adds reconstruct_intlen(). PR Add predictEventProb function for pamm and list objects #252's code builds its own PED grid via combine_df() so the direct path isn't broken, but if anyone passes make_newdata() output as newdata to predictEventProb, there's no intlen and no call to reconstruct_intlen(). After Remove internal columns from make_newdata output #269 merges, the is.ped() branch would need reconstruct_intlen().

• PR Extend arbitrary time point support in add_cif #271 (Extend arbitrary time point support in add_cif): Extends the native add_cif() for arbitrary time points — the same domain as this PR. If both merge, there would be two independent CIF implementations that could diverge. Consider whether predictEventProb should delegate to add_cif internally rather than reimplementing the calculation.

• PR Stabilize cumulative/simulation edge cases and gg_partial ci behavior #268 (Stabilize cumulative/simulation edge cases): Addresses numerical stability in cumulative calculations — directly relevant to the 1e-5 vs 1e-20 issue here.

Verdict

The feature is wanted, but the implementation needs significant rework before merge. Suggested next steps:

1. Fix the 5 bugs above
2. Extract shared helpers and add a pamm_list class
3. Add tests
4. Coordinate with PRs Remove internal columns from make_newdata output #269/Extend arbitrary time point support in add_cif #271 to avoid duplication of CIF logic