This is an excellent question.

The answer, in this particular case, mostly has to do with the outlier time values. (I've

never been convinced that the death at time 999 isn't really a misplaced code for

"missing", actually). If you change the knots used by the spline you can get quite

different values.

For instance, using a smaller data set:

fit1 <- coxph(Surv(tstart, time, status) ~ trt + prior + karno, veteran)

zph1 <- cox.zph(fit1, transform='identity')

plot(zph1[3])

dtime <- unique(veteran$time[veteran$status ==1]) # all of the death times

veteran2 <- survSplit( Surv(time, status) ~ ., data=veteran, cut=dtime)

fit2 <- coxph(Surv(tstart, time, status) ~ trt + prior + karno +

karno:ns(time, df=4), data=veteran2)

tx <- 0:100 * 10 # x positions for plot

ncall <- attr(terms(fit2), "predvars")[[6]]

ty <- eval(ncall, data.frame(time = tx)) %*% coef(fit2)[4:7] + coef(fit2)[3]

lines(tx, ty, col=2)

-------------

Now it looks even worse! The only difference is that the ns() function has chosen a

different set of knots.

The test used by the cox.zph function is based on a score test and is solid. The plot

that it produces uses a smoothed approximation to the variance matrix and is approximate.

So the diagnostic plot will never exactly match an actual fit. In this data set the

outliers exacerbate the issue. To see this try a different time scale.

------------

zph2 <- cox.zph(fit1, transform= sqrt)

plot(zph2[3])

veteran2$stime <- sqrt(veteran2$time)

fit3 <- coxph(Surv(tstart, time, status) ~ trt + prior + karno +

karno:ns(stime, df=4), data=veteran2)

ncall3 <-attr(terms(fit3), "predvars")[[6]]

ty3 <- eval(ncall3, data.frame(stime= sqrt(tx))) %*% coef(fit3)[4:7] + coef(fit3)[3]

lines(sqrt(tx), ty3, col=2)

----------------

The right tail is now better behaved. Eliminating the points >900 makes things even

better behaved.

Terry T.

On 8/8/19 9:07 AM, Ferenci Tamas wrote:

> I was thinking of two possible ways to

> plot a time-varying coefficient in a Cox model.

>

> One is simply to use survival::plot.cox.zph which directly produces a

> beta(t) vs t diagram.

>

> The other is to transform the dataset to counting process format and

> manually include an interaction with time, expanded with spline (to be

> similar to plot.cox.zph). Plotting the coefficient produces the needed

> beta(t) vs t diagram.

>

> I understand that they're slightly different approaches, so I don't

> expect totally identical results, but nevertheless, they approximate

> the very same thing, so I do expect that the results are more or less

> similar.

>

> However:

>

> library( survival )

> library( splines )

>

> data( veteran )

>

> zp <- cox.zph( coxph(Surv(time, status) ~ trt + prior + karno,

> data = veteran ), transform = "identity" )[ 3 ]

>

> veteran3 <- survSplit( Surv(time, status) ~ trt + prior + karno,

> data = veteran, cut = 1:max(veteran$time) )

>

> fit <- coxph(Surv(tstart,time, status) ~ trt + prior + karno +

> karno:ns( time, df = 4 ), data = veteran3 )

> cf <- coef( fit )

> nsvet <- ns( veteran3$time, df = 4 )

>

> plot( zp )

> lines( 0:1000, ns( 0:1000, df = 4, knots = attr( nsvet, "knots" ),

> Boundary.knots = attr( nsvet, "Boundary.knots" ) )%*%cf[

> grep( "karno:ns", names( cf ) ) ] + cf["karno"],

> type = "l", col = "red" )

>

> Where is the mistake? Something must be going on here, because the

> plots are vastly different...

>

> Thank you very much in advance,

> Tamas

[[alternative HTML version deleted]]

______________________________________________

[hidden email] mailing list -- To UNSUBSCRIBE and more, see

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.