# How to perform a likelihood ratio test in a new dataset？

I understood the purpose of likelihood ratio test. However, I am confused by the table below from a published investigation. Clearly the author compare a nested model (VMNS) and a full model (VMNS + tumor differentiation) in C-index, and likelihood ratio test chi square. Now the question is:

1. I know how to do a likelihood test for two nested models in the training set in R like, anova(coxph1,coxph2), but how can I do it in the validation set? Does the following code right? (though I thought it's wrong...)?

coxph1 <- coxph(Surv(time,status) ~ age, data = datatrain)

coxph2 <- coxph(Surv(time,status) ~ smoke + age, data = datatrain)

pred1 <- predict(coxph1,newdata = datatest)

pred2 <- predict(coxph2,newdata = datatest)

coxph1.test <- coxph(Surv(time,status) ~ pred1, data = datatest)

coxph2.test <- coxph(Surv(time,status) ~ pred2, data = datatest)

anova(coxph1.test,coxph2.test)

2. Is there a function to test C-index for cox model in R (third colum in the table)? • Please add a link to the publication from which you got this table. The original publication might contain hints about how these particular authors did these comparisons, which will help provide an answer to your question.
– EdM
Sep 30, 2022 at 19:11
• Thx for your advice. Here is the link "pubmed.ncbi.nlm.nih.gov/34515782". I believe all statements related were located at the last sentence of the first paragraph of Statistical Analysis Oct 1, 2022 at 3:07

As best as I can tell from the linked publication, the models based on the training set were applied, with the same regression coefficients, to the data in the internal and external validation sets.

The problem with your proposed approach is that the coxph() function will try to re-fit models to pred1 and pred2 on the new data. You would get a new regression coefficient for each model, a slope for how each of pred1 and pred2 is associated with log hazard in the new data set. That's useful in some contexts; for example, it's how the validate() function in Harrell's rms package evaluates the optimism in a model. But that's not what you want.

The trick with likelihoods based on a model with pre-defined regression coefficients is to recognize that you just want to evaluate the (partial) likelihood of the data, given that model. AdamO shows, with code, a very simple way to do this with Cox models:

If you simply want the partial likelihood, why not fool R into giving it to you? Simply initialize beta and allow no iterations, then extract the loglik value from the coxph object.

I don't know if that's precisely how the authors of this study did this, but that's the basic idea. An alternative might be to use offset(pred1) and offset(pred2) to force their regression coefficients to be exactly 1, then extract the log likelihoods.

For the C-index calculations, the authors say that they used tools in the compareC package. I don't have experience with that package. The basic R survival package has a concordance() function that can evaluate the C-index for a model applied to new data.

A few warnings.

First, as the authors note, even Harrell doesn't find the C-index useful for comparisons among models. It's a good measure of discrimination for a single model.

Second, a split into training and validation sets isn't the most efficient use of data when there are so few cases. See this post, for example. Resampling from a combined data set is generally better unless there are tens of thousands of cases.

Third, the authors' use of LASSO might have found a set of predictors that work OK, but the same procedure applied to a new data set might well find a different set of predictors. It's not clear that the authors' use of bootstrapping evaluated the entire model-building process including the LASSO predictor selection. That's another advantage of a full resampling-based validation of the model-building process in this situation.