Hi: I didn't see anything on first blush from the mod1 or summary(mod1) objects, but it's not too hard to compute:
> names(mod1) [1] "coefficients" "icoef" "var" [4] "var2" "loglik" "iter" [7] "linear.predictors" "frail" "fvar" [10] "df" "penalty" "pterms" [13] "assign2" "history" "printfun" [16] "scale" "idf" "df.residual" [19] "terms" "means" "call" [22] "dist" "y" # predicted frailties > mod1$frail [1] 3.5445680 0.1082518 -2.2070432 1.9199828 -0.9953027 0.6816612 [7] -1.8924137 1.0421686 -1.5761400 -0.6257329 # Since the degrees of freedom for the frailty term were 8.95 rather than 9, we need to make a small adjustment: with(mod1, var(frail) * (length(frail) - 1)/df[2]) [1] 3.366740 After checking str(mod1), it turns out that it can be extracted from mod1: > str(mod1) List of 23 $ coefficients : Named num 0.987 ..- attr(*, "names")= chr "(Intercept)" $ icoef : num [1:2] 1.678 0.645 $ var : num [1:2, 1:2] 0.33876 -0.00115 -0.00115 0.00655 ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:2] "(Intercept)" "Log(scale)" .. ..$ : chr [1:2] "(Intercept)" "Log(scale)" $ var2 : num [1:2, 1:2] 0.00209 -0.00115 -0.00115 0.00655 $ loglik : num [1:2] -287 -149 $ iter : num [1:2] 5 17 $ linear.predictors: num [1:100] 4.53 4.53 4.53 4.53 4.53 ... $ frail : num [1:10] 3.545 0.108 -2.207 1.92 -0.995 ... $ fvar : num [1:10] 0.353 0.354 0.354 0.353 0.354 ... $ df : num [1:3] 0.00617 8.94974 0.99994 $ penalty : num 19.7 $ pterms : Named num [1:2] 0 2 ..- attr(*, "names")= chr [1:2] "(Intercept)" "frailty.gaussian(Unit)" $ assign2 :List of 3 ..$ (Intercept) : int 1 ..$ frailty.gaussian(Unit): int 2 ..$ sigma : int 2 $ history :List of 1 ..$ frailty.gaussian(Unit):List of 3 .. ..$ theta : num 3.37 #### <<=================== HERE!! :) .. ..$ done : Named logi TRUE .. .. ..- attr(*, "names")= chr "resid" .. ..$ history: num [1:5, 1:4] 1 3 6 3.37 3.37 ... .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. ..$ : NULL .. .. .. ..$ : chr [1:4] "theta" "resid" "fsum" "trace" $ printfun :List of 1 ..$ frailty.gaussian(Unit):function (coef, var, var2, df, history) $ scale : num 0.434 $ idf : num 2 $ df.residual : num 90 $ terms :Classes 'terms', 'formula' length 3 Surv(Time, Status) ~ 1 + frailty.gaussian(Unit) .. ..- attr(*, "variables")= language list(Surv(Time, Status), frailty.gaussian(Unit)) .. ..- attr(*, "factors")= int [1:2, 1] 0 1 .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. ..$ : chr [1:2] "Surv(Time, Status)" "frailty.gaussian(Unit)" .. .. .. ..$ : chr "frailty.gaussian(Unit)" .. ..- attr(*, "term.labels")= chr "frailty.gaussian(Unit)" .. ..- attr(*, "specials")=Dotted pair list of 2 .. .. ..$ strata : NULL .. .. ..$ cluster: NULL .. ..- attr(*, "order")= int 1 .. ..- attr(*, "intercept")= int 1 .. ..- attr(*, "response")= int 1 .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> .. ..- attr(*, "predvars")= language list(Surv(Time, Status), frailty.gaussian(Unit)) .. ..- attr(*, "dataClasses")= Named chr [1:2] "nmatrix.2" "numeric" .. .. ..- attr(*, "names")= chr [1:2] "Surv(Time, Status)" "frailty.gaussian(Unit)" $ means : Named num [1:2] 1 5.5 ..- attr(*, "names")= chr [1:2] "(Intercept)" "frailty.gaussian(Unit)" $ call : language survreg(formula = Surv(Time, Status) ~ 1 + frailty.gaussian(Unit), data = test1) $ dist : chr "weibull" $ y : num [1:100, 1:2] 4.58 4.02 4.47 2.45 4.91 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:100] "1" "2" "3" "4" ... .. ..$ : NULL - attr(*, "class")= chr [1:2] "survreg.penal" "survreg" > mod1$history$frailty.gaussian$theta # same as the calculation (whew!) [1] 3.366740 Basically, you just need to dig a little... HTH, Dennis On Fri, Apr 8, 2011 at 3:52 AM, Dr. Pablo E. Verde < pabloemilio.ve...@uni-duesseldorf.de> wrote: > 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. > [[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.
