1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
|
S-PLUS : Copyright (c) 1988, 1996 MathSoft, Inc.
S : Copyright AT&T.
Version 3.4 Release 1 for Sun SPARC, SunOS 5.3 : 1996
Mayo local startup...
loading survival
loading date
loading slocal.misc
set missing value action to 'na.omit'
set contrasts to ('contr.treatment', 'contr.poly')
stop automatic character -> factor conversion in data frames
Working data will be in .Data
> attach("../.Data")
> dyn.load("../loadmod.o")
> postscript(file="testall.ps")
> options(na.action="na.omit", contrasts="contr.treatment")
> #
> # This data set caused problems for Splus 3.4 due to a mistake in
> # my initial value code. Data courtesy Bob Treder at Statsci
> #
> capacitor <- read.table("data.capacitor", row.names=1,
+ col.names=c("", "days", "event", "voltage"))
> > fitig <- survreg(Surv(days, event)~voltage,
+ dist = "gaussian", data = capacitor)
> summary(fitig)
Call:
survreg(formula = Surv(days, event) ~ voltage, data = capacitor, dist =
"gaussian")
Value Std. Error z p
(Intercept) 1764.9 163.387 10.80 3.36e-27
voltage -53.9 5.545 -9.72 2.56e-22
Log(scale) 4.8 0.105 45.56 0.00e+00
Scale= 121
Gaussian distribution
Loglik(model)= -361.9 Loglik(intercept only)= -420.1
Chisq= 116.33 on 1 degrees of freedom, p= 0
Number of Newton-Raphson Iterations: 5
n= 125
Correlation of Coefficients:
(Intercept) voltage
voltage -0.996
Log(scale) 0.412 -0.384
> > fitix <- survreg(Surv(days, event)~voltage,
+ dist = "extreme", data = capacitor)
> summary(fitix)
Call:
survreg(formula = Surv(days, event) ~ voltage, data = capacitor, dist =
"extreme")
Value Std. Error z p
(Intercept) 2055.59 180.348 11.4 4.28e-30
voltage -62.21 5.967 -10.4 1.88e-25
Log(scale) 4.53 0.108 41.9 0.00e+00
Scale= 92.9
Extreme value distribution
Loglik(model)= -360 Loglik(intercept only)= -427.1
Chisq= 134.25 on 1 degrees of freedom, p= 0
Number of Newton-Raphson Iterations: 6
n= 125
Correlation of Coefficients:
(Intercept) voltage
voltage -0.998
Log(scale) 0.425 -0.420
> > fitil <- survreg(Surv(days, event)~voltage,
+ dist = "logistic", data = capacitor)
> summary(fitil)
Call:
survreg(formula = Surv(days, event) ~ voltage, data = capacitor, dist =
"logistic")
Value Std. Error z p
(Intercept) 1811.56 148.853 12.2 4.48e-34
voltage -55.48 4.986 -11.1 9.39e-29
Log(scale) 4.19 0.117 35.8 2.03e-280
Scale= 66.3
Logistic distribution
Loglik(model)= -360.4 Loglik(intercept only)= -423.7
Chisq= 126.5 on 1 degrees of freedom, p= 0
Number of Newton-Raphson Iterations: 5
n= 125
Correlation of Coefficients:
(Intercept) voltage
voltage -0.996
Log(scale) 0.343 -0.321
> > rm(fitil, fitig, fitix)
> #
> # Good initial values are key to this data set
> # It killed v4 of survreg;
> # data courtesy of Deborah Donnell, Fred Hutchinson Cancer Center
> #
>
> donnell <- scan("data.donnell", what=list(time1=0, time2=0, status=0))
> donnell <- data.frame(donnell)
> > dfit <- survreg(Surv(time1, time2, status, type="interval") ~1, donnell)
> summary(dfit)
Call:
survreg(formula = Surv(time1, time2, status, type = "interval") ~ 1, data =
donnell)
Value Std. Error z p
(Intercept) 2.390 0.804 2.973 0.00295
Log(scale) -0.237 0.346 -0.687 0.49222
Scale= 0.789
Weibull distribution
Loglik(model)= -51 Loglik(intercept only)= -51
Number of Newton-Raphson Iterations: 9
n= 210
Correlation of Coefficients:
(Intercept)
Log(scale) 0.955
> > #
> # Do a contour plot of the donnell data
> #
> npt <- 25
> beta0 <- seq(.4, 2.4, length=npt)
> logsig <- seq(-1.4, 0.41, length=npt)
> donlog <- matrix(0,npt, npt)
> > for (i in 1:npt) {
+ for (j in 1:npt) {
+ fit <- survreg(Surv(time1, time2, status, type="interval") ~1,
+ donnell, init=c(beta0[i],logsig[j]),
+ control=list(maxiter=0));+ donlog[i,j] <- fit$log[1];+ };+ }
> > clev <- -c(51, 51.5, 52:60, 65, 75, 85, 100, 150)
> contour(beta0, logsig, pmax(donlog, -200), levels=clev, xlab="Intercept",
+ ylab="Log(sigma)")
> points(2.39, log(.7885), pch=1, col=2)
> title("Donnell data")
> > #
> # Compute the path of the iteration
> # Step 2 isn't so good, and is followed by 3 iters of step-halving
> #
> niter <- 14
> donpath <- matrix(0,niter+1,2)
> for (i in 0:niter){
+ fit <- survreg(Surv(time1, time2, status, type="interval") ~1,
+ donnell, maxiter=i);+ donpath[i+1,] <- c(fit$coef, log(fit$scale));+ }
Warning messages:
1: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights,
offset, init = init, controlvals = controlvals, dist ....
2: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights,
offset, init = init, controlvals = controlvals, dist ....
3: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights,
offset, init = init, controlvals = controlvals, dist ....
4: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights,
offset, init = init, controlvals = controlvals, dist ....
5: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights,
offset, init = init, controlvals = controlvals, dist ....
6: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights,
offset, init = init, controlvals = controlvals, dist ....
7: Ran out of iterations and did not converge in: survreg.fit(X, Y, weights,
offset, init = init, controlvals = controlvals, dist ....
> points(donpath[,1], donpath[,2])
Points out of bounds X= 3.890545 Y= 0.684751
Points out of bounds X= 2.814988 Y= -0.01247539
Points out of bounds X= 2.51267 Y= -0.182087
> lines(donpath[,1], donpath[,2], col=4)
Lines out of bounds X= 3.890545 Y= 0.684751
Lines out of bounds X= 2.814988 Y= -0.01247539
Lines out of bounds X= 2.51267 Y= -0.182087
Lines out of bounds X= 2.385003 Y= -0.2381598
> > rm(beta0, logsig, niter, fit, npt, donlog, clev)
> #lfit1 <- censorReg(censor(time, status) ~ age + ph.ecog + strata(sex),lung)
> lfit2 <- survreg(Surv(time, status) ~ age + ph.ecog + strata(sex), lung)
> lfit3 <- survreg(Surv(time, status) ~ sex + (age+ph.ecog)*strata(sex), lung)
> > lfit4 <- survreg(Surv(time, status) ~ age + ph.ecog , lung,
+ subset=(sex==1))
> lfit5 <- survreg(Surv(time, status) ~ age + ph.ecog , lung,
+ subset=(sex==2))
> > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
> #aeq(lfit4$coef, lfit1[[1]]$coef)
> #aeq(lfit4$scale, lfit1[[1]]$scale)
> aeq(c(lfit4$scale, lfit5$scale), lfit3$scale )
[1] "Mean relative difference: 1.364018e-07"
> aeq(c(lfit4$scale, lfit5$scale), sapply(lfit1, function(x) x$scale))
Error: Object "lfit1" not found
Dumped
> > #
> # Test out ridge regression and splines
> #
> lfit0 <- survreg(Surv(time, status) ~1, lung)
> lfit1 <- survreg(Surv(time, status) ~ age + ridge(ph.ecog, theta=5), lung)
> lfit2 <- survreg(Surv(time, status) ~ sex + ridge(age, ph.ecog, theta=1), lung)
> lfit3 <- survreg(Surv(time, status) ~ sex + age + ph.ecog, lung)
> > lfit0
Call:
survreg(formula = Surv(time, status) ~ 1, data = lung)
Coefficients:
(Intercept)
6.034903
Scale= 0.7593932
Loglik(model)= -1153.9 Loglik(intercept only)= -1153.9
n= 228
> lfit1
Call:
survreg(formula = Surv(time, status) ~ age + ridge(ph.ecog, theta = 5), data =
lung)
coef se(coef) se2 Chisq DF p
(Intercept) 6.83082 0.42860 0.42860 254.0 1 0.00000
age -0.00783 0.00687 0.00687 1.3 1 0.25000
ridge(ph.ecog) -0.32032 0.08484 0.08405 14.2 1 0.00016
Scale= 0.738
Iterations: 1 outer, 4 Newton-Raphson
Degrees of freedom for terms= 1 1 1 1
Likelihood ratio test=18.6 on 2 df, p=8.73e-05
n=227 (1 observations deleted due to missing)
> lfit2
Call:
survreg(formula = Surv(time, status) ~ sex + ridge(age, ph.ecog, theta = 1),
data = lung)
coef se(coef) se2 Chisq DF p
(Intercept) 6.27163 0.45280 0.45210 191.84 1 0.0e+00
sex 0.40096 0.12371 0.12371 10.50 1 1.2e-03
ridge(age) -0.00746 0.00675 0.00674 1.22 1 2.7e-01
ridge(ph.ecog) -0.33848 0.08329 0.08314 16.51 1 4.8e-05
Scale= 0.731
Iterations: 1 outer, 5 Newton-Raphson
Degrees of freedom for terms= 1 1 2 1
Likelihood ratio test=30 on 3 df, p=1.37e-06
n=227 (1 observations deleted due to missing)
> lfit3
Call:
survreg(formula = Surv(time, status) ~ sex + age + ph.ecog, data = lung)
Coefficients:
(Intercept) sex age ph.ecog
6.273435 0.4010877 -0.007475331 -0.3396365
Scale= 0.7311049
Loglik(model)= -1132.4 Loglik(intercept only)= -1147.4
Chisq= 29.98 on 3 degrees of freedom, p= 1.4e-06
n=227 (1 observations deleted due to missing)
> > > xx <- pspline(lung$age, nterm=3, theta=.3)
> xx <- matrix(unclass(xx), ncol=ncol(xx)) # the raw matrix
> lfit4 <- survreg(Surv(time, status) ~xx, lung)
> lfit5 <- survreg(Surv(time, status) ~age, lung)
> > lfit6 <- survreg(Surv(time, status)~pspline(age, df=2), lung)
> plot(lung$age, predict(lfit6), xlab="Age", ylab="Spline prediction")
> title("Lung Data")
> > lfit7 <- survreg(Surv(time, status) ~ offset(lfit6$lin), lung)
> > lfit4
Call:
survreg(formula = Surv(time, status) ~ xx, data = lung)
Coefficients:
(Intercept) xx1 xx2 xx3 xx4 xx5
13.5507 -7.615118 -7.423983 -7.532781 -7.570687 -14.52685
Scale= 0.7557376
Loglik(model)= -1150.1 Loglik(intercept only)= -1153.9
Chisq= 7.52 on 5 degrees of freedom, p= 0.19
n= 228
> lfit5
Call:
survreg(formula = Surv(time, status) ~ age, data = lung)
Coefficients:
(Intercept) age
6.887117 -0.01360819
Scale= 0.7587492
Loglik(model)= -1151.9 Loglik(intercept only)= -1153.9
Chisq= 3.91 on 1 degrees of freedom, p= 0.048
n= 228
> lfit6
Call:
survreg(formula = Surv(time, status) ~ pspline(age, df = 2), data = lung)
coef se(coef) se2 Chisq DF p
(Intercept) 6.5918 0.63681 0.41853 107.15 1.00 0.000
pspline(age, df = 2), lin -0.0136 0.00687 0.00687 3.94 1.00 0.047
pspline(age, df = 2), non 0.78 1.06 0.400
Scale= 0.756
Iterations: 4 outer, 9 Newton-Raphson
Theta= 0.926
Degrees of freedom for terms= 0.4 2.1 1.0
Likelihood ratio test=5.2 on 1.5 df, p=0.0441 n= 228
> lfit7$coef
(Intercept)
-6.008005e-07
> > rm(lfit1, lfit2, lfit3, lfit4, lfit5, lfit6, lfit7)
> rm(xx, lfit0)
> #
> # Data courtesy of Bercedis Peterson, Duke University.
> # v4 of survreg fails due to 2 groups that have only 1 subject; the coef
> # for them easily gets out of hand. In fact, this data set is my toughest
> # test of the minimizer.
> #
> # A shrinkage model for this coefficient is therefore interesting
>
>
> peterson <- data.frame(
+ scan("data.peterson", what=list(grp=0, time=0, status=0)))
> > fitp <- survreg(Surv(time, status) ~ factor(grp), peterson)
> summary(fitp)
Call:
survreg(formula = Surv(time, status) ~ factor(grp), data = peterson)
Value Std. Error z p
(Intercept) 2.291 0.115 19.92 2.93e-88
factor(grp)2 0.786 0.177 4.44 8.79e-06
factor(grp)3 0.728 0.183 3.97 7.09e-05
factor(grp)4 -1.598 0.218 -7.32 2.48e-13
factor(grp)5 -0.500 0.218 -2.29 2.21e-02
factor(grp)6 0.475 0.170 2.79 5.23e-03
Log(scale) -1.684 0.257 -6.54 6.09e-11
Scale= 0.186
Weibull distribution
Loglik(model)= -26.7 Loglik(intercept only)= -40.7
Chisq= 28.18 on 5 degrees of freedom, p= 3.4e-05
Number of Newton-Raphson Iterations: 8
n= 19
Correlation of Coefficients:
(Intercept) factor(grp)2 factor(grp)3 factor(grp)4 factor(grp)5
factor(grp)2 -0.668
factor(grp)3 -0.683 0.463
factor(grp)4 -0.527 0.352 0.360
factor(grp)5 -0.527 0.352 0.360 0.278
factor(grp)6 -0.617 0.406 0.400 0.325 0.325
Log(scale) -0.364 0.285 0.380 0.192 0.192
factor(grp)6
factor(grp)2
factor(grp)3
factor(grp)4
factor(grp)5
factor(grp)6
Log(scale) 0.083
> > # Now a shrinkage model. Give the group coefficients
> # about 1/2 the scale parameter of the original model, i.e., .18.
> #
> ffit <- survreg(Surv(time, status) ~ frailty(grp, theta=.1), peterson)
> ffit
Call:
survreg(formula = Surv(time, status) ~ frailty(grp, theta = 0.1), data =
peterson)
coef se(coef) se2 Chisq DF p
(Intercept) 2.62 0.172 0.0874 232.0 1.00 0.0000
frailty(grp, theta = 0.1) 10.4 2.15 0.0067
Scale= 0.301
Iterations: 1 outer, 6 Newton-Raphson
Variance of random effect= 0.1 EM likelihood = -11.8
Degrees of freedom for terms= 0.3 2.2 0.7
Likelihood ratio test=13.8 on 1.1 df, p=0.00027 n= 19
> > #
> # Try 3 degrees of freedom Gaussian fit, since there are 6 groups.
> # Compare them to the unconstrained ones. The frailty coefs are
> # on a "sum to 0" constraint rather than "first coef=0", so
> # some conversion is neccessary
> #
> ffit3 <- survreg(Surv(time, status) ~ frailty(grp, df=3, dist="gauss"),
+ peterson)
> print(ffit3)
Call:
survreg(formula = Surv(time, status) ~ frailty(grp, df = 3, dist = "gauss"),
data = peterson)
coef se(coef) se2 Chisq DF p
(Intercept) 2.44 0.223 0.0661 119.8 1 0.00000
frailty(grp, df = 3, dist 16.4 3 0.00096
Scale= 0.251
Iterations: 7 outer, 24 Newton-Raphson
Variance of random effect= 0.197
Degrees of freedom for terms= 0.1 3.0 0.6
Likelihood ratio test=20.1 on 1.7 df, p=2.79e-05 n= 19
> > temp <- mean(c(0, fitp$coef[-1]))
> temp2 <- c(fitp$coef[1] + temp, c(0,fitp$coef[-1]) - temp)
> xx <- rbind(c(nrow(peterson), table(peterson$grp)),
+ temp2,
+ c(ffit3$coef, ffit3$frail))
> dimnames(xx) <- list(c("N", "factor fit", "frailty fit"),
+ c("Intercept", paste("grp", 1:6)))
> signif(xx,2)
Intercept grp 1 grp 2 grp 3 grp 4 grp 5 grp 6
N 19.0 3.000 6.00 6.00 1.00 1.00 2.00
factor fit 2.3 0.018 0.80 0.75 -1.60 -0.48 0.49
frailty fit 2.4 -0.180 0.58 0.55 -0.77 -0.44 0.26
> #
> # All but the first coef are shrunk towards zero.
> #
> rm(ffit, ffit3, temp, temp2, xx, fitp)
> > #
> # Look at predicted values
> #
> ofit1 <- survreg(Surv(futime, fustat) ~ age + ridge(ecog.ps, rx), ovarian)
> > predict(ofit1)
[1] 207.7546 172.7985 358.7725 1426.6414 1353.7225 843.8571 1102.1610
[8] 859.5061 416.3272 1280.4037 820.7276 1882.7133 876.1244 1041.8917
[15] 3477.0123 2622.9581 3761.4852 2207.8635 1362.1943 3113.9504 879.1986
[22] 180.8418 2501.0478 645.2425 555.8297 936.0066
> predict(ofit1, type="response")
[1] 207.7546 172.7985 358.7725 1426.6414 1353.7225 843.8571 1102.1610
[8] 859.5061 416.3272 1280.4037 820.7276 1882.7133 876.1244 1041.8917
[15] 3477.0123 2622.9581 3761.4852 2207.8635 1362.1943 3113.9504 879.1986
[22] 180.8418 2501.0478 645.2425 555.8297 936.0066
> predict(ofit1, type="terms", se=T)
$fit:
age ridge(ecog.ps, rx)
1 -1.37775562 -0.1765498
2 -1.56198696 -0.1765498
3 -0.87785012 -0.1301245
4 0.23871941 0.1336957
5 0.49650010 -0.1765498
6 -0.02255551 -0.1301245
7 -0.06575616 0.1801210
8 -0.31442628 0.1801210
9 -0.68264179 -0.1765498
10 0.08414645 0.1801210
11 -0.05034745 -0.1301245
12 0.51611042 0.1336957
13 -0.29527617 0.1801210
14 -0.07556559 0.1336957
15 1.43981531 -0.1765498
16 1.11151929 -0.1301245
17 1.47203044 -0.1301245
18 0.98566722 -0.1765498
19 0.19249326 0.1336957
20 1.01928857 0.1336957
21 -0.29177342 0.1801210
22 -1.56291591 -0.1301245
23 1.11035170 -0.1765498
24 -0.60115796 0.1801210
25 -0.70389695 0.1336957
26 -0.18273628 0.1336957
attr($fit, "constant"):
[1] 6.890663
$se.fit:
age ridge(ecog.ps, rx)
1 0.356016316 0.1738687
2 0.403622262 0.1738687
3 0.226839188 0.1942911
4 0.061685835 0.1803885
5 0.128297162 0.1738687
6 0.005828413 0.1942911
7 0.016991596 0.1872282
8 0.081248722 0.1872282
9 0.176396749 0.1738687
10 0.021743704 0.1872282
11 0.013009936 0.1942911
12 0.133364529 0.1803885
13 0.076300276 0.1872282
14 0.019526381 0.1803885
15 0.372052731 0.1738687
16 0.287220023 0.1942911
17 0.380377220 0.1942911
18 0.254699460 0.1738687
19 0.049740854 0.1803885
20 0.263387319 0.1803885
21 0.075395153 0.1872282
22 0.403862307 0.1942911
23 0.286918315 0.1738687
24 0.155341078 0.1872282
25 0.181889150 0.1803885
26 0.047219621 0.1803885
> > temp1 <- predict(ofit1, se=T)
> temp2 <- predict(ofit1, type= "response", se=T)
> all.equal(temp2$se.fit, temp1$se.fit* sqrt(exp(temp1$fit)))
[1] "Mean relative difference: Inf"
> #
> # The Stanford data from 1980 is used in Escobar and Meeker
> # t5 = T5 mismatch score
> # Their case numbers correspond to a data set sorted by age
> #
> stanford2 <- read.table("data.stanford",
+ col.names=c("id", "time", "status", "age", "t5"))
> > stanford2$t5 <- ifelse(stanford2$t5 <0, NA, stanford2$t5)
> stanford2 <- stanford2[order(stanford2$age, stanford2$time),]
> stanford2$time <- ifelse(stanford2$time==0, .5, stanford2$time)
> > cage <- stanford2$age - mean(stanford2$age)
> fit1 <- survreg(Surv(time, status) ~ cage + cage^2, stanford2,
+ dist="lognormal")
> fit1
Call:
survreg(formula = Surv(time, status) ~ cage + cage^2, data = stanford2, dist =
"lognormal")
Coefficients:
(Intercept) cage I(cage^2)
6.717596 -0.06190903 -0.003504326
Scale= 2.362866
Loglik(model)= -863.6 Loglik(intercept only)= -868.8
Chisq= 10.5 on 2 degrees of freedom, p= 0.0053
n= 184
> ldcase <- resid(fit1, type="ldcase")
> ldresp <- resid(fit1, type="ldresp")
> print(ldresp)
139 159 181 119 74 120 99
0.1379203 0.145245 0.02628074 0.07320179 0.07624326 0.0399479 0.06328466
108 179 43 134 160 177 153
0.0612898 0.009685606 0.04767553 0.02980549 0.1036051 0.008990546 0.02114946
136 133 176 66 157 114 46
0.0255769 0.1591464 0.008618358 0.03389346 0.01141316 0.01990885 0.02044978
65 184 88 182 180 163
0.02480539 1.085676e-05 0.0547439 0.001786473 0.002574794 0.007654047
84 90 68 48 174 151 125
0.02024456 0.08561197 0.03894985 0.07007566 0.0037674 0.008314653 0.01248552
73 105 117 96 39 38 106
0.01954895 0.01831982 0.01739301 0.01789441 0.02406183 0.02364314 0.04717185
14 123 135 111 83 143 69
0.02051897 0.04763901 0.01663805 0.01367015 0.03204509 0.01857902 0.02058868
27 113 167 156 141 30
0.03896725 0.03775025 0.005091493 0.01528402 0.008682116 0.01746136
144 158 79 102 77 36
0.02593291 0.006620361 0.01375918 0.01547852 0.01786267 0.0233067
183 122 162 121 87 2
3.720795e-05 0.01696469 0.005954799 0.01233286 0.01655939 0.1089489
64 150 85 71 19 21 175
0.06015393 0.007469416 0.016665 0.01893414 0.02645489 0.18433 0.01789942
169 148 138 98 104 103 12
0.004379942 0.007619682 0.009332594 0.014288 0.01445961 0.01449499 0.03404299
89 3 100 55 142 63 168
0.03358406 0.03113308 0.01412657 0.01179741 0.00864158 0.01426955 0.00455403
72 137 10 124 17 94 82
0.01094162 0.009645953 0.01226565 0.01222511 0.01088511 0.01493685 0.0184422
170 149 42 128 67 109 75
0.03988065 0.03038322 0.02127744 0.01439502 0.01285836 0.00894498 0.0199779
26 97 58 178 140 32 126
0.02757124 0.02549339 0.02356049 0.002057497 0.01269584 0.01103393 0.0125303
51 101 29 33 164 60
0.0143023 0.01637414 0.02201027 0.01118993 0.006417554 0.008492277
152 145 112 76 47 118
0.008651511 0.009608663 0.01609214 0.02168279 0.02622512 0.02274276
5 129 31 35 40 130
0.01184996 0.009391146 0.008772106 0.008526052 0.009451662 0.01295997
28 56 91 44 23 37 70
0.01285987 0.01536642 0.02031498 0.02807957 0.01965943 0.01733256 0.009129009
132 9 81 59 127 131
0.009121716 0.009083024 0.01025238 0.01032187 0.01183694 0.01403298
80 20 25 165 24 172 146
0.02363944 0.02181251 0.02723395 0.02043511 0.02019542 0.01152649 0.01265906
86 107 95 116 41 61 155
0.01538527 0.02107502 0.0229847 0.02128395 0.01791007 0.01763098 0.01345059
166 154 4 92 93 62 34
0.01285115 0.0121809 0.01470506 0.02599207 0.03098464 0.03037746 0.02166522
15 173 171 52 110 50
0.01478523 0.007517959 0.008681577 0.01679632 0.02540017 0.03470671
45 53 54 147 115 16 1
0.03229507 0.03017737 0.02416304 0.01870027 0.02172489 0.1164272 0.04257799
6 7 57 78 161 11 8
0.02459122 0.03585529 0.03587691 0.02865161 0.02603297 0.05640971 0.04338251
49 13 22 18
0.03425475 0.06262793 0.1029315 0.1442429
> # The ldcase and ldresp should be compared to table 1 in Escobar and
> # Meeker, Biometrics 1992, p519; the colum they label as (1/2) A_{ii}
>
> plot1 <- function() {
+ # make their figure 1, 2, and 6
+ plot(stanford2$age, stanford2$time, log="y", xlab="Age", ylab="Days",
+ ylim=c(.01, 10^6));+ temp <- predict(fit1, type="response", se.fit=T) ;+ matlines(stanford2$age, cbind(temp$fit, temp$fit-1.96*temp$se.fit,
+ temp$fit+1.96*temp$se.fit),
+ lty=c(1,2,2));+ # these are the wrong CI lines, he plotted std dev, I plotted std err
+ # here are the right ones
+ # Using uncentered age gives different coefs, but makes prediction over an
+ # extended range somewhat simpler
+ refit <- survreg(Surv(time,status)~ age + age^2, stanford2,
+ dist="lognormal");+ plot(stanford2$age, stanford2$time, log="y", xlab="Age", ylab="Days",
+ ylim=c(.01, 10^6), xlim=c(0,75));+ temp2 <- predict(refit, list(age=1:75), type="quantile", p=c(.05, .5, .95));+ matlines(1:75, temp2, lty=c(1,2,2), col=2);+
+ tsplot(ldcase, xlab="Case Number", ylab="(1/2) A");+ title (main="Case weight pertubations");+ tsplot(ldresp, xlab="Case Number", ylab="(1/2) A");+ title(main="Response pertubations");+ }
> > plot1()
Warning: Data values <=0 omitted from logarithmic plot
> #
> # Stanford predictions in other ways
> #
> fit2 <- survreg(Surv(time, status) ~ poly(age,2), stanford2,
+ dist="lognormal")
> > p1 <- predict(fit1, type="response")
> p2 <- predict(fit2, type="response")
> aeq(p1, p2)
[1] T
> > p3 <- predict(fit2, type="terms", se=T)
> p4 <- predict(fit2, type="lp", se=T)
> p5 <- predict(fit1, type="lp", se=T)
> aeq(p3$fit + attr(p3$fit, "constant"), p4$fit)
[1] T
> aeq(p4$fit, p5$fit)
[1] T
> aeq(p3$se.fit, p4$se.fit) #this one should be false
[1] "Mean relative difference: 0.358807"
> aeq(p4$se.fit, p5$se.fit) #this one true
[1] T
> > #
> # Verify that scale can be fixed at a value
> # coefs will differ slightly due to different iteration paths
> tol <- survreg.control()$rel.tolerance
> > # Intercept only models
> fit1 <- survreg(Surv(time,status) ~ 1, lung)
> fit2 <- survreg(Surv(time,status) ~ 1, lung, scale=fit1$scale)
> all.equal(fit1$coef, fit2$coef, tolerance= tol)
[1] T
> all.equal(fit1$loglik, fit2$loglik, tolerance= tol)
[1] T
> > # multiple covariates
> fit1 <- survreg(Surv(time,status) ~ age + ph.karno, lung)
> fit2 <- survreg(Surv(time,status) ~ age + ph.karno, lung,
+ scale=fit1$scale)
> all.equal(fit1$coef, fit2$coef, tolerance=tol)
[1] T
> all.equal(fit1$loglik[2], fit2$loglik[2], tolerance=tol)
[1] T
> > # penalized models
> fit1 <- survreg(Surv(time, status) ~ pspline(age), lung)
> fit2 <- survreg(Surv(time, status) ~ pspline(age), lung, scale=fit1$scale)
> all.equal(fit1$coef, fit2$coef, tolerance=tol)
[1] "Mean relative difference: 0.0002487206"
> all.equal(fit1$loglik[2], fit2$loglik[2], tolerance=tol)
[1] T
> > rm(fit1, fit2, tol)
> > #
> # Test out the strata capabilities
> #
> tol <- survreg.control()$rel.tolerance
> aeq <- function(x,y,...) all.equal(as.vector(x), as.vector(y), ...)
> > # intercept only models
> fit1 <- survreg(Surv(time, status) ~ strata(sex), lung)
> fit2 <- survreg(Surv(time, status) ~ strata(sex) + sex, lung)
> fit3a<- survreg(Surv(time,status) ~1, lung, subset=(sex==1))
> fit3b<- survreg(Surv(time,status) ~1, lung, subset=(sex==2))
> > fit1
Call:
survreg(formula = Surv(time, status) ~ strata(sex), data = lung)
Coefficients:
(Intercept)
6.06217
Scale:
sex=1 sex=2
0.8167547 0.6533025
Loglik(model)= -1152.5 Loglik(intercept only)= -1152.5
n= 228
> fit2
Call:
survreg(formula = Surv(time, status) ~ strata(sex) + sex, data = lung)
Coefficients:
(Intercept) sex
5.49441 0.3801714
Scale:
sex=1 sex=2
0.8084286 0.6355802
Loglik(model)= -1147.1 Loglik(intercept only)= -1152.5
Chisq= 10.9 on 1 degrees of freedom, p= 0.00096
n= 228
> aeq(fit2$scale, c(fit3a$scale, fit3b$scale), tolerance=tol)
[1] T
> aeq(fit2$loglik[2], (fit3a$loglik + fit3b$loglik)[2], tolerance=tol)
[1] T
> aeq(fit2$coef[1] + 1:2*fit2$coef[2], c(fit3a$coef, fit3b$coef), tolerance=tol)
[1] T
> > #penalized models
> fit1 <- survreg(Surv(time, status) ~ pspline(age, theta=.92)+strata(sex), lung)
> fit2 <- survreg(Surv(time, status) ~ pspline(age, theta=.92)+
+ strata(sex) + sex, lung)
> fit1
Call:
survreg(formula = Surv(time, status) ~ pspline(age, theta = 0.92) + strata(sex),
data = lung)
coef se(coef) se2 Chisq DF p
(Intercept) 6.9036 0.8469 0.5688 66.45 1.00 3.3e-16
pspline(age, theta = 0.92 -0.0124 0.0067 0.0067 3.45 1.00 6.3e-02
pspline(age, theta = 0.92 2.53 2.65 4.0e-01
Scale:
sex=1 sex=2
0.807 0.654
Iterations: 1 outer, 3 Newton-Raphson
Theta= 0.92
Degrees of freedom for terms= 0.5 3.6 2.0
Likelihood ratio test=6.54 on 3.1 df, p=0.0937 n= 228
> fit2
Call:
survreg(formula = Surv(time, status) ~ pspline(age, theta = 0.92) + strata(sex) +
sex, data = lung)
coef se(coef) se2 Chisq DF p
(Intercept) 6.3729 0.84471 0.59118 56.92 1.00 4.5e-14
pspline(age, theta = 0.92 -0.0111 0.00666 0.00666 2.77 1.00 9.6e-02
pspline(age, theta = 0.92 2.46 2.68 4.2e-01
sex 0.3686 0.11711 0.11685 9.91 1.00 1.6e-03
Scale:
sex=1 sex=2
0.8 0.636
Iterations: 1 outer, 4 Newton-Raphson
Theta= 0.92
Degrees of freedom for terms= 0.5 3.7 1.0 2.0
Likelihood ratio test=16.8 on 4.2 df, p=0.00245 n= 228
> > age1 <- ifelse(lung$sex==1, lung$age, mean(lung$age))
> age2 <- ifelse(lung$sex==2, lung$age, mean(lung$age))
> fit3 <- survreg(Surv(time,status) ~ pspline(age1, theta=.92) +
+ pspline(age2, theta=.95) + sex + strata(sex), lung,
+ rel.tol=1e-6)
> fit3a<- survreg(Surv(time,status) ~pspline(age, theta=.92), lung,
+ subset=(sex==1))
> fit3b<- survreg(Surv(time,status) ~pspline(age, theta=.95), lung,
+ subset=(sex==2))
> > # relax the tolerance a little, since the above has lots of parameters
> # I still don't exactly match the second group, but very close
> aeq(fit3$scale, c(fit3a$scale, fit3b$scale), tolerance=tol*10)
[1] "Mean relative difference: 0.001270825"
> aeq(fit3$loglik[2], (fit3a$loglik + fit3b$loglik)[2], tolerance=tol*10)
[1] T
> pred <- predict(fit3)
> aeq(pred[lung$sex==1] , predict(fit3a), tolerance=tol*10)
[1] T
> aeq(pred[lung$sex==2], predict(fit3b), tolerance=tol*10)
[1] "Mean relative difference: 0.01158256"
> > > > > #
> # Some tests using the rat data
> #
> rats <- read.table("../testfrail/data.rats",
+ col.names=c("litter", "rx", "time", "status"))
> > rfitnull <- survreg(Surv(time, status) ~1, rats)
> temp <- rfitnull$scale^2 * pi^2/6
> cat("Effective n =", round(temp*(solve(rfitnull$var))[1,1],1), "\n")
Effective n = 65.8
> > rfit0 <- survreg(Surv(time, status) ~ rx , rats)
> print(rfit0)
Call:
survreg(formula = Surv(time, status) ~ rx, data = rats)
Coefficients:
(Intercept) rx
4.983121 -0.2385013
Scale= 0.2637831
Loglik(model)= -242.3 Loglik(intercept only)= -246.3
Chisq= 8 on 1 degrees of freedom, p= 0.0047
n= 150
> > rfit1 <- survreg(Surv(time, status) ~ rx + factor(litter), rats)
> temp <- rbind(c(rfit0$coef, rfit0$scale), c(rfit1$coef[1:2], rfit1$scale))
> dimnames(temp) <- list(c("rfit0", "rfit1"), c("Intercept", "rx", "scale"))
> temp
Intercept rx scale
rfit0 4.983121 -0.2385013 0.2637831
rfit1 4.902437 -0.2189405 0.2025429
> > > rfit2a <- survreg(Surv(time, status) ~ rx +
+ frailty.gaussian(litter, df=13, sparse=F), rats )
> rfit2b <- survreg(Surv(time, status) ~ rx +
+ frailty.gaussian(litter, df=13, sparse=T), rats )
> > rfit3a <- coxph(Surv(time,status) ~ rx +
+ frailty.gaussian(litter, df=13, sparse=F), rats )
> rfit3b <- coxph(Surv(time,status) ~ rx +
+ frailty(litter, df=13, dist="gauss"), rats)
> > temp <- cbind(rfit2a$coef[3:52], rfit2b$frail, rfit3a$coef[2:51], rfit3b$frail)
> dimnames(temp) <- list(NULL, c("surv","surv.sparse","cox","cox.sparse"))
> pairs(temp)
> apply(temp,2,var)/c(rfit2a$scale, rfit2b$scale, 1,1)^2
surv surv.sparse cox cox.sparse
0.1346009 0.1346009 0.1224049 0.1207863
> apply(temp,2,mean)
surv surv.sparse cox cox.sparse
-7.979728e-19 6.938894e-20 -1.096345e-17 1.054712e-17
> > # The parametric model gives the coefficients less variance for the
> # two fits, for the same df, but the scaled results are similar.
> # 13 df is near to the rmle for the rats
>
> rm(temp, rfit2a, rfit2b, rfit3a, rfit3b, rfitnull, rfit0, rfit1)
> > temp <- matrix(scan("data.mpip", skip=23), ncol=13, byrow=T)
> dimnames(temp) <- list(NULL, c("ved", "angina", "education", "prior.mi",
+ "nyha", "rales", "ef", "ecg", "angina2", "futime",
+ "status", "admit", "betab"))
> > mpip <- data.frame(temp)
> lved <- log(mpip$ved + .02)
> > fit1 <- coxph(Surv(futime, status) ~ pspline(lved) + factor(nyha) +
+ rales + pspline(ef), mpip)
> > temp <- predict(fit1, type="terms", se.fit=T)
> yy <- cbind(temp$fit[,4], temp$fit[,4] + 1.96*temp$se[,4],
+ temp$fit[,4] - 1.96*temp$se[,4])
> index <- order(mpip$ef)
> matplot(mpip$ef[index], yy[index,], type="l", lty=c(1,2,2), col=1)
> title(xlab="Ejection Fraction", ylab="Cox model risk",
+ main="Post-Infarction Survival")
> > fit2 <- coxph(Surv(futime, status) ~ lved + factor(nyha) + rales +
+ pspline(ef, df=0), mpip)
Warning messages:
1: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <-
(df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 ....
2: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <-
(df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 ....
3: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <-
(df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 ....
4: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <-
(df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 ....
5: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <-
(df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 ....
6: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <-
(df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 ....
7: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <-
(df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 ....
8: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <-
(df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 ....
9: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <-
(df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 ....
10: Condition has 764 elements: only the first used in: if(n < df + 2) dfc <-
(df - n) + ((df + 1) * df)/2 - 1 else dfc <- -1 + (df + 1 ....
> temp <- predict(fit2, type="terms", se.fit=T)
> yy <- cbind(temp$fit[,4], temp$fit[,4] + 1.96*temp$se[,4],
+ temp$fit[,4] - 1.96*temp$se[,4])
> matplot(mpip$ef[index], yy[index,], type="l", lty=c(1,2,2), col=1)
> title(xlab="Ejection Fraction", ylab="Cox model risk",
+ main="Post-Infarction Survival, AIC")
> > > fit3 <- survreg(Surv(futime, status) ~ lved + factor(nyha) + rales +
+ pspline(ef, df=2), mpip, dist="lognormal")
> temp <- predict(fit3, type="terms", se.fit=T)
> yy <- cbind(temp$fit[,4], temp$fit[,4] + 1.96*temp$se[,4],
+ temp$fit[,4] - 1.96*temp$se[,4])
> matplot(mpip$ef[index], yy[index,], type="l", lty=c(1,2,2), col=1)
> title(xlab="Ejection Fraction", ylab="Log-normal model predictor",
+ main="Post-Infarction Survival")
> q()
Generated postscript file "testall.ps".
|
