Since package gbm requires package survival and can use the Surv() function in the canonical formula interface, and there is also a reference to survfit() in the basehaz.gbm() documentation in the gbm vignette, it might be possible to pull out what you're looking for using the gbm model object. I doubt it, though, based on what I've read so far.

So you might have to go to the source code and make your own function to extract or reconstruct what you need to mimic plot(survfit()). See also the (hidden) helper function reconstructGBMdata().

(I know this isn't the most helpful answer, BTW, but I don't have enough Cross-Validated user history to comment rather than answer.)

I also wanted to comment on using gradient boosting models, specifically for handling missing values. The short answer here is that the gbm algorithm handles missing values explicitly, obviating the need for user-handled imputation. I suggest looking up package rpart if you want to understand the technical details better.

Since package gbm requires package survival and can use the Surv() function in the canonical formula interface, and there is also a reference to survfit() in the basehaz.gbm() documentation in the gbm vignette, it might be possible to pull out what you're looking for using the gbm model object. I doubt it, though, based on what I've read so far.

So you might have to go to the source code and make your own function to extract or reconstruct what you need to mimic plot(survfit()). See also the (hidden) helper function reconstructGBMdata().

I also wanted to comment on using gradient boosting models, specifically for handling missing values. The short answer here is that the gbm algorithm handles missing values explicitly, obviating the need for user-handled imputation. I suggest looking up package rpart if you want to understand the technical details better.