I have the following questions about the variance of the random effects in the survreg() function in the survival package:

1) How can I extract the variance of the random effects after fitting a model?

For example:

set.seed(1007)
x <- runif(100)
m <- rnorm(10, mean = 1, sd =2)
mu <- rep(m, rep(10,10))

test1 <- data.frame(Time = qsurvreg(x, mean = mu, scale= 0.5, distribution = "weibull"),
                    Status = rep(1, 100),
                    Unit = gl(10,10)
                    )

mod1 <- survreg(Surv(Time, Status) ~ 1 + frailty.gaussian(Unit), data = test1)
> mod1
...
                       coef  se(coef) se2    Chisq DF   p      
(Intercept)            0.987 0.582    0.0457  2.87 1.00 9.0e-02
frailty.gaussian(Unit)                       85.26 8.95 1.4e-14

Scale= 0.434 

Iterations: 5 outer, 17 Newton-Raphson
     Variance of random effect= 3.37 
...

It is not clear from the returned list how to get the printed 3.37.

2) What is the meaning of the variance of the random effects if we fit gamma frailities?

For example:

set.seed(1007)
x <- runif(100)
# gamma frailties
m <- rgamma(10, shape = 1, scale = 2) # E(m) = 2, var(m) = 4
mu <- rep(log(m), rep(10,10))

test2 <- data.frame(Time = qsurvreg(x, mean = mu, scale= 0.5, distribution = "weibull"),
                    Status = rep(1, 100),
                    Unit = gl(10,10)
                    )

mod2 <- survreg(Surv(Time, Status) ~ 1 + frailty.gamma(Unit), data = test2)
mod2

...
                    coef se(coef) se2    Chisq  DF   p    
(Intercept)         1.06 0.45     0.0581   5.59 1.00 0.018
frailty.gamma(Unit)                      108.67 8.92 0.000

Scale= 0.434 

Iterations: 10 outer, 33 Newton-Raphson
     Variance of random effect= 1.99   I-likelihood = -95.7 
...

In this case I expected  that the variance of random effects is close to 4.

Thanks! 

Pablo






	[[alternative HTML version deleted]]

