Christopher W. Ryan
2016-Jul-29 20:07 UTC
[R] I'm getting confused with notation and terminology in output from weibull parametric survival model from survreg()
I'm trying to run a Weibull parametric survival model for recurrent
event data, with subject-specific frailties, using survreg() in the
survival package, and I'm having trouble understanding the output and
its notation, and how that translates to some of the books I am using as
references (DF Moore, Applied Survival Analysis Using R; and Kleinbaum
and Klein, Survival Analysis A Self-Learning Text). I understand there
are different notations for different ways of parameterizing a Weibull
or a gamma distribution, and perhaps that's where I am getting hung up.
I also may be confusing "scale" for the Weibull distribution of the
survial times with "scale" for the gamma distribution of the
frailties.
My ultimate goal is to display example survival curves: say, one for
"typically frail" subjects, one for "extra-frail" subjects,
and one for
"not-so-frail" subjects. I'd like to get estimated
"frailty" for each of
my subjects; or at least the distribution of those frailties. Do I need
the parameters of the gamma distribution to do that? If so, how do I
extract them? Or are they readily available in the survreg object?
Here is what I have tried so far:
## create some data similar to my real data, in which
## vast majority of subjects had no event
id <- c(1:600, rep(601:630, each=1), rep(631:650, each=2), rep(651:656,
each=3), rep(677:679, each=4), rep(680, 5))
time <- c(rpois(lambda=800, 600), rpois(lambda=600, length(601:630)),
rpois(lambda=600, length(631:650)*2), rpois(lambda=600,
length(651:656)*3), rpois(lambda=600, length(677:679)*4),
rpois(lambda=600, 5))
event <- c(rep(0, 600), rep(1, (length(id) - 600)))
dd <- data.frame(id=id, time=time, event=event)
dd.2 <- dd[order(id, time), ]
str(dd.2)
table(table(dd.2$id))
# time until censoring, for those without events
summary(subset(dd.2, event==0, select=time))
library(survival)
Surv.1 <- Surv(time, event)
# model without frailties
model.1 <- survreg(Surv.1 ~ 1, data=dd.2, dist="weibull")
# add frailty term
model.2 <- survreg(Surv.1 ~ 1 + frailty(id), data=dd.2,
dist="weibull")
# should be same as above line
model.2.b <- survreg(Surv.1 ~ 1 + frailty(id,
distribution="gamma"),
data=dd.2, dist="weibull")
# I don't know if this is the right way to go about it
a.scale <- model.2$scale
var.X <- model.2$history$frailty$theta
s.shape <- sqrt(var.X/a.scale)
gamma.frail.x <- function(a,s,q){ 1/((s^a) * gamma(a)) * (q^(a-1) *
exp(-(q/s))) }
q <- seq(0.1, 10, by=0.2)
maybe.my.frailties <- gamma.frail.x(a.scale, s.shape, q)))
plot(density(maybe.my.frailties))
## end code
Or, would I be better off changing tactics and using frailtypack?
Thanks for any help. Session info is below, in case it is relevant.
--Chris Ryan
> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] survival_2.39-5
loaded via a namespace (and not attached):
[1] compiler_3.3.1 Matrix_1.2-6 tools_3.3.1 splines_3.3.1
[5] grid_3.3.1 lattice_0.20-33
Dalthorp, Daniel
2016-Jul-29 20:18 UTC
[R] I'm getting confused with notation and terminology in output from weibull parametric survival model from survreg()
The parameterization for Weibull in the 'survival' package corresponds to base R's dweibull, etc. suite as 1/scale --> shape and exp(coef[1]) --> scale On Fri, Jul 29, 2016 at 1:07 PM, Christopher W. Ryan <cryan at binghamton.edu> wrote:> I'm trying to run a Weibull parametric survival model for recurrent event > data, with subject-specific frailties, using survreg() in the survival > package, and I'm having trouble understanding the output and its notation, > and how that translates to some of the books I am using as references (DF > Moore, Applied Survival Analysis Using R; and Kleinbaum and Klein, Survival > Analysis A Self-Learning Text). I understand there are different notations > for different ways of parameterizing a Weibull or a gamma distribution, and > perhaps that's where I am getting hung up. I also may be confusing "scale" > for the Weibull distribution of the survial times with "scale" for the > gamma distribution of the frailties. > > My ultimate goal is to display example survival curves: say, one for > "typically frail" subjects, one for "extra-frail" subjects, and one for > "not-so-frail" subjects. I'd like to get estimated "frailty" for each of my > subjects; or at least the distribution of those frailties. Do I need the > parameters of the gamma distribution to do that? If so, how do I extract > them? Or are they readily available in the survreg object? > > Here is what I have tried so far: > > ## create some data similar to my real data, in which > ## vast majority of subjects had no event > id <- c(1:600, rep(601:630, each=1), rep(631:650, each=2), rep(651:656, > each=3), rep(677:679, each=4), rep(680, 5)) > time <- c(rpois(lambda=800, 600), rpois(lambda=600, length(601:630)), > rpois(lambda=600, length(631:650)*2), rpois(lambda=600, length(651:656)*3), > rpois(lambda=600, length(677:679)*4), rpois(lambda=600, 5)) > event <- c(rep(0, 600), rep(1, (length(id) - 600))) > dd <- data.frame(id=id, time=time, event=event) > dd.2 <- dd[order(id, time), ] > str(dd.2) > table(table(dd.2$id)) > # time until censoring, for those without events > summary(subset(dd.2, event==0, select=time)) > > library(survival) > Surv.1 <- Surv(time, event) > > # model without frailties > model.1 <- survreg(Surv.1 ~ 1, data=dd.2, dist="weibull") > > # add frailty term > model.2 <- survreg(Surv.1 ~ 1 + frailty(id), data=dd.2, dist="weibull") > > # should be same as above line > model.2.b <- survreg(Surv.1 ~ 1 + frailty(id, distribution="gamma"), > data=dd.2, dist="weibull") > > # I don't know if this is the right way to go about it > a.scale <- model.2$scale > var.X <- model.2$history$frailty$theta > s.shape <- sqrt(var.X/a.scale) > > gamma.frail.x <- function(a,s,q){ 1/((s^a) * gamma(a)) * (q^(a-1) * > exp(-(q/s))) } > q <- seq(0.1, 10, by=0.2) > > maybe.my.frailties <- gamma.frail.x(a.scale, s.shape, q))) > plot(density(maybe.my.frailties)) > > ## end code > > > Or, would I be better off changing tactics and using frailtypack? > > Thanks for any help. Session info is below, in case it is relevant. > > --Chris Ryan > > > > > sessionInfo() > R version 3.3.1 (2016-06-21) > Platform: x86_64-w64-mingw32/x64 (64-bit) > Running under: Windows 7 x64 (build 7601) Service Pack 1 > > locale: > [1] LC_COLLATE=English_United States.1252 > [2] LC_CTYPE=English_United States.1252 > [3] LC_MONETARY=English_United States.1252 > [4] LC_NUMERIC=C > [5] LC_TIME=English_United States.1252 > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] survival_2.39-5 > > loaded via a namespace (and not attached): > [1] compiler_3.3.1 Matrix_1.2-6 tools_3.3.1 splines_3.3.1 > [5] grid_3.3.1 lattice_0.20-33 > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >-- Dan Dalthorp, PhD USGS Forest and Rangeland Ecosystem Science Center Forest Sciences Lab, Rm 189 3200 SW Jefferson Way Corvallis, OR 97331 ph: 541-750-0953 ddalthorp at usgs.gov [[alternative HTML version deleted]]