source("http://192.38.117.59/~linearpredictors/datafiles/readPbc3.R")
library(survival)

## Anova gives likelihood ratio tests for linearity
## Note that the P-values in the table are Wald tests
## These can be obtained (for the 7th and 8th coefficents simultaneously) as:

WaldLinearity <- function(model){
  beta <- model$coef[7:8]
  V <- model$var[7:8, 7:8]
  W <- beta %*% solve(V) %*% beta
  list(W = W,
       P = 1-pchisq(W, df = 2)
       )
}

(model1 <- coxph(Surv(followup, status != 0) ~ tment + alb + log(bili) + log(alkph) +
                log(asptr) + age, data = pbc3))

(model2 <- update(model1, . ~ . + ifelse(alb > 30, alb - 30, 0) + ifelse(alb > 40, alb - 40,0)))

anova(model1,model2)
WaldLinearity(model2)


(model3 <- update(model1, . ~ . + ifelse(log(bili) > log(10), log(bili) - log(10), 0) + ifelse(log(bili) > log(40), log(bili) - log(40),0)))

anova(model1,model3)
WaldLinearity(model3)

(model4 <- update(model1, . ~ . + ifelse(log(alkph) > log(500), log(alkph) - log(500), 0) + ifelse(log(alkph) > log(1200), log(alkph) - log(1200),0)))

anova(model1,model4)
WaldLinearity(model4)

(model5 <- update(model1, . ~ . + ifelse(log(asptr) > log(60), log(asptr) - log(60), 0) + ifelse(log(asptr) > log(120), log(asptr) - log(120),0)))

anova(model1,model5)
WaldLinearity(model5)

(model6 <- update(model1, . ~ . + ifelse(age > 50, age - 50, 0) + ifelse(age > 60, age - 60,0)))

anova(model1,model6)
WaldLinearity(model6)