random effects Cox model

I am trying to plot the dose-response relation between pm2.5 exposure and hypertension incidence. I fitted a random effects cox model, in which a random effect was added for the study centers. Then I used the plotHR function to plot the dose-response relation. But I met an error. Below is the example R code I used.

library(survival)
library(coxme) 
library(splines)
library(Greg)

data("eortc")
eortc$age<-rnorm(2323,40,10)

efit1 <- coxph(Surv(y, uncens) ~ ns(age) + trt , eortc)
plotHR(efit1,term = 1,xlim = c(30,70))

efit2 <- coxme(Surv(y, uncens) ~ ns(age) + trt + (1|center), eortc)
plotHR(efit2,term = 1,xlim = c(30,70))

I can plot efit1 using plotHR, but I met an error when plot efit2 in which I added a random effect in the cox model. Anyone knows how to solve this problem? Thanks!

This show what predict.coxme can return for such a model. I'm guessing it won't be totally satisfying, but I claim that it is useful because it demonstrates what mixed models are like "underneath the hood". Each center has a separate prediction but since there is no claim that these have any particular distribution, efforts to aggregate them at their means is considered statistically suspect. This plots the exponentiated linear predictors for the model involving a "spline" fit with no knots (which gives you a bunch of lines). I'm coloring by trt status: black for "0" and red for "1"

plot( eortc$age, exp( predict(efit2)), col=1+(eortc$trt==1) )

The "average" difference between the trt==0 and trt==1 groups does show up and is consistent with the measured treatment effect of exp(0.705) -> [1] 2.023847, and the "age" effect which in the model was not significant is a very shallow linear rise.