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
|
R version 2.14.0 Under development (unstable) (2011-04-10 r55401)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> options(na.action=na.exclude) # preserve missings
> options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
> library(survival)
Loading required package: splines
>
> #
> # Fit the models found in Wei et. al.
> #
> wfit <- coxph(Surv(stop, event) ~ (rx + size + number)* strata(enum) +
+ cluster(id), bladder, method='breslow')
> wfit
Call:
coxph(formula = Surv(stop, event) ~ (rx + size + number) * strata(enum) +
cluster(id), data = bladder, method = "breslow")
coef exp(coef) se(coef) robust se z p
rx -0.5176 0.596 0.3158 0.3075 -1.683 0.0920
size 0.0679 1.070 0.1012 0.0853 0.796 0.4300
number 0.2360 1.266 0.0761 0.0721 3.274 0.0011
rx:strata(enum)enum=2 -0.1018 0.903 0.5043 0.3265 -0.312 0.7600
rx:strata(enum)enum=3 -0.1823 0.833 0.5579 0.3916 -0.465 0.6400
rx:strata(enum)enum=4 -0.1332 0.875 0.6581 0.4968 -0.268 0.7900
size:strata(enum)enum=2 -0.1440 0.866 0.1680 0.1119 -1.287 0.2000
size:strata(enum)enum=3 -0.2792 0.756 0.2086 0.1511 -1.847 0.0650
size:strata(enum)enum=4 -0.2711 0.763 0.2515 0.1856 -1.460 0.1400
number:strata(enum)enum=2 -0.0984 0.906 0.1193 0.1144 -0.861 0.3900
number:strata(enum)enum=3 -0.0662 0.936 0.1298 0.1167 -0.567 0.5700
number:strata(enum)enum=4 0.0928 1.097 0.1466 0.1175 0.790 0.4300
Likelihood ratio test=29.4 on 12 df, p=0.00344 n= 340, number of events= 112
>
> # Check the rx coefs versus Wei, et al, JASA 1989
> rx <- c(1,4,5,6) # the treatment coefs above
> cmat <- diag(4); cmat[1,] <- 1; #contrast matrix
> wfit$coef[rx] %*% cmat # the coefs in their paper (table 5)
[,1] [,2] [,3] [,4]
[1,] -0.5176209 -0.6194404 -0.6998771 -0.6507935
> t(cmat) %*% wfit$var[rx,rx] %*% cmat # var matrix (eqn 3.2)
[,1] [,2] [,3] [,4]
[1,] 0.09455501 0.06017669 0.05677331 0.0437777
[2,] 0.06017669 0.13242834 0.13011557 0.1160420
[3,] 0.05677331 0.13011557 0.17235879 0.1590865
[4,] 0.04377770 0.11604200 0.15908650 0.2398112
>
> # Anderson-Gill fit
> fita <- coxph(Surv(start, stop, event) ~ rx + size + number + cluster(id),
+ bladder2, method='breslow')
> summary(fita)
Call:
coxph(formula = Surv(start, stop, event) ~ rx + size + number +
cluster(id), data = bladder2, method = "breslow")
n= 178, number of events= 112
coef exp(coef) se(coef) robust se z Pr(>|z|)
rx -0.45979 0.63142 0.19996 0.25801 -1.782 0.07474 .
size -0.04256 0.95833 0.06903 0.07555 -0.563 0.57317
number 0.17164 1.18726 0.04733 0.06131 2.799 0.00512 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rx 0.6314 1.5837 0.3808 1.047
size 0.9583 1.0435 0.8264 1.111
number 1.1873 0.8423 1.0528 1.339
Concordance= 0.634 (se = 0.03 )
Rsquare= 0.09 (max possible= 0.994 )
Likelihood ratio test= 16.77 on 3 df, p=0.000787
Wald test = 11.76 on 3 df, p=0.008256
Score (logrank) test = 18.57 on 3 df, p=0.0003355, Robust = 11.44 p=0.009588
(Note: the likelihood ratio and score tests assume independence of
observations within a cluster, the Wald and robust score tests do not).
>
> # Prentice fits. Their model 1 a and b are the same
> fit1p <- coxph(Surv(stop, event) ~ rx + size + number, bladder2,
+ subset=(enum==1), method='breslow')
> fit2pa <- coxph(Surv(stop, event) ~ rx + size + number, bladder2,
+ subset=(enum==2), method='breslow')
> fit2pb <- coxph(Surv(stop-start, event) ~ rx + size + number, bladder2,
+ subset=(enum==2), method='breslow')
> fit3pa <- coxph(Surv(stop, event) ~ rx + size + number, bladder2,
+ subset=(enum==3), method='breslow')
> #and etc.
> fit1p
Call:
coxph(formula = Surv(stop, event) ~ rx + size + number, data = bladder2,
subset = (enum == 1), method = "breslow")
coef exp(coef) se(coef) z p
rx -0.5176 0.596 0.3158 -1.639 0.1000
size 0.0679 1.070 0.1012 0.671 0.5000
number 0.2360 1.266 0.0761 3.102 0.0019
Likelihood ratio test=9.66 on 3 df, p=0.0216 n= 85, number of events= 47
> fit2pa
Call:
coxph(formula = Surv(stop, event) ~ rx + size + number, data = bladder2,
subset = (enum == 2), method = "breslow")
coef exp(coef) se(coef) z p
rx -0.42421 0.654 0.4022 -1.0547 0.29
size -0.12503 0.882 0.1171 -1.0679 0.29
number 0.00199 1.002 0.0938 0.0212 0.98
Likelihood ratio test=2.02 on 3 df, p=0.569 n= 46, number of events= 29
> fit2pb
Call:
coxph(formula = Surv(stop - start, event) ~ rx + size + number,
data = bladder2, subset = (enum == 2), method = "breslow")
coef exp(coef) se(coef) z p
rx -0.25911 0.772 0.4051 -0.6396 0.52
size -0.11636 0.890 0.1192 -0.9759 0.33
number -0.00571 0.994 0.0967 -0.0591 0.95
Likelihood ratio test=1.27 on 3 df, p=0.735 n= 46, number of events= 29
> fit3pa
Call:
coxph(formula = Surv(stop, event) ~ rx + size + number, data = bladder2,
subset = (enum == 3), method = "breslow")
coef exp(coef) se(coef) z p
rx -0.8985 0.407 0.554 -1.623 0.10
size 0.0850 1.089 0.209 0.407 0.68
number -0.0172 0.983 0.128 -0.134 0.89
Likelihood ratio test=4.16 on 3 df, p=0.245 n= 27, number of events= 22
> rm(rx, cmat, wfit, fita, fit1p, fit2pa, fit2pb, fit3pa)
>
|
