[R] I'm getting confused with notation and terminology in output from weibull parametric survival model from survreg()

[Christopher W. Ryan]

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