# Do-it-yourself Cox
library(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)

# Taking a closer look at the score function
eta <- rep(0, nrow(X))
haz <- as.numeric(exp(eta))
rsk <- rev(cumsum(rev(haz)))
P0 <- outer(haz, rsk, '/')
P0[upper.tri(P0)] <- 0
P0[301:312,301:312]
P[301:312,301:312]
Data[306:312, colnames(X)]
(P %*% d)[301:312]
PP0 <- P0[,which(d==1)]
PP <- P[,which(d==1)]
PP[290:312,140:144]

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

# Ties
fit <- coxph(Surv(time, status!=0) ~ trt + stage + hepato + bili, pbc)
coef(fit)
rbind(as.numeric(b), coef(fit))