fix a bug and return individual score function under the null

[zhangh12]

1. imat should be a symmetric matrix. Latest version is not.

2. logrank test is equivalent to the score test of cox ph model. inv(null_imat) %*% null_score can be used to compute the covariance between logrank test and other test statistics (e.g. log hazard ratio, or estimates of other models based on the same data). I hope this individual level score under the null can be returned by coxph as well.

More context about 2).

Many estimates can be approximated as inv(hessian) %% (sum of individual scores). The logrank test is the score test based on the Cox PH model under the null. If individual score under the null is available, we can define one-sided logrank test based on S = inv(null hessian) %% (sum of null individual score), and the testing statistic is S' %% inv(null hessian) %% cov(null score) %% inv(null hessian) %% S (two sided) or S^2 diag(inv(null hessian) %% cov(null score) %*% inv(null hessian)) (one-sided).

Imagining that another model (M1) is also fitted on the same data. The estimates can be approximated as inv(hessian) %*% (sum of individual scores). With the null score of the Cox PH model, we can compute the covariance between estimates of M1 and the one-sided logrank statistic.

[therneau]

Before changing, I'd like an example where imat is not symmetric. Also, define what you mean by "individual level score".

[zhangh12]

Before changing, I'd like an example where imat is not symmetric. Also, define what you mean by "individual level score".

It turns out that imat becomes symmetric after the line "sctest = temp; / score test */" when beta and imat return to original scale:

    imat[j][i] *= scale[i]*scale[j];
    imat[i][j] = imat[j][i];           /* <------------- */

So technically this is not a bug.

Score is the first order derivative of log partial likelihood. It can be presented as the sum over samples, I call each of those terms as "individual level" score (i.e. S_i (beta)). Usually resid(..., type = 'score') can return S_i, but evaluate them at beta's estimate. I am requesting returning S_i at beta = null value (in current implementation, null value = initial value) as well.

The logrank test is the score test of Cox PH model under the null. S_i(beta = null) is useful to compute the covariance between logrank statistic and other estimate. For example, we can compute Cov(logrank, log(HR)) through S_i(beta = null) and S_i(beta = estimate)

My codes implement S_i(beta) for Breslow and Efron's methods for handling ties.

[therneau]

The imat in the code is not "technically" correct, it is completely correct. If you were to research the cholesky decomposition more fully (cholesky[235].c chsolve[235].c) you would find that there a good reasons to only work on 1/2 the matrix. That is why the symmetry is restored only after the iteration is finished. Changing that code in the middle of the C loop is unnecessary and potentially dangerous. I think you have jumped to an unwarrarnted conclusion about the desirablilty of symmetry.

With respect to the score residuals at beta =0, you can get them easily now via two calls: fit0 <- coxph(Surv. ........, iter=0) r0 <- resid(fit0, type='score')

The formula you have in the prior comment is not the score residual, by the way. See equation 8 in the validation pdf. The sum is correct, but not the individual terms. The resid routine gives the correct answer.

You are asking that I add a new n x p matrix to the default output of coxph, without any evidence that it will be used by anyone else but yourself. What is the practical utility of the correlation that you propose? (I'm not aware of any papers, but I don't keep up with the literature). It is on the face a very strange statistic, at least as I try to grasp it. Any additions to a fundamental function like coxph first require a debate about the potential utility of such a change. What makes this addition statistically compelling?

Last, any such change requires that there also be an addition to the test suite that validates the results. This means more than code to save results and verify that you get the same answer as last time, I mean validation in the sense of the validation vignette.

[zhangh12]

Hi Terry,

Yes you are right the survival package did it right and I am sorry for inaccuracy. I just notice that survival is a "recommended" package; I have never noticed CRAN has the "Priority" category :-) I have figured out a way to extract that "score" information, and I agree that I don't have enough evidence to justify adding it to the current implementation. I may share the paper once it is accepted for your reference if of interested. Thank you.

Han

The imat in the code is not "technically" correct, it is completely correct. If you were to research the cholesky decomposition more fully (cholesky[235].c chsolve[235].c) you would find that there a good reasons to only work on 1/2 the matrix. That is why the symmetry is restored only after the iteration is finished. Changing that code in the middle of the C loop is unnecessary and potentially dangerous. I think you have jumped to an unwarrarnted conclusion about the desirablilty of symmetry.

With respect to the score residuals at beta =0, you can get them easily now via two calls: fit0 <- coxph(Surv. ........, iter=0) r0 <- resid(fit0, type='score')

The formula you have in the prior comment is not the score residual, by the way. See equation 8 in the validation pdf. The sum is correct, but not the individual terms. The resid routine gives the correct answer.

You are asking that I add a new n x p matrix to the default output of coxph, without any evidence that it will be used by anyone else but yourself. What is the practical utility of the correlation that you propose? (I'm not aware of any papers, but I don't keep up with the literature). It is on the face a very strange statistic, at least as I try to grasp it. Any additions to a fundamental function like coxph first require a debate about the potential utility of such a change. What makes this addition statistically compelling?

Last, any such change requires that there also be an addition to the test suite that validates the results. This means more than code to save results and verify that you get the same answer as last time, I mean validation in the sense of the validation vignette.