# Do-it-yourself Cox require(survival) Data <- subset(pbc, !is.na(trt)) Data <- Data[order(Data$time),] d <- Data$status != 0 n <- length(d) X <- with(Data, cbind(trt, stage, hepato)) X <- scale(X, scale=FALSE) b <- rep(0, ncol(X)) for (i in 1:5) { if (!all(is.finite(b))) stop("Algorithm failed!") eta <- X%*%b haz <- as.numeric(exp(eta)) rsk <- rev(cumsum(rev(haz))) P <- outer(haz, rsk, '/') P[upper.tri(P)] <- 0 W <- -P %*% diag(d) %*% t(P) diag(W) <- diag(P %*% diag(d) %*% t(1-P)) b <- solve(t(X) %*% W %*% X) %*% t(X) %*% (d - P %*% d) + b print(b) } # Same as fit <- coxph(Surv(time, status!=0) ~ trt + stage + hepato, Data) coef(fit) # Now add bilirubin d <- Data$status != 0 n <- length(d) X <- with(Data, cbind(trt, stage, hepato, bili)) X <- scale(X, scale=FALSE) b <- rep(0, ncol(X)) for (i in 1:8) { if (!all(is.finite(b))) stop("Algorithm failed!") eta <- X%*%b haz <- as.numeric(exp(eta)) rsk <- rev(cumsum(rev(haz))) P <- outer(haz, rsk, '/') P[upper.tri(P)] <- 0 W <- -P %*% diag(d) %*% t(P) diag(W) <- diag(P %*% diag(d) %*% t(1-P)) b <- solve(t(X) %*% W %*% X) %*% t(X) %*% (d - P %*% d) + b print(b) } # Now add step-halving d <- Data$status != 0 n <- length(d) X <- with(Data, cbind(trt, stage, hepato, bili)) X <- scale(X, scale=FALSE) b <- rep(0, ncol(X)) for (i in 1:20) { if (!all(is.finite(b))) stop("Algorithm failed!") eta <- X%*%b haz <- as.numeric(exp(eta)) rsk <- rev(cumsum(rev(haz))) P <- outer(haz, rsk, '/') P[upper.tri(P)] <- 0 W <- -P %*% diag(d) %*% t(P) diag(W) <- diag(P %*% diag(d) %*% t(1-P)) bb <- solve(t(X) %*% W %*% X) %*% t(X) %*% (d - P %*% d) + b b <- .5*bb + .5*b print(b) } fit <- coxph(Surv(time, status!=0) ~ trt + stage + hepato + bili, Data) coef(fit) fit$iter # Ties fit <- coxph(Surv(time, status!=0) ~ trt + stage + hepato + bili, pbc) coef(fit) require(xtable) xtable(rbind(as.numeric(b), coef(fit)), digits=5)