library(survival)
if(!requireNamespace('emplik', quietly=TRUE)){
    install.packages('emplik', dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(emplik)

data(myeloma)
myeloma = data.frame(myeloma)
colnames(myeloma) = c('time', 'vstatus', 'logBUN', 'HGB', 'platelet', 'age', 'logWBC', 'FRAC', 'logPBM', 'protein', 'calcium')
print(head(myeloma))

res = coxph(Surv(time, vstatus) ~ ., data=myeloma)
print(res)

## Lets find any outliers by looking at the residuals:
##
par(mfrow=c(2, 1))
plot(res$residuals, type='p', pch=19, xlab='sample index', ylab='residual', main='martingale residuals')
abline(v=48, col='red')
grid()
plot(resid(res, type='deviance'), type='p', pch=19, xlab='sample index', ylab='residual', main='deviance residuals')
grid()
par(mfrow=c(1, 1))

## Print this outlier:
##
print('Outlying patient:')
print(myeloma[48,])
col_means = colMeans(myeloma)
col_stds = apply(myeloma, 2, sd)
t_stat = (myeloma[48, ] - col_means)/col_stds
print('Outlying patients T-stat:')
print(t_stat)

print('Range of time:')
print(range(myeloma$time))