psm
Parametric Survival Model
psm
is a modification of Therneau's survreg
function for
fitting the accelerated failure time family of parametric survival
models. psm
uses the rms
class for automatic
anova
, fastbw
, calibrate
, validate
, and
other functions. Hazard.psm
, Survival.psm
,
Quantile.psm
, and Mean.psm
create S functions that
evaluate the hazard, survival, quantile, and mean (expected value)
functions analytically, as functions of time or probabilities and the
linear predictor values.
For the print
method, format of output is controlled by the
user previously running options(prType="lang")
where
lang
is "plain"
(the default), "latex"
, or
"html"
.
The residuals.psm
function exists mainly to compute normalized
(standardized) residuals and to censor them (i.e., return them as
Surv
objects) just as the original failure time variable was
censored. These residuals are useful for checking the underlying
distributional assumption (see the examples). To get these residuals,
the fit must have specified y=TRUE
. A lines
method for these
residuals automatically draws a curve with the assumed standardized
survival distribution. A survplot
method runs the standardized
censored residuals through npsurv
to get Kaplan-Meier estimates,
with optional stratification (automatically grouping a continuous
variable into quantiles) and then through survplot.npsurv
to plot
them. Then lines
is invoked to show the theoretical curve. Other
types of residuals are computed by residuals
using
residuals.survreg
.
Usage
psm(formula=formula(data),
data=parent.frame(), weights,
subset, na.action=na.delete, dist="weibull",
init=NULL, scale=0,
control=survreg.control(),
parms=NULL,
model=FALSE, x=FALSE, y=TRUE, time.inc, …)# S3 method for psm
print(x, correlation=FALSE, digits=4, coefs=TRUE,
title, …)
Hazard(object, …)
# S3 method for psm
Hazard(object, …) # for psm fit
# E.g. lambda <- Hazard(fit)
Survival(object, …)
# S3 method for psm
Survival(object, …) # for psm
# E.g. survival <- Survival(fit)
# S3 method for psm
Quantile(object, …) # for psm
# E.g. quantsurv <- Quantile(fit)
# S3 method for psm
Mean(object, …) # for psm
# E.g. meant <- Mean(fit)
# lambda(times, lp) # get hazard function at t=times, xbeta=lp
# survival(times, lp) # survival function at t=times, lp
# quantsurv(q, lp) # quantiles of survival time
# meant(lp) # mean survival time
# S3 method for psm
residuals(object, type=c("censored.normalized",
"response", "deviance", "dfbeta",
"dfbetas", "working", "ldcase", "ldresp", "ldshape", "matrix"), …)
# S3 method for residuals.psm.censored.normalized
survplot(fit, x, g=4, col, main, …)
# S3 method for residuals.psm.censored.normalized
lines(x, n=100, lty=1, xlim,
lwd=3, …)
# for type="censored.normalized"
Arguments
- formula
an S statistical model formula. Interactions up to third order are supported. The left hand side must be a
Surv
object.- object
a fit created by
psm
. Forsurvplot
with residuals frompsm
,object
is the result ofresiduals.psm
.- fit
a fit created by
psm
- data,subset,weights,dist,scale,init,na.action,control
see
survreg
.- parms
a list of fixed parameters. For the \(t\)-distribution this is the degrees of freedom; most of the distributions have no parameters.
- model
set to
TRUE
to include the model frame in the returned object- x
set to
TRUE
to include the design matrix in the object produced bypsm
. For thesurvplot
method,x
is an optional stratification variable (character, numeric, or categorical). Forlines.residuals.psm.censored.normalized
,x
is the result ofresiduals.psm
. Forprint
it is the result ofpsm
.- y
set to
TRUE
to include theSurv()
matrix- time.inc
setting for default time spacing. Used in constructing time axis in
survplot
, and also in make confidence bars. Default is 30 if time variable hasunits="Day"
, 1 otherwise, unless maximum follow-up time \(< 1\). Then max time/10 is used astime.inc
. Iftime.inc
is not given and max time/defaulttime.inc
is \(> 25\),time.inc
is increased.- correlation
set to
TRUE
to print the correlation matrix for parameter estimates- digits
number of places to print to the right of the decimal point
- coefs
specify
coefs=FALSE
to suppress printing the table of model coefficients, standard errors, etc. Specifycoefs=n
to print only the firstn
regression coefficients in the model.- title
a character string title to be passed to
prModFit
- …
other arguments to fitting routines, or to pass to
survplot
fromsurvplot.residuals.psm.censored.normalized
. Passed to the genericlines
function forlines
.- times
a scalar or vector of times for which to evaluate survival probability or hazard
- lp
a scalar or vector of linear predictor values at which to evaluate survival probability or hazard. If both
times
andlp
are vectors, they must be of the same length.- q
a scalar or vector of probabilities. The default is .5, so just the median survival time is returned. If
q
andlp
are both vectors, a matrix of quantiles is returned, with rows corresponding tolp
and columns toq
.- type
type of residual desired. Default is censored normalized residuals, defined as (link(Y) - linear.predictors)/scale parameter, where the link function was usually the log function. See
survreg
for other types.- n
number of points to evaluate theoretical standardized survival function for
lines.residuals.psm.censored.normalized
- lty
line type for
lines
, default is 1- xlim
range of times (or transformed times) for which to evaluate the standardized survival function. Default is range in normalized residuals.
- lwd
line width for theoretical distribution, default is 3
- g
number of quantile groups to use for stratifying continuous variables having more than 5 levels
- col
vector of colors for
survplot
method, corresponding to levels ofx
(must be a scalar if there is nox
)- main
main plot title for
survplot
. If omitted, is the name or label ofx
ifx
is given. Usemain=""
to suppress a title when you specifyx
.
Details
The object survreg.distributions
contains definitions of properties
of the various survival distributions.
psm
does not trap singularity errors due to the way survreg.fit
does matrix inversion. It will trap non-convergence (thus returning
fit$fail=TRUE
) if you give the argument failure=2
inside the
control
list which is passed to survreg.fit
. For example, use
f <- psm(S ~ x, control=list(failure=2, maxiter=20))
to allow up to
20 iterations and to set f$fail=TRUE
in case of non-convergence.
This is especially useful in simulation work.
Value
psm
returns a fit object with all the information survreg
would store as
well as what rms
stores and units
and time.inc
.
Hazard
, Survival
, and Quantile
return S-functions.
residuals.psm
with type="censored.normalized"
returns a
Surv
object which has a special attribute "theoretical"
which is used by the lines
routine. This is the assumed standardized survival function as a function
of time or transformed time.
See Also
rms
, survreg
,
residuals.survreg
, survreg.object
,
survreg.distributions
,
pphsm
, survplot
, survest
,
Surv
,
na.delete
,
na.detail.response
, datadist
,
latex.psm
, GiniMd
, prModFit
,
ggplot.Predict
, plot.Predict
Examples
# NOT RUN {
n <- 400
set.seed(1)
age <- rnorm(n, 50, 12)
sex <- factor(sample(c('Female','Male'),n,TRUE))
dd <- datadist(age,sex)
options(datadist='dd')
# Population hazard function:
h <- .02*exp(.06*(age-50)+.8*(sex=='Female'))
d.time <- -log(runif(n))/h
cens <- 15*runif(n)
death <- ifelse(d.time <= cens,1,0)
d.time <- pmin(d.time, cens)
f <- psm(Surv(d.time,death) ~ sex*pol(age,2),
dist='lognormal')
# Log-normal model is a bad fit for proportional hazards data
anova(f)
fastbw(f) # if deletes sex while keeping age*sex ignore the result
f <- update(f, x=TRUE,y=TRUE) # so can validate, compute certain resids
validate(f, B=10) # ordinarily use B=300 or more
plot(Predict(f, age, sex)) # needs datadist since no explicit age, hosp.
# Could have used ggplot(Predict(...))
survplot(f, age=c(20,60)) # needs datadist since hospital not set here
# latex(f)
S <- Survival(f)
plot(f$linear.predictors, S(6, f$linear.predictors),
xlab=expression(X*hat(beta)),
ylab=expression(S(6,X*hat(beta))))
# plots 6-month survival as a function of linear predictor (X*Beta hat)
times <- seq(0,24,by=.25)
plot(times, S(times,0), type='l') # plots survival curve at X*Beta hat=0
lam <- Hazard(f)
plot(times, lam(times,0), type='l') # similarly for hazard function
med <- Quantile(f) # new function defaults to computing median only
lp <- seq(-3, 5, by=.1)
plot(lp, med(lp=lp), ylab="Median Survival Time")
med(c(.25,.5), f$linear.predictors)
# prints matrix with 2 columns
# fit a model with no predictors
f <- psm(Surv(d.time,death) ~ 1, dist="weibull")
f
pphsm(f) # print proportional hazards form
g <- survest(f)
plot(g$time, g$surv, xlab='Time', type='l',
ylab=expression(S(t)))
f <- psm(Surv(d.time,death) ~ age,
dist="loglogistic", y=TRUE)
r <- resid(f, 'cens') # note abbreviation
survplot(npsurv(r ~ 1), conf='none')
# plot Kaplan-Meier estimate of
# survival function of standardized residuals
survplot(npsurv(r ~ cut2(age, g=2)), conf='none')
# both strata should be n(0,1)
lines(r) # add theoretical survival function
#More simply:
survplot(r, age, g=2)
options(datadist=NULL)
# }