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.Verde@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]]