Covariance (Cox regression)

library("survival")
test1 <- data.frame(time=c(4,3.1,1.1,1,2.1,2,3), status=c(1,1,1,0,1,1,0),
                    x=c(0,2,1,1,1,0,0), sex=c(0,0,0,0,1,1,1))
fitc <- coxph(Surv(time, status) ~ x + sex, test1)
dfbeta <- residuals(fitc, type = "dfbeta")
scorer <- residuals(fitc, type = "score")
fitr <- coxph(Surv(time, status) ~ x + sex, test1, robust = TRUE)
fitc$var
t(dfbeta) %*% dfbeta
fitr$var
fitc$var %*% t(scorer) %*% scorer %*% fitc$var


# observed information matrix (Im)
test2 <- test1[order(-test1$time),]
time <- as.matrix(test2$time)
event <- as.matrix(test2$status)
x <- as.matrix(test2[,c("x","sex")])
n <- dim(x)[1]
p <- dim(x)[2]
beta <- as.matrix(fitc$coefficients, nrow = p)

xbar1 <- rep(0, n)
xbar2 <- array(rep(0, p*n), dim = c(p, n))
xbar3 <- array(rep(0, p*p*n), dim = c(p, p, n))
Im <- array(rep(0, p*p), dim = c(p, p))
ll <- 0

for (i in 1:n) {
  exb <- exp(x[i,]%*%beta)
  if (i == 1) {
    xbar1[1] <- exb
    for (r in 1:p) {
      xbar2[r,1] <- x[1,r]*exb
      for (s in 1:r) {
        xbar3[r,s,1] <- x[1,r]*x[1,s]*exb
      }
    }
  } else {
    xbar1[i] <- xbar1[i-1] + exb
    for (r in 1:p) {
      xbar2[r,i] <- xbar2[r,i-1] + x[i,r]*exb
      for (s in 1:r) {
        xbar3[r,s,i] <- xbar3[r,s,i-1] + x[i,r]*x[i,s]*exb
      }
    }
  }

  if (abs(event[i] - 1.0) <= .Machine$double.eps) {
    ll <- ll + x[i,]%*%beta - log(xbar1[i])
    for (r in 1:p) {
      for (s in 1:r) {
        Im[r,s] <- Im[r,s] + xbar3[r,s,i]/xbar1[i] - xbar2[r,i]*xbar2[s,i]/xbar1[i]^2
        if (r != s) {
          Im[s,r] <- Im[r,s]
        }
      }
    }
  }
}

modelvar <- solve(Im)
modelvar
fitc$var