# 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)