I wrote these functions to calculate survival probabilities at various levels of the random effect estimate from a coxme.object. The code is heavily adapted from survfit.coxph(). It should work, but all errors and inelegant hacks I have introduced are certainly my own doing. Comments/improvements welcome, enjoy!

survProbs.coxme.R

#Internal Functions
MYagsurv <- function(y, x, wt, risk, survtype=3, vartype=3) {
    nvar <- ncol(as.matrix(x))
    status <- y[, ncol(y)]
    dtime <- y[, ncol(y) - 1]
    death <- (status == 1)
    time <- sort(unique(dtime))
    nevent <- as.vector(rowsum(wt * death, dtime))
    ncens <- as.vector(rowsum(wt * (!death), dtime))
    wrisk <- c(wt * risk)  #Had to add c() to remove the dimnames so that the multiplication later would work
    
    rcumsum <- function(x) rev(cumsum(rev(x)))
    nrisk <- rcumsum(rowsum(wrisk, dtime))
    irisk <- rcumsum(rowsum(wt, dtime))
    
    if (ncol(y) == 2) {
        temp2 <- rowsum(wrisk * x, dtime)
        xsum <- apply(temp2, 2, rcumsum)
    } else {
        delta <- min(diff(time))/2
        etime <- c(sort(unique(y[, 1])), max(y[, 1]) + delta)
        indx <- approx(etime, 1:length(etime), time, method = "constant", 
            rule = 2, f = 1)$y
        esum <- rcumsum(rowsum(wrisk, y[, 1]))
        nrisk <- nrisk - c(esum, 0)[indx]
        irisk <- irisk - c(rcumsum(rowsum(wt, y[, 1])), 0)[indx]
        xout <- apply(rowsum(wrisk * x, y[, 1]), 2, rcumsum)
        xin <- apply(rowsum(wrisk * x, dtime), 2, rcumsum)
        xsum <- xin - (rbind(xout, 0))[indx, , drop = F]
    }

    #Calculate number of deaths and number of times
    ndeath <- rowsum(status, dtime)
    ntime <- length(time)
    
    #Final calcs, send all of this stuff to the proper C subroutine
    if (survtype == 1) {
        indx <- (which(status == 1))[order(dtime[status == 1])]
        km <- .C(survival:::Cagsurv4, as.integer(ndeath), as.double(risk[indx]), 
            as.double(wt[indx]), as.integer(ntime), as.double(nrisk), 
            inc = double(ntime))
    }
    
    if (survtype == 3 || vartype == 3) {
      xsum2 <- rowsum((wrisk * death) * x, dtime)
      erisk <- rowsum(wrisk * death, dtime)
  
      tsum <- .C(survival:::Cagsurv5, as.integer(length(nevent)), as.integer(nvar), 
                 as.integer(ndeath), as.double(nrisk), as.double(erisk), 
                 as.double(xsum), as.double(xsum2), sum1 = double(length(nevent)), 
                 sum2 = double(length(nevent)), xbar = matrix(0, 
                 length(nevent), nvar))
    }
    
    #Calc hazard rates, etc and wrap up results
    haz <- switch(survtype, nevent/nrisk, nevent/nrisk, nevent * 
        tsum$sum1)
    varhaz <- switch(vartype, nevent/(nrisk * ifelse(nevent >= 
        nrisk, nrisk, nrisk - nevent)), nevent/nrisk^2, nevent * 
        tsum$sum2)
    xbar <- switch(vartype, (xsum/nrisk) * haz, (xsum/nrisk) * 
        haz, nevent * tsum$xbar)

    result <- list(n = nrow(y), time = time, n.event = nevent, 
        n.risk = irisk, n.censor = ncens, hazard = haz, cumhaz = cumsum(haz), 
        varhaz = varhaz, ndeath = ndeath, xbar = apply(matrix(xbar, 
        ncol = nvar), 2, cumsum))
    
    if (survtype == 1) {
      result$surv <- cumprod(km$inc)
    }
    else {
      result$surv <- exp(-result$cumhaz)
    }
    
    return(result)
}

MYrisk <- function(object, newdata, rfx=0) {
  dMat <- survival:::model.matrix.coxph(object, newdata)
  xmeans <- object$means
    names(xmeans) <- names(fixef(object))
  dMat <- scale(dMat, center=xmeans, scale=FALSE)
  newrisk <- exp(dMat %*% fixef(object) + rfx)
  return(newrisk)
}

predSurv <- function(survlist, newrisk) {
  n.survlist <- length(survlist)
  strata.names <- names(survlist) 
  n.newrisk <- nrow(newrisk)
  
  for (i in 1:n.survlist) {
    survlist[[i]]$surv <- as.data.frame(outer(survlist[[i]]$surv, newrisk, "^"))
      colnames(survlist[[i]]$surv) <- paste("Surv_newdata_row", 1:n.newrisk, sep="")
    survlist[[i]]$cumhaz <- as.data.frame(outer(survlist[[i]]$cumhaz, newrisk, "*"))
      colnames(survlist[[i]]$cumhaz) <- paste("CumHaz_newdata_row", 1:n.newrisk, sep="")
    survlist[[i]]$newrisk <- newrisk
      dimnames(survlist[[i]]$newrisk)[[1]] <- list(paste("newdata_row", 1:n.newrisk, sep=""))
  }
  
  return(survlist)
}


survProbs <- function(survHat, newtime) {
  out <- NULL
  for (i in 1:length(survHat)) {
    if (any(newtime == survHat[[i]]$time)) {
      idx <- survHat[[i]]$time==newtime
      out[[i]] <- data.frame(
        Time=survHat[[i]]$time[idx],
        `Survival Probability`=survHat[[i]]$surv[idx,1],
        `Cumulative Incidence`=1-survHat[[i]]$surv[idx,1]
      )
    } else{
      btwnrows <- c(max(survHat[[i]]$time[survHat[[i]]$time<newtime]), 
                    min(survHat[[i]]$time[survHat[[i]]$time>newtime]))
      idx <- survHat[[i]]$time %in% btwnrows
      interpDF <- data.frame(
        time=survHat[[i]]$time[idx],
        surv=survHat[[i]]$surv[idx,1]
      )
      out[[i]] <- data.frame(
        Time=newtime,
        `Survival Probability`=approx(x=interpDF[,"time"], y=interpDF[,"surv"], xout=newtime)$y
      )
      out[[i]]$`Cumulative Incidence`=1-out[[i]][,2]
    }
  }
  return(out)
}

#Final Function
survProbs.coxme <- function(object, stratVar=NULL, rfx=0, groupFctr=NULL, level=NULL, newdata, newtime) {
  #prep values for MYagsurv()
  y <- object$y
  x <- object$x
  if (is.null(object$weights)) {
    wt <- rep(1, nrow(y))
  } else wt <- object$weights
  survtype <- vartype <- switch(object$ties, efron=3, breslow=2)
  risk <- exp(x %*% fixef(object))
  #do MYagsurv() by strata
  if (is.null(stratVar)) {
    stratvar <- wt
  }
  if (is.factor(stratVar)) {
     ustrata <- levels(stratVar)
  } else ustrata <- sort(unique(stratVar))
  nstrata <- length(ustrata)
  survlist <- vector("list", nstrata)
  names(survlist) <- ustrata
  for (i in 1:nstrata) {
    indx <- which(stratVar == ustrata[i])
    survlist[[i]] <- MYagsurv(y[indx, , drop = F], x[indx, , drop = F], wt[indx], risk[indx], survtype, vartype)
  }
  #Calc newrisk for user-supplied newdata
  if (!is.null(level)) {
    if (!is.null(groupFctr)) {
      rfx <- as.numeric(ranef(object)[[groupFctr]][level])
    } else rfx <- as.numeric(ranef(object)[[1]][level])
  }
  newrisk <- MYrisk(object, newdata, rfx)
  #Calculate survival function based on covariate profile in newdata
  survHat <- predSurv(survlist, newrisk)
  #Find value of survival functionat newtime
  survProbs <- survProbs(survHat, newtime)
    names(survProbs) <- names(survlist)
  print(survProbs)
  return(list(survList=survlist, survListHat=survHat, survProbs=survProbs))
}