Estimating Cox Model parameters with and without coxph function from survival package

Reading data


library(survival)
data(veteran)
time.var <- veteran[,3] 
status.var <- veteran[,4] 
covariate.matrix <- veteran[,c(1,5:8)] 
covariate.names <- names(veteran[,c(1,5:8)])
data <- cbind.data.frame(time.var,status.var,covariate.matrix)

head(data)

  time.var status.var trt karno diagtime age prior
1       72          1   1    60        7  69     0
2      411          1   1    70        5  64    10
3      228          1   1    60        3  38     0
4      126          1   1    60        9  63    10
5      118          1   1    70       11  65    10
6       10          1   1    20        5  49     0

Writing log partial-likelihood for optimization (very simple: it can't not handle tied observations)


cox.LPL <- function(init.param, time.var, status.var,covariate.matrix){
 options(warn = -1)
 time <- time.var
 status <- status.var
 x <- covariate.matrix
 data <- cbind.data.frame(time,status,x)
 data <- data[with(data, order(-time)), ]
 risk.score <- c(data.matrix(data[,-c(1:2)]) %*% init.param)
 PL1 <- sum(risk.score*data$status)
 PL2 <- -sum(log(cumsum(exp(risk.score)))[data$status==1])
 return(-(PL1+PL2))
}


Test run with data to obtain negative log likelihood value


cox.LPL(rep(0,ncol(covariate.matrix)),time.var,status.var,covariate.matrix)


Optimization and coxph function values being compared after running with same dataset


res1 <- nlm(function(x) cox.LPL(x,time.var,status.var,covariate.matrix), rep(0,ncol(covariate.matrix)))$estimate


res2 <- coef(coxph(Surv(time.var,status.var)~data[,3]+data[,4]+data[,5]+data[,6]+data[,7])) # may need to change the numbers depending on how many covariates are there in the dataset that you want to include in the model


res <- round(cbind(res1,res2),5)
dimnames(res) <- list(covariate.names, c("LogPL.optimization", "coxph"))


res



         LogPL.optimization    coxph
trt                 0.19768  0.19305
karno              -0.03409 -0.03408
diagtime            0.00152  0.00172
age                -0.00393 -0.00388
prior              -0.00758 -0.00776

Acknowledgement: original package, teaching link