Cox 比例风险中的偏对数似然函数由以下公式给出 其中是我们观察到的事件的观察次数（通常有次观察，因此次观察被删失）并且是时间的风险集，定义为：。

$${}_{p}\ell(\beta) = \sum\limits_{i=1}^{K}X_i'\beta - \sum\limits_{i=1}^{K}\log\Big(\sum\limits_{l\in \mathscr{R}(t_i)}^{}e^{X_l'\beta}\Big),$$ $K$$n$$K-n$$\mathscr{R}(t_i)$$t_i$$\mathscr{R}(t_i) = : \{X_j: t_j >= t_i, j = 1, \dots, n \}$

我正在尝试为给定的向量和输入数据集实现计算此部分对数似然函数的函数。我认为首先按观察到的时间对数据进行排序是很聪明的，这样较高的行数表示较高的生存时间。对于这种形式的数据，我准备了 2 个解释变量的实现：$\beta$

full_cox_loglik <- function(beta1, beta2, x1, x2, censored){
  sum(rev(censored)*(beta1*rev(x1) + beta2*rev(x2) -
                       log(cumsum(exp(beta1*rev(x1) + beta2*rev(x2))))))
}

其中beta1是变量的系数x1，beta2是变量的系数，x2并且censored是指示观察是否有事件（然后）或被删失（然后是）的向量。$1$$0$

小心我还准备了第二个更长的实现，它通常适用于每个维度——它也以排序格式获取数据。这假设dCox是一个具有解释变量和删失列的数据框，beta是一个系数向量，并status_number指示具有有关删失信息的列的数量：

library(foreach)
partial_coxph_loglik <- function (dCox, beta, status_number) {
  n <- nrow(dCox)
  foreach(i=1:n) %dopar% {
    sum(dCox[i,status_number]*(dCox[i,-status_number]*beta))
  } %>% unlist -> part1


  foreach(i=1:n) %dopar% {
    exp(sum(dCox[i,-status_number]*beta))
  } %>% unlist -> part2


  foreach(i=1:n) %dopar% {
    part1[i] - dCox[i, status_number]*(log(sum(part2[i:n])))
  } %>% unlist -> part3

  sum(part3)
}

部分对数似然在哪里更易读（对于该实现）：

$${}_{p}\ell(\beta) = \sum\limits_{i=1}^{K} \text{part1}_i - \sum\limits_{i=1}^{K}\log\Big(\text{part2}_i\Big).$$

对于这个实现，我试图计算 Cox 比例模型的部分对数似然值，这些数据是从设置为 的真实参数生成的。$\beta$beta=c(2,2)

但是我收到的结果告诉我，部分对数似然的最大值不在这一点上，beta=c(2,2)而是离这一点很远。

可以为来自 Weibull 分布的 Cox 比例风险模型准备模拟生存数据，并在此处解释方法https://stats.stackexchange.com/a/135129/49793。基于该解决方案的类似实现如下（适用于 2 个解释变量）：

set.seed(456)
dataCox <- function(N, lambda, rho, x, beta, censRate){

  # real Weibull times
  u <- runif(N)
  Treal <- (- log(u) / (lambda * exp(x %*% beta)))^(1 / rho)

  # censoring times
  Censoring <- rexp(N, censRate)

  # follow-up times and event indicators
  time <- pmin(Treal, Censoring)
  status <- as.numeric(Treal <= Censoring)

  # data set
  data.frame(id=1:N, time=time, status=status, x=x)
}

x <- matrix(sample(0:1, size = 2000, replace = TRUE), ncol = 2)

dataCox(10^3, lambda = 5, rho = 1.5, x, beta = c(2,2), censRate = 0.2) -> dCox

因此，对于这个实现和模拟数据，我计算了betafrom c(0,0)to的部分对数似然值，c(2,2)并收到了这样的结果：

library(dplyr)
dCox %>%
  dplyr::arrange(time) -> dCoxArr
beta1 <- seq(0,2,0.05)
beta2 <- seq(0,2,0.05)
res <- numeric(length(beta1))
for(i in 1:length(beta1)){
  full_cox_loglik(beta1[i], beta2[i], dCoxArr$x.1,
                  dCoxArr$x.2, dCoxArr$status ) -> res[i]

}

library(ggplot2)
qplot(beta1, res)

res2 <- numeric(length(beta1))
for(i in 1:length(beta1)){
  partial_coxph_loglik(dCoxArr[, c(4,5,3)],
                       c(beta1[i],beta2[i]), 3 ) -> res2[i]
  cat("\r", i, "\r")
}


qplot(beta1, res2)

看起来部分对数似然函数的最大值不在beta = c(2,2)我生成数据的点。所以现在出现了我的问题？我在哪里犯错了？在数据生成中？在部分对数似然实现中？还是别的地方？