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

## Note that the likelihood tests differ slightly from those in the book,
## as the latter are calculated using the Breslow method to handle ties,
## rather than Efron which is the default method in R.

## Also, be aware that some of the models with terms of powers 2 and 3
## have problems with singularities.

## fit the basic Cox model with linear effect of bilirubin

linear <- coxph(Surv(followup,status!=0)~ bili,
                data = subset(pbc3, followup > 0))

## We are interested in including bilirubin raised to the powers of
powers <- c(-3, -2, -1, -0.5, 0, 0.5, 2, 3)
## where 0 is the log-tranformation and not a constant polynomial

## Now, fit all the two-covariate Cox models from the table:

for (q1 in powers){
  if (q1 == 0){
    ## fit model with bili and log(bili)
    twocovariates <- coxph(Surv(followup,status!=0)~ bili + log(bili),
                           data = subset(pbc3, followup > 0))
  }
  else{
    ## fit model with bili and bili^q1
    twocovariates <- coxph(Surv(followup,status!=0)~ bili + I(bili^q1),
                           data = subset(pbc3, followup > 0)
                           )
  }  
  
  ## finally calculate the differences in -2logL
  difference <- 2*summary(twocovariates)$loglik[2] - 2*summary(linear)$loglik[2]

  ## and print this to screen in a readable manner
  cat("q1 = ",q1,": \t", difference,"\n", sep="")
}

## And now to the three-covariate Cox-models
for (q1 in powers){
  for (q2 in powers[powers < q1]){

    if (q1 != 0){
      q1term <- pbc3$bili^q1
    }
    else q1term <- log(pbc3$bili)
    
    if (q2 != 0){
      q2term <- pbc3$bili^q2
    }
    else q2term <- log(pbc3$bili)

    rm(threecovariates, difference)
    threecovariates <- coxph(Surv(followup,status!=0)~ bili + q1term + q2term,
                           data = subset(pbc3, followup >0))
    
    ## differences in -2logL
    difference <- 2*summary(threecovariates)$loglik[2] - 2*summary(linear)$loglik[2]
    
    ## print this to screen in a readable manner
    cat("q1 = ",q1, ",\t q2 = ", q2,"\t", round(difference, digits =3),"\n") 
  }
}