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]]
______________________________________________
R-help@r-project.org mailing list
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.
