predict.flexrsurv: Predictions for a relative survival model

Description

Predict linear predictors, hazard and cumulative hazard for a model fitted by flexrsuv

Usage

# S3 method for flexrsurv
predict(object, newdata= NULL,
   type = c("lp", "link", "risk", "hazard", "hazardrate",
             "rate", "loghazard", "log", "lograte",
             "cumulative.rate", "cumulative.hazard", "cumulative", "cum",
             "survival", "surv", "netsurv"),
   se.fit=FALSE,  		
   ci.fit = FALSE,
   level = .95,
   na.action=na.pass, ...)

Value

a vector or a list containing the predictions (element "fit") and their standard errors (element "se.fit") if the se.fit option is TRUE.

Arguments

object: the results of a flexrsurv fit.
newdata: Optional new data at which to do predictions. If absent predictions are for the data frame used in the original fit.
type: the type of predicted value. Choices are the linear predictor ("lp", "log", "loghazard", "lograte"), the hazard ("rate", "hazard", "hazardrate", "risk") or the cumulative hazard ("cum", "cumulative.hazard", "cumulative").
se.fit: if TRUE, pointwise standard errors are produced for the predictions (not available for cumulative hazard).
ci.fit: if TRUE, confidence intervale are produced for the predictions.
level: Confidence level.
na.action: function determining what should be done with missing values in newdata. The default is to predict NA.
...: For future methods

Details

For cumulative hazard, the cumulative hazard is computed from 0 until the given end time. The cumulative hazard is computed using the same numerical integration method as the one used to fit the model.

Examples

Run this code


# \donttest{

if (requireNamespace("relsurv", quietly = TRUE)) {

	# data from package relsurv
	data(rdata, package="relsurv")
	
	# rate table from package relsurv
	data(slopop, package="relsurv")
	
	
	# get the death rate at event (or end of followup) from slopop for rdata
	rdata$iage <- findInterval(rdata$age*365.24+rdata$time, attr(slopop, "cutpoints")[[1]])
	rdata$iyear <- findInterval(rdata$year+rdata$time, attr(slopop, "cutpoints")[[2]])
	therate <- rep(-1, dim(rdata)[1])
	for( i in 1:dim(rdata)[1]){
	  therate[i] <- slopop[rdata$iage[i], rdata$iyear[i], rdata$sex[i]]
	}
	
	rdata$slorate <- therate
	
	# change sex coding
	rdata$sex01 <- rdata$sex -1
	# centering age
	rdata$agec <- rdata$age- 60
	
	# fit a relative survival model with a non linear effect of age
	fit <- flexrsurv(Surv(time,cens)~sex01+NLL(age, Knots=60, Degree=3,
	                                           Boundary.knots = c(24, 95)), 
	                 rate=slorate, data=rdata,
	                 knots.Bh=1850,  # one interior knot at 5 years
	                 degree.Bh=3,
	                 Spline = "b-spline", 
	                 initbyglm=TRUE, 
	                 int_meth= "BOOLE",
	                 step=50
	                 )
	summary(fit, correlation=TRUE)
	
	
	
	newrdata <- rdata
	newrdata$age <- rep(60, length(rdata$age))
	newrdata$sex <- factor(newrdata$sex, labels=c("m", "f"))
	
	linpred <- predict(fit, newdata=newrdata, type="lp", ci.fit=TRUE)
	predhazard <- predict(fit, newdata=newrdata, type="hazard" , ci.fit=TRUE)
	predcumhazard <- predict(fit, newdata=newrdata, type="cum", ci.fit=TRUE)
	
	require(ggplot2)
	tmp <- cbind(newrdata, linpred)
	glp <- ggplot(tmp, aes(time, colour=sex))
	glp + geom_ribbon(aes(ymin = lwr, ymax = upr, fill=sex)) +
	   geom_line(aes(y=fit)) +
	   scale_fill_manual(values = alpha(c("blue", "red"), .3))
	
	
	tmp <- cbind(newrdata, predhazard)
	glp <- ggplot(tmp, aes(time, colour=sex))
	glp + geom_ribbon(aes(ymin = lwr, ymax = upr, fill=sex)) +
	   geom_line(aes(y=fit))  +
	   scale_fill_manual(values = alpha(c("blue", "red"), .3))
	
	
	tmp <- cbind(newrdata, predcumhazard)
	glp <- ggplot(tmp, aes(time, colour=sex))
	glp + geom_ribbon(aes(ymin = lwr, ymax = upr, fill=sex)) +
	   geom_line(aes(y=fit)) +
	   scale_fill_manual(values = alpha(c("blue", "red"), .3))
}
# }