Dear S-Plus users.
>From Dr Renaud Lancelot I got a quick reply which actually solved my
problem:
> >
> > Dear S-Plus users.
> >
> > Is there any S-Plus function for estimating parametric survival models
> > with frailty, for instance a Weibull-model with gamma-frailty like in
> > Aalen, Bjertness & Sønju : Statistics in Medicine 14, 1819-1829 (1995)?
Lancelot's reply:
>
> It is implemented in S+ 2000 r2 / Win (it probably comes from Pr
> Therneau's Surv5 library: please check). See on-line help file for
> frailty():
>
> DESCRIPTION
> Fit a penalized factor variable in a survival model; random effects or
> frailties are one instance of this.
>
> USAGE
> frailty(x, distribution="gamma", ...)
>
> REQUIRED ARGUMENTS
> x a variable describing discrete groups.
>
> OPTIONAL ARGUMENTS
> distribution character string specifying the distribution of the random
> effect. Currently, the legal values are "gamma", "gaussian", and "t".
> ... other arguments to the specific frailty function.
>
> VALUE
> a variable with attached attributes, so as to be recognized as a
> penalized term by coxph or survReg.
>
> DETAILS
> this routine calls frailty.gamma, frailty.gaussian, or frailty.t to do
> the actual work.
>
> SEE ALSO
> coxph, survReg, frailty.gamma, frailty.gaussian, frailty.t.
>
> EXAMPLES
> coxph(Surv(time, status) ~ age + sex + frailty(inst), data=lung,
> na.action=na.exclude)
>
While I was already aware of this function, my confusion arose from the
fact that I did not get an estimate of the frailty variance as the
standard output from applying this function:
*** Parametric Survival ***
Call:
survReg(formula = Surv(Obsdays, Restored == 3, type = "right") ~ Tunnty
+ Caract +
frailty(Patientn), data = oppeg5, na.action = na.exclude, dist =
"weibull",
scale = 0, control = list(maxiter = 30, rel.tolerance = 1e-005,
failure
= 1))
Value Std. Error z p
(Intercept) 9.442 0.696 13.564 6.54e-042
Tunnty -0.571 0.193 -2.956 3.12e-003
Caract -0.321 0.358 -0.897 3.70e-001
Log(scale) -1.435 0.147 -9.747 1.90e-022
Scale= 0.238
Weibull distribution
Loglik(model)= -174.4 Loglik(intercept only)= -226.9
Chisq= 105.05 on 45.5 degrees of freedom, p= 1.3e-006
Number of Newton-Raphson Iterations: 6 26
n= 302
Correlation of Coefficients:
(Intercept) Tunnty Caract
Tunnty -0.478
Caract -0.828 -0.007
Log(scale) 0.290 -0.145 -0.108
However, the following did the trick (gave me the frailty variance):
> weibfrail <- survReg(formula = Surv(Obsdays, Restored == 3, type = "right") ~
> Tunnty + Caract +
+ frailty(Patientn), data = oppeg5, na.action = na.exclude, dist =
"weibull",
+ scale = 0, control = list(maxiter = 30, rel.tolerance = 1e-005,
failure = 1))
> weibfrail
Call:
survReg(formula = Surv(Obsdays, Restored == 3, type = "right") ~ Tunnty
+ Caract +
frailty(Patientn), data = oppeg5, na.action = na.exclude, dist =
"weibull",
scale = 0, control = list(maxiter = 30, rel.tolerance = 1e-005,
failure = 1))
coef se(coef) se2 Chisq DF p
(Intercept) 9.442 0.696 0.587 183.98 1.0 0.0e+000
Tunnty -0.571 0.193 0.185 8.74 1.0 3.1e-003
Caract -0.321 0.358 0.279 0.80 1.0 3.7e-001
frailty(Patientn) 144.91 44.4 1.3e-012
Scale= 0.238
Iterations: 6 outer, 26 Newton-Raphson
Variance of random effect= 1.18 EM likelihood = -43
Degrees of freedom for terms= 0.7 0.9 0.6 44.4 1.0
Likelihood ratio test=105 on 45.5 df, p=1.34e-006 n= 302
>
I suggest the frailty parameter should be included in the standard
summary output!
Best regards
Geir E. Eide
>
>
> --
> Dr Renaud Lancelot, vétérinaire
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
|