@@ -11,9 +11,10 @@ predict.coxph <- function(object, newdata,
1111 type <- match.arg(type )
1212 if (type == " survival" ) {
1313 survival <- TRUE
14- type <- " expected" # this is to stop lots of "or" statements
14+ type <- " expected" # survival and expecte have nearly the same code path
1515 }
1616 else survival <- FALSE
17+ if (type == " expected" ) reference <- " sample" # a common ref is easiest
1718
1819 n <- object $ n
1920 Terms <- object $ terms
@@ -162,8 +163,14 @@ predict.coxph <- function(object, newdata,
162163 pred <- se <- double(nrow(mf2 ))
163164 newx <- newx - rep(object $ means , each = nrow(newx ))
164165 newrisk <- c(exp(newx %*% object $ coef ) + newoffset )
166+ if (ncol(y ) == 3 && survival ) { #
167+ t0 <- unname(min(y [,1 ])) # the start of the survival curve
168+ # simpler is all(newy[,1] == t0), but
169+ # use of all.equal allows for roundoff error in newdata
170+ if (! isTRUE(all.equal(as.vector(newy [,1 ]), rep(t0 , nrow(newy )))))
171+ stop(" predicted survival must be from the start of the curve" )
165172 }
166-
173+ }
167174 survtype <- ifelse(object $ method == ' efron' , 3 ,2 )
168175 for (i in ustrata ) {
169176 indx <- which(oldstrat == i )
@@ -173,66 +180,60 @@ predict.coxph <- function(object, newdata,
173180 afit.n <- length(afit $ time )
174181 if (missing(newdata )) {
175182 # In this case we need se.fit, nothing else
176- j1 <- approx(afit $ time , 1 : afit.n , y [indx ,1 ], method = ' constant' ,
177- f = 0 , yleft = 0 , yright = afit.n )$ y
183+ j1 <- findInterval(y [indx ,1 ], afit $ time )
178184 chaz <- c(0 , afit $ cumhaz )[j1 + 1 ]
179185 varh <- c(0 , cumsum(afit $ varhaz ))[j1 + 1 ]
180186 xbar <- rbind(0 , afit $ xbar )[j1 + 1 ,,drop = F ]
181187 if (ncol(y )== 2 ) {
182188 dt <- (chaz * x [indx ,]) - xbar
183189 se [indx ] <- sqrt(varh + rowSums((dt %*% object $ var ) * dt )) *
184190 risk [indx ]
185- }
191+ }
186192 else {
187- j2 <- approx(afit $ time , 1 : afit.n , y [indx ,2 ], method = ' constant' ,
188- f = 0 , yleft = 0 , yright = afit.n )$ y
189- chaz2 <- c(0 , afit $ cumhaz )[j2 + 1 ]
190- varh2 <- c(0 , cumsum(afit $ varhaz ))[j2 + 1 ]
191- xbar2 <- rbind(0 , afit $ xbar )[j2 + 1 ,,drop = F ]
193+ j2 <- findInterval(y [indx ,2 ], afit $ time )
194+ chaz2 <- c(0 , afit $ cumhaz )[j2 + 1L ]
195+ varh2 <- c(0 , cumsum(afit $ varhaz ))[j2 + 1L ]
196+ xbar2 <- rbind(0 , afit $ xbar )[j2 + 1L ,,drop = F ]
192197 dt <- (chaz * x [indx ,]) - xbar
193198 v1 <- varh + rowSums((dt %*% object $ var ) * dt )
194199 dt2 <- (chaz2 * x [indx ,]) - xbar2
195200 v2 <- varh2 + rowSums((dt2 %*% object $ var ) * dt2 )
196201 se [indx ] <- sqrt(v2 - v1 )* risk [indx ]
197- }
198202 }
203+ }
199204
200205 else {
201206 # there is new data
202207 use.x <- TRUE
203208 indx2 <- which(newstrat == i )
204- j1 <- approx(afit $ time , 1 : afit.n , newy [indx2 ,1 ],
205- method = ' constant' , f = 0 , yleft = 0 , yright = afit.n )$ y
209+ j1 <- findInterval(newy [indx2 ,1 ], afit $ time )
206210 chaz <- c(0 , afit $ cumhaz )[j1 + 1 ]
207211 pred [indx2 ] <- chaz * newrisk [indx2 ]
208212 if (se.fit ) {
209213 varh <- c(0 , cumsum(afit $ varhaz ))[j1 + 1 ]
210214 xbar <- rbind(0 , afit $ xbar )[j1 + 1 ,,drop = F ]
211- }
215+ }
212216 if (ncol(y )== 2 ) {
213217 if (se.fit ) {
214218 dt <- (chaz * newx [indx2 ,]) - xbar
215219 se [indx2 ] <- sqrt(varh + rowSums((dt %*% object $ var ) * dt )) *
216220 newrisk [indx2 ]
217- }
218221 }
222+ }
219223 else {
220- j2 <- approx(afit $ time , 1 : afit.n , newy [indx2 ,2 ],
221- method = ' constant' , f = 0 , yleft = 0 , yright = afit.n )$ y
222- chaz2 <- approx(- afit $ time , afit $ cumhaz , - newy [indx2 ,2 ],
223- method = " constant" , rule = 2 , f = 0 )$ y
224- chaz2 <- c(0 , afit $ cumhaz )[j2 + 1 ]
224+ j2 <- findInterval(newy [indx2 ,2 ], afit $ time )
225+ chaz2 <- c(0 , afit $ cumhaz )[j2 + 1L ]
225226 pred [indx2 ] <- (chaz2 - chaz ) * newrisk [indx2 ]
226-
227+
227228 if (se.fit ) {
228- varh2 <- c(0 , cumsum(afit $ varhaz ))[j1 + 1 ]
229- xbar2 <- rbind(0 , afit $ xbar )[j1 + 1 ,,drop = F ]
229+ varh2 <- c(0 , cumsum(afit $ varhaz ))[j2 + 1L ]
230+ xbar2 <- rbind(0 , afit $ xbar )[j2 + 1L ,,drop = F ]
230231 dt <- (chaz * newx [indx2 ,]) - xbar
231232 dt2 <- (chaz2 * newx [indx2 ,]) - xbar2
232233
233234 v2 <- varh2 + rowSums((dt2 %*% object $ var ) * dt2 )
234235 v1 <- varh + rowSums((dt %*% object $ var ) * dt )
235- se [indx2 ] <- sqrt(v2 - v1 )* risk [indx2 ]
236+ se [indx2 ] <- sqrt(v2 - v1 )* newrisk [indx2 ]
236237 }
237238 }
238239 }
0 commit comments