7
$\begingroup$

I try to figure out how to describe my continuous variable. Unfortunately, I did not understand all of the statistics. I would really appreciate I you guys would help me out here. To better illustrate my problem, I wrote the following example.

library(rms)
library(survival)

data(pbc)
d <- pbc
rm(pbc, pbcseq)
d$status <- ifelse(d$status != 0, 1, 0)

dd = datadist(d)
options(datadist='dd')

# linear model
f1 <- cph(Surv(time, status) ~  albumin, data=d)
p1 <- Predict(f1, fun=exp)
(a1 <- anova(f1))
Function(f1)
plot(p1, anova=a1, pval=TRUE, ylab="Hazard Ratio")

# rcs model
f2 <- cph(Surv(time, status) ~  rcs(albumin, 4), data=d)
p2 <- Predict(f2, fun=exp)
(a2 <- anova(f2))
Function(f2)
plot(p2, anova=a2, pval=TRUE, ylab="Hazard Ratio")

# minimal CI width
p1$diff <- p1$upper-p1$lower
min(p1$diff) # = 0.002321521
p1[which(p1$diff==min(p1$diff)),]$albumin # = 3.494002
describe(d$albumin) # mean = 3.497

p2$diff <- p2$upper-p2$lower
min(p2$diff) # = 0.2039817
p2[which(p2$diff==min(p2$diff)),]$albumin # = 3.502447
describe(d$albumin) # mean = 3.497

# both models in a single figure
p <- rbind(linear.model=p1, rcs.model=p2)
library(ggplot2)
df <- data.frame(albumin=p$albumin, yhat=p$yhat, lower=p$lower, upper=p$upper, predictor=p$.set.)
(g <- ggplot(data=df, aes(x=albumin, y=yhat, group=predictor, color=predictor)) + geom_line(size=1))
(g <- g + geom_ribbon(data=df, aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0))
(g <- g + theme_bw())
(g <- g + xlab("Albumin"))
(g <- g + ylab("Hazard Ratio"))
(g <- g + theme(axis.line = element_line(color='black', size=1)))
(g <- g + theme(axis.ticks = element_line(color='black', size=1)))
(g <- g + theme( plot.background = element_blank() ))
(g <- g + theme( panel.grid.minor = element_blank() ))
(g <- g + theme( panel.border = element_blank() ))
  1. Why shows the plot of the linear model (p1) not a straight line?
  2. How can I plot the models f1 and f2 in the same figure?
  3. How can I compare the models f1 and f2 to investigate which models fits the data better? ... like anova() for coxph in the survival package
  4. Why is the minimal CI width near the mean of albumin more pronounce in the linear (f1) model?
  5. What does the P value in the plots mean? How do I have to interpret the output of anova(...)

pkot of f1

plot of f2

combined plot

Update #1 Following the answer from Harrell, I updated the code above showing how to combine spline plots of two predictors in a single figure. One last question: How can I compare the two rms models like anova(m1, m2) of the survival package as shown below?

> m1 <- coxph(Surv(time, status) ~ albumin, data=d)
> m2 <- coxph(Surv(time, status) ~ pspline(albumin), data=d)
> anova(m1, m2) # compare models
Analysis of Deviance Table
 Cox model: response is  Surv(time, status)
 Model 1: ~ albumin
 Model 2: ~ pspline(albumin)
   loglik  Chisq Df P(>|Chi|)
1 -975.61                    
2 -973.26 4.6983 11    0.9449
> summary(m1)
Call:
coxph(formula = Surv(time, status) ~ albumin, data = d)

  n= 418, number of events= 186 

           coef exp(coef) se(coef)      z Pr(>|z|)    
albumin -1.4695    0.2300   0.1714 -8.574   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

        exp(coef) exp(-coef) lower .95 upper .95
albumin      0.23      4.347    0.1644    0.3219

Concordance= 0.688  (se = 0.023 )
Rsquare= 0.147   (max possible= 0.992 )
Likelihood ratio test= 66.6  on 1 df,   p=3.331e-16
Wald test            = 73.51  on 1 df,   p=0
Score (logrank) test = 72.38  on 1 df,   p=0

UPDATE #2 I think I just answered my "one last question" by myself (see below). I hope this does not show correct accidentally. I would think that I can compare models from cph and coxph that way, can't I? Is the way of calculating the degrees of freedom df correct?

> # using coxph from survival
> m1 <- coxph(Surv(time, status) ~  albumin, data=d)
> m2 <- coxph(Surv(time, status) ~  albumin + age, data=d)
> # loglik = a vector of length 2 containing the log-likelihood with the initial values and with the final values of the coefficients.
> m1$loglik[2]
[1] -975.6126
> m2$loglik[2]
[1] -973.2272
> (df <- abs(length(m1$coefficients) - length(m2$coefficients)))
[1] 1
> (LR <- 2 * (m2$loglik[2] - m1$loglik[2]))
[1] 4.770787
> pchisq(LR, df, lower=FALSE)
[1] 0.02894659
> anova(m2, m1)
Analysis of Deviance Table
 Cox model: response is  Surv(time, status)
 Model 1: ~ albumin + age
 Model 2: ~ albumin
   loglik  Chisq Df P(>|Chi|)  
1 -973.23                      
2 -975.61 4.7708  1   0.02895 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> m1 <- cph(Surv(time, status) ~  albumin, data=d)
> m2 <- cph(Surv(time, status) ~  albumin + age, data=d)
> # loglik = a vector of length 2 containing the log-likelihood with the initial values and with the final values of the coefficients.
> m1$loglik[2]
[1] -975.6126
> m2$loglik[2]
[1] -973.2272
> (df <- abs(length(m1$coefficients) - length(m2$coefficients)))
[1] 1
> (LR <- 2 * (m2$loglik[2] - m1$loglik[2]))
[1] 4.770787
> pchisq(LR, df, lower=FALSE)
[1] 0.02894659

UPDATE #3 I changed the example following the kind answer from DWin as follows. This way the degrees of freedom should be calculated properly:

library(Hmisc)
library(rms)
library(ggplot2)
library(gridExtra)

data(pbc)
d <- pbc
rm(pbc, pbcseq)
d$status <- ifelse(d$status != 0, 1, 0)

### log likelihood test using a coxph model
m1 <- coxph(Surv(time, status) ~  albumin, data=d)
m2 <- coxph(Surv(time, status) ~  albumin + age, data=d)
# loglik = a vector of length 2 containing the log-likelihood with the initial values and with the final values of the coefficients.
m1$loglik[2]
m2$loglik[2]
(df <- abs(sum(anova(m1)$Df, na.rm=TRUE) - sum(anova(m2)$Df, na.rm=TRUE)))
(LR <- 2 * (m2$loglik[2] - m1$loglik[2])) # the most parsimonious models have to be first
pchisq(LR, df, lower=FALSE)
anova(m2, m1)

### log likelihood test using a cph model
dd = datadist(d)
options(datadist='dd')
m3 <- cph(Surv(time, status) ~  albumin, data=d)
m4 <- cph(Surv(time, status) ~  albumin + age, data=d)
# loglik = a vector of length 2 containing the log-likelihood with the initial values and with the final values of the coefficients.
m3$loglik[2]
m4$loglik[2]
(df <- abs(print(anova(m3)[, "d.f."])[['TOTAL']] - print(anova(m4)[, "d.f."])[['TOTAL']]))
(LR <- 2 * (m4$loglik[2] - m3$loglik[2])) # the most parsimonious models have to be first
pchisq(LR, df, lower=FALSE)
Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ The answer to the first question is "because you specified fun=exp in Predict". $\endgroup$
    – Glen_b
    Commented Jul 25, 2015 at 23:53

2 Answers 2

Reset to default
6
$\begingroup$

I don't have time to answer all of your questions. To answer two of them: the $P$-value is for the test of all coefficients pertaining to a given predictor. In other words it is the global influence of all terms whether linear or nonlinear or interactions involving them. This is a test of association, where the null hypothesis is no association (flatness). This particular calculation uses the Wald $\chi^2$ test, which is a generalization of the $z$-test.

You'll get better output with ggplot(Predict(...)) instead of plot(). For either you can combine the results from Predict() using the rbind function. Type ?rbind.Predict for help.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks for your answer! Regarding the P- value, does the test have a specific name? Something like student’s t-test? Do you know one or more examples in medical literature where this P values was stated? $\endgroup$
    – Gurkenhals
    Commented Jul 26, 2015 at 14:51
  • $\begingroup$ I expanded my answer to clarify. $\endgroup$
    – Frank Harrell
    Commented Jul 26, 2015 at 16:19
  • $\begingroup$ As far as I understand the P-value of the Wald χ2 test to test the null hypothesis of no association, the basic message of the P value is similar to the "P value of a HR" (do not know how to name it correctly) as shown in the output following summary(m1) of UPDATE #1. In both cases I would interpret a P value < 0.05 that the predictor is associated with the outcome ... is that correct? $\endgroup$
    – Gurkenhals
    Commented Jul 27, 2015 at 14:15
1
$\begingroup$

The "p-value" could be called an "unlikelihood statistics". It is supposed to be the extent to which the data disagrees with the Null or "alternative model". TYu are actually performing a "likelihood ratio test" when you compare two (nested) models. Although you have chosen to use the pspline function which defaults to a fairly high number of knots, it probably doesn't matter much in this case. In a real data analysis it might make sense to limit the number of knots. I find that more than 4 knots to a spline term seems to add very little in additional meaning. Both the difference in deviance or ( -2*(LL1-LL0) ) contribute to a t-statistic and the associated p-value. The more knots, then higher the hurdle over which the likelihood ratio test must jump to achieve conventional levels of "significance".

I'm having a bit of difficulty with the anova on the pspline model versus the non-spline model. If you only have a difference of (-975.61 - -973.26) I don't know why the chi-square statistic is 4.6983. Seems as though it should be a bit over 2. At any rate the df=11 is going to prevent you from concluding tat the spline model is significantly better since the critical value of chisq for a difference of 10 df is around 20.

Answer to Q in comment.: The df's of the anova-object from the coxph models is

anova(m1,m2)$Df

The df's of the f2 object can be seen with:

a2[, "d.f."]

(To examine the R objects you can either print them or use str )

Improve this answer
$\endgroup$
9
  • $\begingroup$ Oh, my example in UPDATE #1 is wrong as you pointed out that I compare not nested models with the likelihood ratio test. Howevery, In UPDATE #2 I do better. However again, since the Df in UPDATE #1 is 11, the Df seems not to be calculated as the differences of predictors as I do in UPDATE #2. Leeding me to the question: How do I get the degrees of freedom for the likelihood ratio test? $\endgroup$
    – Gurkenhals
    Commented Jul 27, 2015 at 3:41
  • $\begingroup$ The df for the LRT is the difference in df for the two models. My criticism was using a spline function that allowed too many knots. I also use Harrell's rms package and when I do so, I generally set the max number of knots to 4 using the rcs function. 2 or 3 knots usually suffices. Interpreting the results I get with 4 knots is not particularly useful. $\endgroup$
    – DWin
    Commented Jul 27, 2015 at 4:00
  • $\begingroup$ Okay thanks. I understand that. However, do you know how I can calculate the difference of degree of freedoms between two models? $\endgroup$
    – Gurkenhals
    Commented Jul 27, 2015 at 4:28
  • $\begingroup$ It's not clear what you don't understand. Take the df of model1 and subtract the df of model2. That the df of for the chi-square statistic. $\endgroup$
    – DWin
    Commented Jul 27, 2015 at 5:37
  • 1
    $\begingroup$ Please spend a bit more time studying the software. The d.f. is clearly printed by the L.R. chi-square statistic when you print the model fit object. d.f. are also printed when you use lrtest(fit1 fit2). And remember that the d.f. for the test of linearity against a spline function is equal to the number of knots minus 2 when using a restricted cubic spline function. $\endgroup$
    – Frank Harrell
    Commented Jul 27, 2015 at 15:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.