[R] survreg and cluster and frailty
Andrew Beckerman
a.beckerman at sheffield.ac.uk
Fri Jan 4 14:30:59 CET 2008
Dear All - Any help/advice appreciated....
OSX 10.5.1, R 2.6.1, survival 2.34
I am trying to fit an interval censored model with cluster or frailty
term (first 20 lines of the data set are below).
-------MODEL--------
m1<-survreg(Surv(start,stop,event,type="interval")~tree+hag+entcir
+depth+d2n+d2f+div+cluster(nestide),ss)
-----DETAILS-------
The independent variables are habitat characteristics, and we are
modelling parrot survival, with nestide as nest from which more than
one parrot babies come from.... like a litter.
-----PROBLEM------
summary(m1) does not appear to produce a robust estimate for the
terms, for this model, or the examples in ?cluster.
------SOLUTION???---------
This link, indicating this potential problem, suggests a workaround....
http://tolstoy.newcastle.edu.au/R/help/06/03/22598.html
is it legitimate to apply this to the above model, for example, via:
2*pnorm(summary(m1)$table[,1]/summary(m1)$table[,2])
-------MOREOVER, attempting the same model with a frailty term
produces some interesting problems.... such as:
(1) the example for frailty does not seem to work at all, and defaults
to model without frailty
frailty.model <- survreg(Surv(time, status) ~ rx + frailty(litter),
rats )
Warning message:
In survpenal.fit(X, Y, weights, offset, init = init, controlvals =
control, :
Inner loop failed to coverge for iterations 2
> summary(frailty.model)
Call:
survreg(formula = Surv(time, status) ~ rx + frailty(litter),
data = rats)
Value Std. Error z p
(Intercept) 5.130 0.2198 23.34 1.76e-120
rx -0.208 0.0713 -2.92 3.48e-03
Log(scale) -1.666 0.1309 -12.73 4.11e-37
Scale= 0.189
Weibull distribution
Loglik(model)= -202.9 Loglik(intercept only)= -246.3
Chisq= 86.73 on 39.5 degrees of freedom, p= 2.1e-05
Number of Newton-Raphson Iterations: 10 71
n= 150
(2) Very simple parrot model - note that inner loop has failed....
what does this mean?
survreg(Surv(start,stop,event,type="interval")~div+frailty(nestide),ss)
Call:
survreg(formula = Surv(start, stop, event, type = "interval") ~
div + frailty(nestide), data = ss)
coef se(coef) se2 Chisq DF p
(Intercept) 3.66 0.76 0.358 23.21 1 1.5e-06
div -1.24 1.06 0.476 1.36 1 2.4e-01
frailty(nestide) 695.25 26 0.0e+00
Scale= 0.241
Iterations: 10 outer, 106 Newton-Raphson
Variance of random effect= 0.492 I-likelihood = -75
Degrees of freedom for terms= 0.2 0.2 26.0 1.0
Likelihood ratio test=175 on 25.4 df, p=0 n= 92
Warning message:
In survpenal.fit(X, Y, weights, offset, init = init, controlvals =
control, :
Inner loop failed to coverge for iterations 2 5
(3) Slightly more complex model - notice inner loop problem, and
complete loss of div as a coherent variable.... again, any ideas why?
> survreg(Surv(start,stop,event,type="interval")~div+hag+d2f+d2n+tree
+frailty(nestide),ss)
Call:
survreg(formula = Surv(start, stop, event, type = "interval") ~
div + hag + d2f + d2n + tree + frailty(nestide), data = ss)
coef se(coef) se2 Chisq DF p
(Intercept) 3.342919 7.88e-01 2.53e-01 17.99 1.0 2.2e-05
div 0.00e+00 0.00e+00 1.0
hag -0.000411 4.94e-04 1.50e-04 0.69 1.0 4.0e-01
d2f -0.000282 8.18e-05 2.49e-05 11.88 1.0 5.7e-04
d2n -0.000299 1.28e-04 3.45e-05 5.49 1.0 1.9e-02
tree 0.068125 4.03e-01 1.19e-01 0.03 1.0 8.7e-01
frailty(nestide) 249.60 23.9 0.0e+00
Scale= 0.239
Iterations: 10 outer, 82 Newton-Raphson
Variance of random effect= 0.649 I-likelihood = -63.3
Degrees of freedom for terms= 0.1 NaN 0.1 0.1 0.1 0.1 23.9 1.0
Likelihood ratio test=176 on NaN df, p=NaN n= 92
Warning message:
In survpenal.fit(X, Y, weights, offset, init = init, controlvals =
control, :
Inner loop failed to coverge for iterations 2
Many thanks,
Andrew
FIRST 20 Lines of data
> ss[1:20,]
start stop event tree hag entcir entarea depth nestid d2n d2f
1 1 8 1 1 170 45 192 78 1 850 7350
2 1 8 1 1 230 115 611 47 2 520 6810
3 1 8 1 1 150 64 156 70 3 590 6870
4 22 57 1 1 410 58 270 100 4 3000 1310
5 22 57 1 1 410 58 270 100 4 3000 1310
6 22 57 1 1 410 58 270 100 4 3000 1310
7 22 57 1 1 410 58 270 100 4 3000 1310
8 22 64 1 1 140 60 432 134 5 187 3090
9 29 64 1 1 140 60 432 134 5 187 3090
10 29 64 1 1 140 60 432 134 5 187 3090
11 29 64 1 1 140 60 432 134 5 187 3090
12 15 71 1 1 120 56 260 77 6 490 5920
13 15 64 1 1 120 56 260 77 6 490 5920
14 15 57 1 1 120 56 260 77 6 490 5920
15 15 100 0 1 170 43 168 104 7 409 3090
16 8 100 0 1 80 128 658 64 8 3360 3520
17 8 29 1 1 80 128 658 64 8 3360 3520
18 15 29 1 1 80 128 658 64 8 3360 3520
19 15 36 1 1 80 128 658 64 8 3360 3520
20 8 100 0 1 170 56 152 156 9 350 6040
div nestide
1 0.6027874 1
2 0.6121574 2
3 0.8094382 3
4 0.7353479 4
5 0.7353479 4
6 0.7353479 4
7 0.7353479 4
8 0.8645524 5
9 0.8645524 5
10 0.8645524 5
11 0.8645524 5
12 0.6458459 6
13 0.6458459 6
14 0.6458459 6
15 0.2722860 7
16 0.8397355 8
17 0.8397355 8
18 0.8397355 8
19 0.8397355 8
20 0.7203888 9
---------------------------------------------------------------------------------
Dr. Andrew Beckerman
Department of Animal and Plant Sciences, University of Sheffield,
Alfred Denny Building, Western Bank, Sheffield S10 2TN, UK
ph +44 (0)114 222 0026; fx +44 (0)114 222 0002
http://www.beckslab.staff.shef.ac.uk/
More information about the R-help
mailing list