# Calculating fraction of new information in Cox proportional hazard models

**URL:** <https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773>\
**Category:** data analysis\
**Tags:** survival-analysis, r, prediction\
**Created:** [May 18, 2023, 3:35am UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773 "2023-05-18T03:35:37Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ayawad](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/a/5e9695/32.png) [@Ayawad](https://discourse.datamethods.org/u/Ayawad)\
**Post date:** [May 18, 2023, 3:35am UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/1 "2023-05-18T03:35:38Z")

</div>

Hello,

In continuation to the topic we discussed today in our workshp, in regard to the blog post about calculation of the added predictive information by a biomarker:

> **[Statistical Thinking - Statistically Efficient Ways to Quantify Added...](https://www.fharrell.com/post/addvalue/)**
>
> Researchers have used contorted, inefficient, and arbitrary analyses to demonstrated added value in biomarkers, genes, and new lab measurements. Traditional statistical measures have always been up to the task, and are more powerful and more...

The example mentioned in the blog post, calculates the added predictive information of a biomarker in a Binary Logistic Model.  
The example given, tries to assess the added information of a new biomarker (cholesterol) to a base model.

The formula to calculate the added information is:  
**fraction of new information= 1- variance(pre)/variance(post)**

I will attach here the code used to calculate the **variance** of the model with the biomarker:

```auto
g <- lrm(sigdz ~ rcs(age,4) * sex + rcs(choleste,4) + rcs(age,4) %ia%
         rcs(choleste,4), data=acath) # This is the basemodel+biomarker
post <- predict(g, type='fitted') # post-test probability
variance <- round(var(post), 3)

```

> **Summary**
>
> This text will be hidden

I tried to do the same calculations for a Cox proportional hazard model, would like to get feedback on the code/approach if that is possible:

I am attaching a hypothetical example, I am using a different package (that is the one I used at the time, would be happy to learn other ways of doing the calculations using the rms package)

```auto
Model <- coxph(Surv(surv_time, event) ~ biomarker+ variable1+variable2, data =mydata) # The base model+biomarker

# In to the predict the survival probability I am using exp(-expected), ref: https://cran.r-project.org/web/packages/survival/survival.pdf (page 86-87)

model_prediction <- predict(Model , type="expected")  
model_prob <- exp(-model_predicton )

model_variance <- round(var(model_prob), 3)

```

> **Summary**
>
> This text will be hidden

I think my question is, how can I calculate the probability for the outcome in cox model?

I hope I am being clear enough!

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [May 18, 2023, 1:04pm UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/2 "2023-05-18T13:04:59Z")

</div>

I’m not clear that expected survival is the right metric, and why \exp(-x) is used.

A general issue is choosing a scale for the calculation of the variance of predictions. For binary logistic model I’ve seen more expected results by computing on the probability scale, i.e., running the expit transformation on X\hat{\beta}. For the Cox model it is very easy to compute the variance on X\hat{\beta} which is also nice in not requiring you to choose a time. For a probability scale you have to choose a time u and estimate \Pr(Y \geq u | X) for every subject. This is easy with the `rms` package `cph` and `survest` functions and also with `Predict(..., time=)`.

It’s worth also taking the simpler approach of using R^{2}\_\text{adj} using [this](https://hbiostat.org/bib/r2.html).

---

<div class="post-metadata">

**Author:** ![Ayawad](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/a/5e9695/32.png) [@Ayawad](https://discourse.datamethods.org/u/Ayawad)\
**Post date:** [June 21, 2023, 12:49am UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/3 "2023-06-21T00:49:36Z")

</div>

Thank you for the reply.  
I am revisiting the post for a paper I am working on.  
We are trying to compare the performance of two different biomarkers.  
The biomarkers are supposed to be measuring the same things.

The comparison is of the cox proportional hazard models.  
I calculated the traditional metrics (c statisic, NRI, IDI)  
I wanted to calculate the fraction of added information and produce the back to back histogram.

I noticed the link you provided in your reply for R^2 adj is not working, is there a way to link it again?

Thank you

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [June 21, 2023, 12:31pm UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/4 "2023-06-21T12:31:40Z")

</div>

I fixed the link. Please try again.

---

<div class="post-metadata">

**Author:** ![Ayawad](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/a/5e9695/32.png) [@Ayawad](https://discourse.datamethods.org/u/Ayawad)\
**Post date:** [June 22, 2023, 4:00pm UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/5 "2023-06-22T16:00:17Z")

</div>

Thank you! it works.  
I have few more questions to make sure I am applying the metrics correctly, if that is ok  
(I am sorry in advance if the questions are basic)

From my understanding, and from the blog post, there are different metrics that can be used to asses the fraction of new information (other than the traditionally used ones):

1. The liklehood ratio (the frequenst metric)- _I managed to calculate this._

2. The adj R^2 (related to LR ratio) it allows comparison between two or more different models.

3. The variance in pre vs post test probability and visualising with histograms: **this is where I am a bit struggling.**  
I took your advice in using **survest** function to compute the survival probabilities. I was wondering about choosing the time u, how do I make that decision. (do I choose a unified time point for all the subjects, is it the median/mean time?, or do I provide a vector with the follow up time for each subject)

Again, I am sorry if my questions are too basic, I am still learning.

Thank you

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [June 22, 2023, 4:02pm UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/6 "2023-06-22T16:02:53Z")

</div>

There is no unique way of choosing them time so I’d do it for several possible times, separately.

---

<div class="post-metadata">

**Author:** ![Ayawad](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/a/5e9695/32.png) [@Ayawad](https://discourse.datamethods.org/u/Ayawad)\
**Post date:** [June 22, 2023, 5:15pm UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/7 "2023-06-22T17:15:21Z")

</div>

Thank you for clarifying this!

Best,  
Aya

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [June 22, 2023, 6:31pm UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/8 "2023-06-22T18:31:52Z")

</div>

An advantage of sticking to the linear predictor scale is that it is time-invariant. But I seen a logistic regression example where the variance on the logit scale was a bit counter-intuitive.

---

<div class="post-metadata">

**Author:** ![Ayawad](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/a/5e9695/32.png) [@Ayawad](https://discourse.datamethods.org/u/Ayawad)\
**Post date:** [June 22, 2023, 7:57pm UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/9 "2023-06-22T19:57:57Z")

</div>

how can I calculate the variance while sticking to the linear predictor scale?  
h(t|X) = h0(t) \* exp(b1_X1 + b2_ X2 + … + bp\*Xp)

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [June 23, 2023, 1:29pm UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/10 "2023-06-23T13:29:32Z")

</div>

The linear predictor is what’s inside \exp() so doesn’t require t or h(t).

---

<div class="post-metadata">

**Author:** ![JiaqiLi](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/j/74df32/32.png) [@JiaqiLi](https://discourse.datamethods.org/u/JiaqiLi)\
**Post date:** [December 11, 2023, 9:01am UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/11 "2023-12-11T09:01:30Z")

</div>

Hi Prof. Harrell,

Sorry for asking some simple questions, does the “X^β” refer to the β parameter in cox model?  
For example, β for AGE in the following cox model.  
`Model <- cph(Surv(py, d) ~ AGE,data = dt)`

How can I obtain the variance of β? The following code seems not correct.

```auto
predict(Model , type="expected") # for cph

```

In general linear models, it seems that a higher R-squared value is better, such as above 0.5. However, in Cox models, the pseudo R-squared values I calculate through the cph are usually very low, such as 0.2 or 0.1 or even lower. Is a higher pseudo R2 better like that in general linear models,

And, I have another fundamental question. Is the R2 in cph the adjusted R2?

Thank you.

```auto
                       Model Tests Discrimination    
                                              Indexes    
Obs 618 LR chi2 111.51 R2 0.208    
Events 85 d.f. 18 R2(18,618)0.140    
Center 7.9482 Pr(> chi2) 0.0000 R2(18,85)0.667    
                 Score chi2 117.34 Dxy 0.602    
                 Pr(> chi2) 0.0000

```

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [December 11, 2023, 12:46pm UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/12 "2023-12-11T12:46:32Z")

</div>

My preferred pseudo adjusted (for no. of predictors) R^2 is the 0.667 here. This and the D\_{xy} are quite high. Times to events are usually hard to predict. There are two adjusted R^2 in the output. See [here](https://hbiostat.org/bib/r2.html).

To get the variance-covariance matrix of \beta use `vcov(fit object)`. To get the variance of the linear predictor X\beta use `var(predict(fit object))`.

---

<div class="post-metadata">

**Author:** ![JiaqiLi](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/j/74df32/32.png) [@JiaqiLi](https://discourse.datamethods.org/u/JiaqiLi)\
**Post date:** [December 12, 2023, 4:13am UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/13 "2023-12-12T04:13:42Z")

</div>

Thank you very much, Prof. Harrell.

I want to confirm my understanding。 Is the following expression correct?  
R2(18,618) is the `R2(p, n): adjusted for estimating p regression coefficients (non-intercepts) `

~~R2(18,85) is the `R2(p, m): use effective instead of actual sample size, don’t penalize for overfitting`~~  
!!! This was a mistake I made while copying, and Prof. Harrel corrected it. The correct expression is  
R2(18,85) is the `R2(p, m): effective sample size adjusted for estimated regression coefficients` !!!

About the terminology of R2(0.667) that you recommend, is “Pseudo adjusted R-squared based on an effective sample size” an appropriate expression?  
(I want to make sure that I use the correct terminology when referring to the R2)

---

<div class="post-metadata">

**Author:** ![JiaqiLi](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/j/74df32/32.png) [@JiaqiLi](https://discourse.datamethods.org/u/JiaqiLi)\
**Post date:** [December 12, 2023, 4:57am UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/14 "2023-12-12T04:57:02Z")

</div>

> [@f2harrell](#):
>
> variance-covariance

Prof. Harrell,

If possible, could you please help me review my analysis results and confirm whether my interpretations are correct?

I am currently examining whether there is additional value in predicting diabetes by observing the levels of triglycerides in non-fasting and in fasting blood samples.  
I have two datasets—one for non-fasting and another for fasting blood samples.  
I have constructed two simple models to observe the impact of triglyceride levels.

The following is the non-fasting code. The fasting code is omitted due to limited space. By the way, No. of fasting is 5830, which is lower than non fasting (7217).

```auto
discb <- cph(
        Surv(dm_py, diabetes) ~ AGE + SEXC + strat(AREA) +
            BMIndex + # not use rcs() for age and BMI because no apparent nonlinear trends
            rcs(TIMEm, 4) * rcs(GLU, 4) # TIMEm is time after meals(mins)
        ,
        x = TRUE,
        y = TRUE,
        data = dtnfp
    )
# Obs 7217 LR chi2 576.28 R2 0.097    
# Events 819 d.f. 18 R2(18,7217)0.074    
# Center 3.6175 Pr(> chi2) 0.0000 R2(18,819)0.494    
# Score chi2 862.85 Dxy 0.494    
# Pr(> chi2) 0.0000                  

discn <- cph(
        Surv(dm_py, diabetes) ~ AGE + SEXC + strat(AREA) +
            rcs(TIMEm, 4) * rcs(GLU, 4) +
            rcs(TIMEm, 4) * rcs(TGmmol, 4),
        x = TRUE,
        y = TRUE,
        data = dtnfp
    )
# Obs 7217 LR chi2 678.72 R2 0.113    
# Events 819 d.f. 30 R2(30,7217)0.086    
# Center 2.7732 Pr(> chi2) 0.0000 R2(30,819)0.547    
# Score chi2 987.62 Dxy 0.524    
# Pr(> chi2) 0.0000                        

discb
discn

var(predict(discb))
var(predict(discn))

```

| Table. Model discrimination from including serum triglycerides. | | |
| --- | --- | --- |
| | Model 1(age+ sex + bmi +glu) | Model 2 (further + tg) |
| Nonfasting | | |
| Likelihood ratio test statistic | 576.3 | 678.7 |
| fraction of new information | Ref. | 15.1 |
| Variance-covariance | 0.63 | 0.77 |
| fraction of new information | Ref. | 18.2 |
| Pseudo adjusted R-squared | 0.49 | 0.55 |
| Somers’ D | 0.49 | 0.52 |
| Fasting | | |
| Likelihood ratio test statistic | 812.2 | 867.7 |
| fraction of new information | Ref. | 6.4 |
| Variance-covariance | 1.01 | 1.15 |
| fraction of new information | Ref. | 12.2 |
| Pseudo adjusted R-squared | 0.64 | 0.66 |
| Somers’ D | 0.58 | 0.59 |

For the above results, my interpretation is:

1. non-fasting glucose combined with triglycerides provides a better prediction of diabetes.
2. Fasting glucose is sufficient for a good prediction of diabetes, and adding triglycerides has a relatively slight improvement.

R2 in fasting in model 1 is 0.64, which is higher than that (0.49) in non-fasting. Can these R2s be compared?  
Also, can these variance-covariance (0.6/0.7 in non-fasting and 1.0/1.1 in fasting) be compared?

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [December 12, 2023, 2:27pm UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/15 "2023-12-12T14:27:17Z")

</div>

> [@JiaqiLi](#):
>
> I want to confirm my understanding。 Is the following expression correct?  
> R2(18,618) is the `R2(p, n): adjusted for estimating p regression coefficients (non-intercepts) `
> 
> R2(18,85) is the `R2(p, m): use effective instead of actual sample size, don’t penalize for overfitting`
> 
> About the terminology of R2(0.667) that you recommend, is “Pseudo adjusted R-squared based on an effective sample size” an appropriate expression?  
> (I want to make sure that I use the correct terminology when referring to the R2)

That’s all correct except that any R^{2}\_\text{adj} that adjusts for the regression degrees of freedom is penalizing for overfitting. Not as much as AIC would (it has a 2p penalty instead of p), but trying to get an unbiased estimate of R^2 just like the ordinary R^{2}\_\text{adj} does with ordinary linear regression.

---

<div class="post-metadata">

**Author:** ![JiaqiLi](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/j/74df32/32.png) [@JiaqiLi](https://discourse.datamethods.org/u/JiaqiLi)\
**Post date:** [December 13, 2023, 12:19am UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/16 "2023-12-13T00:19:50Z")

</div>

> [@f2harrell](#):
>
> That’s all correct except that any R^{2}_\text{adj}R2adjR^{2}_\text{adj} that adjusts for the regression degrees of freedom is penalizing for overfitting. Not as much as AIC would (it has a 2p2p2p penalty instead of ppp), but trying to get an unbiased estimate of R^2R2R^2 just like the ordinary R^{2}_\text{adj}R2adjR^{2}_\text{adj} does with ordinary linear regression.

I understood my mistake.  
I made an error when copying the explanation for R²(p, m). R²(p, m) represents the effective sample size adjusted for estimated regression coefficients. I mistakenly copied the definition of R²(m). I will add comments on this mistake in my previous response.  
Thank you very much!!

I also want to confirm the following (sorry for being verbose): Is it valid to compare the R2 calculated on different datasets (fasting and nonfasting) for the same model? What should be considered when comparing R2 for model evaluation?

Some metrics, such as AUC, are not suitable for comparison, Likelihood ratio can be used for comparing nested models, but it seems that comparing the same model across different datasets may not be appropriate. R2 is a good metric for comparison in most situation.

Is my understanding correct?"

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [December 13, 2023, 1:02pm UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/17 "2023-12-13T13:02:25Z")

</div>

Comparing any R^2 measure on different datasets is risky, as R^2 will be lower when predictors have a narrower range. The same goes for other measures such as c and D\_{xy}. It is still useful but you have to take the comparisons with a big grain of salt.

---

<div class="post-metadata">

**Author:** ![JiaqiLi](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/j/74df32/32.png) [@JiaqiLi](https://discourse.datamethods.org/u/JiaqiLi)\
**Post date:** [December 14, 2023, 12:27am UTC](https://discourse.datamethods.org/t/calculating-fraction-of-new-information-in-cox-proportional-hazard-models/6773/18 "2023-12-14T00:27:56Z")

</div>

I see. I will avoid comparing them in different datasets. Thank you for your kind assistance. Your comments are very valuable.