Return covariance matrix for the predictions using `predict.coxph`

[darentsai]

Dear maintainers,

The predict method of coxph object is able to produce standard errors for the predictied values if the se.fit option is TRUE. I further want the covariance matrix for the predictions.

I try to replace the last return part of predict.coxph function

https://github.com/therneau/survival/blob/e8de7b732daa9a82c1ceaad93fd68ee7545c7736/R/predict.coxph.R#L333-L334

with

https://github.com/darentsai/survival/blob/2d3467e28d4f0e32049d1c972c5ac437fdca46a5/R/predict.coxph.R#L333-L336

    if (se.fit) {
      if (full.cov) list(fit=pred, vcov.fit = newx %*% object$var %*% t(newx))
      else list(fit=pred, se.fit=se)
    } else pred

and the output includes a plausible covariance matrix (element vcov.fit). I have checked that the square root of the diagonal of vcov.fit is exactly se.fit.

Is it feasible for predict.coxph to support the covariance matrix in the future?

Best regards,

[therneau]

The number of predictions is potentially very large, and returning an n x n matrix could easily use up all the memory. It would have to be an option, would only apply to linear predictors (not risk, etc). What is the use case for it?