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 <cr...@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@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
ddalth...@usgs.gov
[[alternative HTML version deleted]]
______________________________________________
R-help@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.
