I constructed the same cox regression models by using cph in rms package and coxph in R, but when I compared the two models with BIC, I got 4086.559 for coxph, 4114.43 for cph. But the two models both had the same AIC. As I can see the parameters involved in BIC calculation, the two shouldn't be different. Why did that happen?
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
3
As the Wikipedia entry on AIC notes, with $k$ the number of fitted parameters and $n$ the number of observations:
The formula for the Bayesian information criterion (BIC) is similar to the formula for AIC, but with a different penalty for the number of parameters. With AIC the penalty is $2k$, whereas with BIC the penalty is $\ln(n) k$.
So if the AIC is the same between a coxph and a cph model on the same data, the difference must be in the definition of the number of observations, $n$. Should that be the total number of cases, or the number of events? That's the substantive statistical issue behind this question.
It turns out that the nobs() function, used by BIC() to retrieves the number of cases from an object, returns the total number of cases from a cph object, versus the number of events for coxph. A simple example from the survival vignette:
> cfit1.coxph<- coxph(Surv(time, status) ~ age + sex + wt.loss, data=lung)
> cfit1.coxph
Call:
coxph(formula = Surv(time, status) ~ age + sex + wt.loss, data = lung)
coef exp(coef) se(coef) z p
age 0.0200882 1.0202913 0.0096644 2.079 0.0377
sex -0.5210319 0.5939074 0.1743541 -2.988 0.0028
wt.loss 0.0007596 1.0007599 0.0061934 0.123 0.9024
Likelihood ratio test=14.67 on 3 df, p=0.002122
n= 214, number of events= 152
(14 observations deleted due to missingness)
> nobs(cfit1.coxph)
[1] 152
That's 214 cases, 152 events, and nobs() returns the number of events from the coxph object. In contrast:
> cfit1.cph<- cph(Surv(time, status) ~ age + sex + wt.loss, data=lung)
> nobs(cfit1.cph)
[1] 214
So with nobs() you get the total number of cases from a cph object. AICs are the same, BICs higher for cph as in this question:
> AIC(cfit1.coxph)
[1] 1352.112
> AIC(cfit1.cph)
[1] 1352.112
> BIC(cfit1.coxph)
[1] 1361.183
> BIC(cfit1.cph)
[1] 1362.21
The difference in BIC values is (within rounding) what you would expect with 3 parameters fitted:
> BIC(cfit1.coxph)-BIC(cfit1.cph)
[1] -1.026287
> 3*log(nobs(cfit1.coxph))-3*log(nobs(cfit1.cph))
[1] -1.026286
As you should only be doing BIC comparisons on the same data set, either stick to one of coxph() or cph(), or recalculate BIC to use a single nobs() value of choice: do you want case number or event number?
-
$\begingroup$ Thanks so much for your answer! I need to use the cph in rms package to test for the non-linearity of BMI/weight and breast cancer risk. The results showed that BMI could have linear association with breast cancer risk while weight not. So, for BMI, I used coxph package to construct the model---there is no need to include a spline term for BMI, while for weight, I need to include weight as a spline term in the model by using cph. Besides, weight showed better model fit than BMI. So you can see I need to compare the models using coxph package and cph. $\endgroup$– ZhoufengCommented Jul 20, 2021 at 22:11
-
$\begingroup$ Can I recalculate BIC in cph by using the total number of events instead of the total number of observations to make it consistent with the BIC calculation in coxph package? Or is it more reasonable to use AIC rather than BIC when selecting models in this scenario? $\endgroup$– ZhoufengCommented Jul 20, 2021 at 22:15
-
1$\begingroup$ @Zhoufeng the choice between AIC and BIC is essentially how strongly you want to penalize additional predictors. If possible, stick with
cphfor all modeling--cphcan do most of whatcoxphcan be do. If there's somecoxphanalysis that you can't do withcph, then you will be fine re-calculating the BIC for thecphmodel from the number of events. Do be careful, however, that things like missing data aren't confusing your model comparisons: you need to have the exact same cases used for all models after missing data are (somewhat silently) removed. $\endgroup$– EdMCommented Jul 21, 2021 at 15:14