[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