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.

Reply via email to