[R] Interpreting the example given by Prof Frank Harrell in {Design} validate.cph

Frank Harrell f.harrell at vanderbilt.edu
Fri Feb 25 15:23:54 CET 2011


Here's the way I would explore this, and some of the code is made more tidy. 
Note that also you could vectorize your simulation.  I have used set.seed
multiple times to make bootstrap samples the same across runs.  -Frank

. . .
if (data[i, 3] == 4) data[i, 5] <- sample(c(0, 1), 1, prob=c(.06,  .94))} 

d <- data.frame(tumor=factor(data[,1]), ecog=factor(data[,2]),
rx=factor(data[,3]), os=data[,4], censor=data[,5])
S <- with(d, Surv(os, censor))

## Check collinearity of rx with other predictors
lrm(rx ~ tumor*ecog, data=d)
## What is the marginal strength of rx (assuming PH)?
cph(S ~ rx, data=d)
## What is partial effect of rx (assuming PH)?
anova(cph(S ~ tumor + ecog + rx, data=d))
## What is combined partial effect of tumor and ecog adjusting for rx?
anova(cph(S ~ tumor + ecog + strat(rx), data=d), tumor, ecog) ## nothing but
noise
## What is their effect not adjusting for rx
cph(S ~ tumor + ecog, data=d)  ## huge

f <- cph(S ~ tumor + ecog, x=TRUE, y=TRUE, surv=TRUE, data=d)
set.seed(1)
validate(f, B=100, dxy=TRUE)
w <- rep(1, 1000)   #  only one stratum, doesn't change model
f <- cph(S ~ tumor + ecog + strat(w), x=TRUE, y=TRUE, surv=TRUE, data=d)
set.seed(1)
validate(f, B=100, dxy=TRUE, u=60)
## identical to last validate except for -Dxy

f <- cph(S ~ tumor + ecog + strat(rx), x=TRUE, y=TRUE, surv=TRUE,
time.inc=60, data=d)
set.seed(1)
validate(f, B=100)  ## no predictive ability
set.seed(1)
validate(f, B=100, dxy=TRUE, u=60)
## Only Dxy indicates some predictive information; large in abs. value
## than model ignoring rx (0.3842 vs. 0.3177)




-----
Frank Harrell
Department of Biostatistics, Vanderbilt University
-- 
View this message in context: http://r.789695.n4.nabble.com/Interpreting-the-example-given-by-Prof-Frank-Harrell-in-Design-validate-cph-tp3316820p3324516.html
Sent from the R help mailing list archive at Nabble.com.



More information about the R-help mailing list