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
|
temp <- matrix(scan("data.mpip", skip=23), ncol=13, byrow=T)
dimnames(temp) <- list(NULL, c('ved', 'angina', 'education', 'prior.mi',
'nyha', 'rales', 'ef', 'ecg', 'angina2', 'futime',
'status', 'admit', 'betab'))
mpip <- data.frame(temp)
lved <- log(mpip$ved + .02)
fit1 <- coxph(Surv(futime, status) ~ pspline(lved) + factor(nyha) +
rales + pspline(ef), mpip)
temp <- predict(fit1, type='terms', se.fit=T)
yy <- cbind(temp$fit[,4], temp$fit[,4] + 1.96*temp$se[,4],
temp$fit[,4] - 1.96*temp$se[,4])
index <- order(mpip$ef)
matplot(mpip$ef[index], yy[index,], type='l', lty=c(1,2,2), col=1)
title(xlab="Ejection Fraction", ylab="Cox model risk",
main="Post-Infarction Survival")
fit2 <- coxph(Surv(futime, status) ~ lved + factor(nyha) + rales +
pspline(ef, df=0), mpip)
temp <- predict(fit2, type='terms', se.fit=T)
yy <- cbind(temp$fit[,4], temp$fit[,4] + 1.96*temp$se[,4],
temp$fit[,4] - 1.96*temp$se[,4])
matplot(mpip$ef[index], yy[index,], type='l', lty=c(1,2,2), col=1)
title(xlab="Ejection Fraction", ylab="Cox model risk",
main="Post-Infarction Survival, AIC")
fit3 <- survreg(Surv(futime, status) ~ lved + factor(nyha) + rales +
pspline(ef, df=2), mpip, dist='lognormal')
temp <- predict(fit3, type='terms', se.fit=T)
yy <- cbind(temp$fit[,4], temp$fit[,4] + 1.96*temp$se[,4],
temp$fit[,4] - 1.96*temp$se[,4])
matplot(mpip$ef[index], yy[index,], type='l', lty=c(1,2,2), col=1)
title(xlab="Ejection Fraction", ylab="Log-normal model predictor",
main="Post-Infarction Survival")
|
