TABLE OF CONTENTS
S3Methods/predict.splinesurv [ Functions ]
NAME
predict.splinesurv --- prediction method for splinesurv objects
FUNCTION
Compute curve estimates, or linear predictor/risk estimates for splinesurv.
This routine mainly works recursively. If called with iter=NULL, it calls itself with a subset of iterations multiple times and aggregates the results.
SYNOPSIS
4033 predict.splinesurv <- function(object, type = c("hazard", "survival", "lp", "risk", "frailty"), 4034 x = NULL, newdata = NULL, iter = NULL, fn = mean, alpha = NULL, npost = 100, ...)
INPUTS
See package documentation
SOURCE
4037 { 4038 type <- match.arg(type) 4039 fit <- object; haz <- 1 4040 ntimes <- 100 4041 # check for valid newdata 4042 if((type == "hazard" | type == "survival") & !is.null(newdata)) if(dim(newdata)[1] > 1) 4043 stop("newdata may only have one row") 4044 # if iter=NULL and a curve is required, call the routine again with a subset of iters 4045 if((type == "hazard" | type == "frailty" | type == "survival") & is.null(iter)){ 4046 # get the set of "x" values at which the curve will be predicted 4047 if(type == "hazard" | type == "survival") { 4048 if(is.null(x) | is.character(x)) 4049 x <- seq(from = min(fit$data$time), to = max(fit$data$time), length = ntimes) 4050 } 4051 if(type == "frailty" & is.null(x)) { 4052 knots <- fit$history$frailty.spline.knots[1, ] 4053 knots <- knots[knots>-Inf] 4054 bounds <- range(knots) 4055 x <- seq(from = bounds[1], to = bounds[2], length = ntimes) 4056 } 4057 # get the set of iterations that will be used 4058 ngooditer <- min(fit$control$maxiter - fit$control$burnin, npost) 4059 iters <- round(seq(from = fit$control$burnin + 1, to = fit$control$maxiter, 4060 length = ngooditer)) 4061 # matrix for storage of predictions 4062 preds <- matrix(0, length(x), ngooditer) 4063 i <- 1 4064 # call predict for a single iteration (after this if statement) 4065 for(iter in iters) { 4066 thispred <- predict(fit, type, x, newdata, iter, ...) 4067 preds[, i] <- thispred[, 2] 4068 i <- i + 1 4069 } 4070 # aggregate the set of predictions 4071 pred <- apply(preds, 1, fn, na.rm = TRUE) 4072 if(!is.null(alpha)) quantiles <- t(apply(preds, 1, quantile, 4073 probs = c(alpha / 2, 1 - alpha / 2), na.rm = TRUE)) else quantiles <- NULL 4074 out <- cbind(x, pred, quantiles) 4075 colnames(out) <- c(colnames(thispred), colnames(quantiles)) 4076 return(as.data.frame(out)) 4077 } 4078 4079 # similar to above, but handles "risk" as well 4080 if(type == "hazard" | type == "survival" | (type == "risk" & !is.null(x))){ 4081 if(is.null(x) | is.character(x)) times <- seq(from = min(fit$data$time), 4082 to = max(fit$data$time), length = ntimes) else times <- x 4083 # construct a temp hazard curve 4084 hazard <- fit$hazard 4085 hazard$x <- times; hazard$haspar <- haspar(hazard); 4086 hazard$hasspline <- hasspline(hazard); hazard$spline.norm <- FALSE 4087 if(is.null(iter)){ 4088 # if iter=NULL, call the routine for a subset of iters as above 4089 ngooditer <- min(fit$control$maxiter - fit$control$burnin, npost) 4090 iters <- round(seq(from = fit$control$burnin, to = fit$control$maxiter, 4091 length = ngooditer)) 4092 preds <- matrix(0, length(times), ngooditer) 4093 i <- 1 4094 for(iter in iters) { 4095 thispred <- predict(fit, type, times, newdata, iter, ...) 4096 preds[, i] <- thispred[, 2] 4097 i <- i + 1 4098 } 4099 pred <- apply(preds, 1, fn, na.rm = TRUE) 4100 if(!is.null(alpha)) quantiles <- apply(preds, 1, quantile, 4101 probs = c(alpha / 2, 1 - alpha / 2), na.rm = TRUE) else quantiles <- NULL 4102 out <- cbind(times, pred, t(quantiles)) 4103 colnames(out) <- c(colnames(thispred), rownames(quantiles)) 4104 return(as.data.frame(out)) 4105 }else{ 4106 # continue constructing a temporary hazard RCurve 4107 hazard$spline.par <- fit$history$hazard.spline.par[iter, ] 4108 hazard$spline.knots <- fit$history$hazard.spline.knots[iter, ] 4109 hazard$param.par <- fit$history$hazard.param.par[iter, ] 4110 hazard$spline.par <- hazard$spline.par[hazard$spline.par>-Inf] 4111 hazard$spline.knots <- hazard$spline.knots[hazard$spline.knots>-Inf] 4112 hazard$weight <- fit$history$hazard.weight[iter, ] 4113 } 4114 # set up the basis and evaluate the parametric and spline components 4115 hazard <- makesplinebasis(hazard, quick = TRUE) 4116 hazard <- evalparametric(hazard) 4117 hazard <- evalspline(hazard, quick = TRUE) 4118 haz <- hazard$y 4119 # if newdata is nonzero, compute the frailty and covariate multiplier for each row 4120 if(!is.null(newdata)){ 4121 risk <- predict(fit, "risk", NULL, newdata, iter)$risk[1] 4122 haz <- haz * risk 4123 clusterind <- attr(fit$terms, "special")$cluster 4124 if(!is.null(clusterind)) cluster <- newdata[[clusterind]] else cluster <- NULL 4125 if(!is.null(cluster) && !is.na(cluster)) 4126 frail <- fit$history$frailty[iter, cluster] else frail <- 1 4127 haz <- haz * frail 4128 } 4129 if(type == "hazard") return(data.frame(time = times, hazard = haz)) 4130 # compute survival curve by cumulative summing 4131 if(type == "survival"){ 4132 dx <- c(diff(times), 0) 4133 surv <- exp(-cumsum(haz * dx)) 4134 return(data.frame(time = times, survival = surv)) 4135 } 4136 } 4137 # for lp and risk, evaluate the model terms and do the prediction 4138 if(type == "lp" | type == "risk") 4139 { 4140 if(is.null(newdata)) { 4141 vars <- model.matrix(fit$terms, model.frame(fit))[, -1, drop = FALSE] 4142 }else{ 4143 temp <- cbind(data.frame(i = 0, j = 0, time = 0, delta = 0, newdata)) 4144 vars <- model.matrix(fit$terms, model.frame(fit, data = temp))[, -1, drop = FALSE] 4145 } 4146 if(is.null(iter)) 4147 coef <- fit$posterior.mean$coefficients 4148 else 4149 coef <- fit$history$coefficients[iter, ] 4150 # compute lp 4151 lp <- as.numeric(vars%*%coef) 4152 if(type == "lp") { 4153 lp <- data.frame(lp) 4154 try(rownames(lp) <- if(is.null(newdata)) rownames(fit$data) else 4155 rownames(newdata), silent = TRUE) 4156 return(lp) 4157 } 4158 } 4159 # for risk, given lp and hazard, multiply additionally by the frailty 4160 if(type == "risk") 4161 { 4162 pred <- exp(lp) 4163 if(is.null(newdata)) { 4164 Ji <- table(fit$data$i) 4165 if(is.null(iter)) 4166 frailty <- rep(fit$posterior.mean$frailty, Ji) 4167 else 4168 frailty <- rep(fit$history$frailty[iter, ], Ji) 4169 pred <- frailty * pred 4170 } 4171 risk <- outer(haz, pred, "*") 4172 try(colnames(risk) <- if(is.null(newdata)) rownames(fit$data) else 4173 rownames(newdata), silent = TRUE) 4174 if(any(duplicated(colnames(risk)))) colnames(risk) <- NULL 4175 if(dim(risk)[1] == 1){ 4176 risk <- as.data.frame(t(risk)) 4177 colnames(risk) <- "risk" 4178 }else{ 4179 risk <- data.frame(time = times, risk) 4180 } 4181 return(risk) 4182 } 4183 # for frailty curve 4184 if(type == "frailty") 4185 { 4186 # construct a temp frailty curve that can be evaluated 4187 frailty <- fit$frailty 4188 if(is.null(x)) x <- seq(from = max(0, min(frailty$spline.knots)), 4189 to = max(frailty$spline.knots), length = ntimes) 4190 frailty$x <- x; frailty$haspar <- haspar(frailty); 4191 frailty$hasspline <- hasspline(frailty); frailty$spline.norm <- TRUE 4192 if(is.null(iter)){ 4193 frailty$spline.par <- fit$posterior.mean$frailty.spline.par 4194 frailty$param.par <- fit$posterior.mean$frailty.param.par 4195 frailty$weight <- fit$posterior.mean$frailty.weight 4196 if(type == "frailty") return(data.frame(x = x, density = 1)) 4197 }else{ 4198 frailty$spline.par <- fit$history$frailty.spline.par[iter, ] 4199 frailty$spline.knots <- fit$history$frailty.spline.knots[iter, ] 4200 frailty$param.par <- fit$history$frailty.param.par[iter, ] 4201 frailty$spline.par <- frailty$spline.par[frailty$spline.par>-Inf] 4202 frailty$spline.knots <- frailty$spline.knots[frailty$spline.knots>-Inf] 4203 frailty$weight <- fit$history$frailty.weight[iter, ] 4204 } 4205 frailty <- makesplinebasis(frailty, quick = TRUE) 4206 frailty <- evalparametric(frailty) 4207 frailty <- evalspline(frailty, quick = TRUE) 4208 density <- frailty$y 4209 return(data.frame(x = x, density = density)) 4210 } 4211 }