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
|
survfit.cph <- function(formula, newdata, se.fit=TRUE, conf.int=.95,
individual=FALSE, type=NULL, vartype=NULL,
conf.type=c('log', 'log-log', 'plain', 'none'),
id, ...) {
object <- formula
Call <- match.call()
Call[[1]] <- as.name("survfit") ## nicer output for the user
censor <- TRUE
ftype <- object$type
if (! length(ftype)) {
## Use the appropriate one from the model
w <- c("exact", "breslow", "efron")
survtype <- match(object$method, w)
}
else {
w <- c("kalbfleisch-prentice", "aalen", "efron",
"kaplan-meier", "breslow", "fleming-harrington",
"greenwood", "tsiatis", "exact")
survtype <- match(match.arg(type, w), w)
survtype <- c(1,2,3,1,2,3,2,2,1)[survtype]
}
vartype <- if(! length(vartype)) survtype
else {
w <- c("greenwood", "aalen", "efron", "tsiatis")
vt <- match(match.arg(vartype, w), w)
if(vt == 4) 2 else vt
}
if (! se.fit) conf.type <- "none"
else conf.type <- match.arg(conf.type)
xpres <- length(object$means) > 0
y <- object[['y']] # require exact name match
type <- attr(y, 'type')
if(! length(y)) stop('must use y=TRUE with fit')
if(xpres) {
X <- object[['x']]
if(! length(X)) stop('must use x=TRUE with fit')
n <- nrow(X)
xcenter <- object$means
X <- X - rep(xcenter, rep.int(n, ncol(X)))
}
else {
n <- nrow(y)
X <- matrix(0, nrow=n, ncol=1)
}
strata <- object$strata ###
strata.pres <- length(strata) > 0
if(! length(strata)) strata <- rep(0, n)
offset <- object$offset
if(! length(offset)) offset <- rep(0., n)
weights <- object$weights
if(! length(weights)) weights <- rep(1., n)
missid <- missing(id)
if (! missid) individual <- TRUE
else if (missid && individual) id <- rep(0, n)
else id <- NULL
if (individual && type != "counting")
stop("The individual option is only valid for start-stop data")
## Compute confidence limits for survival based on -log survival,
## constraining to be in [0,1]; d = std.error of cum hazard * z value
ciupper <- function(surv, d) ifelse(surv==0, 0, pmin(1, surv*exp(d)))
cilower <- function(surv, d) ifelse(surv==0, 0, surv*exp(-d))
risk <- rep(exp(object$linear.predictors), length=n)
## need to center offset??
## coxph.fit centered offset inside linear predictors
if(missing(newdata)) {
X2 <- if(xpres) matrix(0., nrow=1, ncol=ncol(X)) else
matrix(0., nrow=1, ncol=1)
rq <- ro <- NULL
newrisk <- 1
}
else {
if (length(object$frail))
stop("The newdata argument is not supported for sparse frailty terms")
X2 <- predictrms(object, newdata, type='x', expand.na=FALSE)
## result with type='x' has attributes strata and offset which may be NULL
rq <- attr(X2, 'strata')
ro <- attr(X2, 'offset')
n2 <- nrow(X2)
if(length(rq) && any(levels(rq) %nin% levels(strata)))
stop('new dataset has strata levels not found in the original')
if(! length(rq)) rq <- rep(1, n2)
ro <- if(length(ro)) ro - mean(offset) else rep(0., n2)
X2 <- X2 - rep(xcenter, rep.int(n2, ncol(X2)))
newrisk <- exp(matxv(X2, object$coefficients) + ro)
}
y2 <- NULL
if (individual) {
if(missing(newdata))
stop("The newdata argument must be present when individual=TRUE")
isS <- sapply(newdata, is.Surv)
if(sum(isS) != 1)
stop("newdata must contain exactly one Surv object when individual=TRUE")
y2 <- newdata[[which(isS)]]
warning('some aspects of individual=TRUE not yet implemented. Try survfit.coxph.')
}
g <- survfitcoxph.fit(y, X, weights, X2, risk, newrisk, strata,
se.fit, survtype, vartype,
if(length(object$var)) object$var else
matrix(0, nrow=1, ncol=1),
id=id, y2=y2, strata2=rq)
if(strata.pres) {
if (is.matrix(g$surv)) {
nr <- nrow(g$surv) #a vector if newdata had only 1 row
indx1 <- split(1:nr, rep(1:length(g$strata), g$strata))
rows <- indx1[as.numeric(rq)] #the rows for each curve
indx2 <- unlist(rows) #index for time, n.risk, n.event, n.censor
indx3 <- as.integer(rq) #index for n and strata
for(i in 2:length(rows)) rows[[i]] <- rows[[i]]+ (i-1)*nr #linear subscript
indx4 <- unlist(rows) #index for surv and std.err
temp <- g$strata[indx3]
names(temp) <- row.names(X2) #row.names(mf2)
new <- list(n = g$n[indx3],
time= g$time[indx2],
n.risk= g$n.risk[indx2],
n.event=g$n.event[indx2],
n.censor=g$n.censor[indx2],
strata = temp,
surv= g$surv[indx4],
cumhaz = g$cumhaz[indx4])
if (se.fit) new$std.err <- g$std.err[indx4]
g <- new
}
}
## Insert type so that survfit.cph produces object like survfit.coxph
g$type <- type
if (! censor) {
kfun <- function(x, keep) {
if (is.matrix(x)) x[keep,, drop=FALSE]
else if (length(x) == length(keep)) x[keep] else x
}
keep <- g$n.event > 0
if(length(g$strata)) {
w <- factor(rep(names(g$strata), g$strata), levels=names(g$strata))
g$strata <- c(table(w[keep]))
}
g <- lapply(g, kfun, keep)
}
if (se.fit) {
zval <- qnorm(1 - (1 - conf.int)/2, 0, 1)
if (conf.type=='plain') {
u <- g$surv + zval* g$std.err * g$surv
z <- g$surv - zval* g$std.err * g$surv
g <- c(g, list(upper=pmin(u,1), lower=pmax(z,0),
conf.type='plain', conf.int=conf.int))
}
if (conf.type=='log')
g <- c(g, list(upper=ciupper(g$surv, zval * g$std.err),
lower=cilower(g$surv, zval * g$std.err),
conf.type='log', conf.int=conf.int))
if (conf.type=='log-log') {
who <- (g$surv==0 | g$surv==1) #special cases
xx <- ifelse(who, .1, g$surv) #avoid some "log(0)" messages
u <- exp(-exp(log(-log(xx)) + zval * g$std.err/log(xx)))
u <- ifelse(who, g$surv + 0 * g$std.err, u)
z <- exp(-exp(log(-log(xx)) - zval*g$std.err/log(xx)))
z <- ifelse(who, g$surv + 0 * g$std.err, z)
g <- c(g, list(upper=u, lower=z,
conf.type='log-log', conf.int=conf.int))
}
}
g$requested.strata <- rq
g$call <- Call
class(g) <- c('survfit.cph', 'survfit.cox', 'survfit')
g
}
|
