Case study II: Leukemia data

This document presents an analysis of Near Redundancy and Practical non-identifiability using the Leukemia data set from the R package spBayesSurv. For more information see Rubio et al. (2019) and Espindola, Montoya, and Rubio (2023).

See the GitHub repository: Near Redundancy and Practical non-identifiability: Survival Models

Data preparation

########################################################################
# Case study II: Leukemia data set from the spBayesSurv R package
########################################################################

# See README file for more information
rm(list = ls())

# Routines
# https://github.com/FJRubio67/NRPNISurv
source("routines.R")

# Required packages (from CRAN)
# install.packages('PackageName')
library(knitr)
library(survival)
library(spBayesSurv)
library(numDeriv)

# Required packages (from GitHub)
#library(devtools)
#install_github("FJRubio67/HazReg")
library(HazReg)


#----------------------------------------------------------------------------------------------
# Leukemia data
#---------------------------------------------------------------------------------------------

data(LeukSurv)
?LeukSurv
head(LeukSurv)

##   time cens    xcoord    ycoord age sex   wbc   tpi district
## 1    1    1 0.2050717 0.4972437  61   0  13.3 -1.96        9
## 2    1    1 0.2855568 0.8489526  76   0 450.0 -3.39        7
## 3    1    1 0.1764057 0.7364939  74   0 154.0 -4.95        7
## 4    1    1 0.2447630 0.2105843  79   1 500.0 -1.40       24
## 5    1    1 0.3274531 0.9073870  83   1 160.0 -2.59        7
## 6    1    1 0.6383682 0.3627343  81   1  30.4  0.03       11

dim(LeukSurv)

## [1] 1043    9

#******************
# Data preparation
#******************

n <- dim(LeukSurv)[1]  # number of individuals
# Design matrices
x <- as.matrix(cbind(scale(LeukSurv$age), LeukSurv$sex, scale(LeukSurv$wbc), scale(LeukSurv$tpi) ))
colnames(x) <- cbind("std age", "sex", "wbc", "tpi")
xt <- as.matrix(cbind(scale(LeukSurv$age), scale(LeukSurv$wbc), scale(LeukSurv$tpi)))

# Required quantities
status <- as.vector(LeukSurv$cens)
times <- as.vector(LeukSurv$time)/365.24 # in years

Maximum likelihood estimation

#----------------------------------------------------------------------------------------------
# Maximum likelihood estimation
# See : https://github.com/FJRubio67/HazReg
#----------------------------------------------------------------------------------------------

# PGW model with no covariates
OPTPGW0 <- GHMLE(init = c(0,0,0), times = times, status = status, 
                  hstr = "baseline", dist = "PGW", method = "nlminb", maxit = 10000)

# PGW-PH model
OPTPGWPH <- GHMLE(init = c(0,0,0, rep(0,ncol(x))), times = times, status = status, 
     hstr = "PH", dist = "PGW", des = x, method = "nlminb", maxit = 10000)

# PGW-AFT model
OPTPGWAFT <- GHMLE(init = c(0,0,0,rep(0,ncol(x))), times = times, status = status, 
        hstr = "AFT", dist = "PGW", des = x, method = "nlminb", maxit = 10000)

# PGW-GH model
OPTPGWGH <- GHMLE(init = c(0,0,0, rep(0, ncol(xt)), rep(0,ncol(x))), times = times, status = status, 
            hstr = "GH", dist = "PGW", des = x, des_t = xt, method = "nlminb", maxit = 10000)

# MLE PGW no covariates
MLE.PGW0 <- exp(OPTPGW0$OPT$par[1:3])
kable(MLE.PGW0, digits = 3)

x
0.082
0.826
3.160

# MLE PGW-PH structure
MLE.PGWPH <- c(exp(OPTPGWPH$OPT$par[1:3]), OPTPGWPH$OPT$par[-c(1:3)])
kable(MLE.PGWPH, digits = 3)

x
0.139
0.815
2.570
0.536
0.064
0.217
0.097

# MLE PGW-AFT structure
MLE.PGWAFT <- c(exp(OPTPGWAFT$OPT$par[1:3]), OPTPGWAFT$OPT$par[-c(1:3)])
kable(MLE.PGWAFT, digits = 3)

x
0.093
1.014
3.504
1.028
0.081
0.484
0.215

# MLE PGW-GH structure
MLE.PGWGH <- c(exp(OPTPGWGH$OPT$par[1:3]), OPTPGWGH$OPT$par[-c(1:3)])
kable(MLE.PGWGH, digits = 3)

x
0.095
1.006
3.474
0.911
0.898
0.413
0.979
0.077
0.681
0.303

Model selection

# AIC
AIC.PGW0 <- 2*OPTPGW0$OPT$objective + 2*length(MLE.PGW0)
AIC.PGWGH <- 2*OPTPGWGH$OPT$objective + 2*length(MLE.PGWGH)
AIC.PGWAFT <- 2*OPTPGWAFT$OPT$objective + 2*length(MLE.PGWAFT)
AIC.PGWPH <- 2*OPTPGWPH$OPT$objective + 2*length(MLE.PGWPH)

# Best model: PGW-GH
c(AIC.PGW0, AIC.PGWPH, AIC.PGWAFT, AIC.PGWGH)

## [1] 1836.589 1586.483 1539.478 1537.239

Detecting Near Redundancy

######################################################################
# KL Divergence criterion
######################################################################
# KL Divergence
dKL <- minKL(MLE.PGWGH[1:3], c(1,1))$min.KL
dKL

## [1] 0.06162896

# Redundancy constant
k = length(MLE.PGWGH)
M = 0.05
ne <- n - 0.5*sum(!status)

UKL <- 0.5*k*M*length(MLE.PGWGH)*log(ne)/ne

# Test
(dKL < UKL)

## [1] FALSE

######################################################################
# Hellinger distance criterion
######################################################################

# Hellinger Distance
DH <- minHell(MLE.PGWGH[1:3], c(1,1))$min.dHell
DH

## [1] 0.1007712

# Hellinger criterion
kappa <- 0.05
UH <- log( 1 - (1-2*kappa)^2)/(2*log(1-DH^2))

# Test
(ne < UH)

## [1] FALSE

######################################################################
# Hessian Method
######################################################################

# Hessian/Eigenvalues
HESS.PGW  <- -hessian(OPTPGWGH$log_lik, x = OPTPGWGH$OPT$par)
  
eigen.val <- eigen(HESS.PGW)$values

ref.val <- abs(as.vector(eigen.val))/max(abs(as.vector(eigen.val)))

sev <- sort(ref.val)

sev

##  [1] 0.003533086 0.004102972 0.005640710 0.007719617 0.040147119 0.097920316
##  [7] 0.188310565 0.259869672 0.281508018 1.000000000

# Upper bound
Un <- 0.001

# Test
(sev[1] < Un)

## [1] FALSE

Bootstrap for distance-based criteria and Hessian method

#######################################################################################
# Bootstrap
#######################################################################################

# Number of Boostrap samples 
B <- 1000

# Bootstrap MLEs
MLE.B <- matrix(0, ncol = length(MLE.PGWGH), nrow = B)
neB <- UKLB <- UHB <- DHB<- DKLB <- vector()
indHess <- vector()

for(i in 1:B){
  ind <- sample(1:n, replace = T)
  OPTB <- GHMLE(init = c(0,0,0, rep(0, ncol(xt)), rep(0,ncol(x))), times = times[ind], status = status[ind], 
         hstr = "GH", dist = "PGW", des = x[ind,], des_t = xt[ind,], method = "nlminb", maxit = 10000)
  MLE.B[i,] <- c(exp(OPTB$OPT$par[1:3]), OPTB$OPT$par[-c(1:3)])
  neB[i] <- n - 0.5*sum(!status[ind])
  
  UKLB[i] <- 0.5*k*M*length(MLE.B[i,])*log(neB[i])/neB[i]
  
  # Hessian/Eigenvalues
  HESSB  <- -hessian(OPTB$log_lik, x = OPTB$OPT$par)
  
  eigen.val <- eigen(HESSB)$values
  
  ref.val <- abs(as.vector(eigen.val))/max(abs(as.vector(eigen.val)))
  
  sev <- sort(ref.val)
  
  # Test
  indHess[i] <- as.numeric(sev[1] < Un)
}


######################################################################
# Bootstrap KL Divergence criterion
######################################################################

# Bootstrap KL Divergence

for(i in 1:B){
  DKLB[i] <- minKL(MLE.B[i,1:3], c(1,1))$min.KL
}

# Bootstrap probability of KL divergence criterion

indKLB <- (DKLB < UKLB) 

mean(indKLB)

## [1] 0

######################################################################
# Bootstrap Hellinger distance criterion
######################################################################

# Hellinger Distance
for(i in 1:B){
DHB[i] <- minHell(MLE.B[i,1:3], c(1,1))$min.dHell

# Hellinger criterion
UHB[i] <- log( 1 - (1-2*kappa)^2)/(2*log(1-DHB[i]^2))
}

# Bootstrap probability of Hellinger distance criterion
indHB <- (neB < UHB)

mean(indHB)

## [1] 0

######################################################################
# Bootstrap Hessian method
######################################################################


# Bootstrap probability of KL divergence criterion

mean(indHess)

## [1] 0

Detecting Practical Non-Identifiability

#######################################################################################
# Profile likelihood
#######################################################################################
p0 <- ncol(xt)
p1 <- ncol(x)

p <- length(MLE.PGWGH) # number of parameters
ML <- OPTPGWGH$OPT$objective # Maximum likelihood


# Profile likelihood function for parameter "ind"
prof.lik <- function(par1, ind){
  
  tempf <- function(par){
    tempv <- rep(0,p)
    tempv <- replace(x = tempv, c(1:p)[-ind] , par)
    tempv <- replace(x = tempv, ind , par1)
    out0 <- OPTPGWGH$log_lik(tempv)
    return(out0)
  } 
  
  out0 <-  -nlminb(OPTPGWGH$OPT$par[-ind],tempf, control = list(iter.max = 10000))$objective + ML
  out1 <-  -nlminb(OPTPGWGH$OPT$par[-ind],tempf, control = list(iter.max = 10000))$objective + ML
  out <- min(out0, out1)
  
  out2 <- ifelse(exp(out)<=1, exp(out), 0)
  return(out2)
}


# Profile likelihoods

# Profile likelihood of Parameter 1
prof1 <- Vectorize(function(par) prof.lik(log(par),1))
curve(prof1,0.05,0.2 , n = 200, lwd = 2, xlab = expression(sigma), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 2
prof2 <- Vectorize(function(par) prof.lik(log(par),2))
curve(prof2,0.8,1.3 , n = 200, lwd = 2, xlab = expression(nu), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 3
prof3 <- Vectorize(function(par) prof.lik(log(par),3))
curve(prof3,2,5 , n = 200, lwd = 2, xlab = expression(gamma), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 4
prof4 <- Vectorize(function(par) prof.lik(par,4))
curve(prof4,0.2,1.6, n = 200, lwd = 2, xlab = expression(alpha[1]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 5
prof5 <- Vectorize(function(par) prof.lik(par,5))
curve(prof5,0.25, 1.75 , n = 200, lwd = 2, xlab = expression(alpha[2]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 6
prof6 <- Vectorize(function(par) prof.lik(par,6))
curve(prof6,-0.5, 1.25 , n = 200, lwd = 2, xlab = expression(alpha[3]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 7
prof7 <- Vectorize(function(par) prof.lik(par,7))
curve(prof7,0.6,1.4, n = 200, lwd = 2, xlab = expression(beta[1]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 8
prof8 <- Vectorize(function(par) prof.lik(par,8))
curve(prof8,-0.15, 0.35 , n = 200, lwd = 2, xlab = expression(beta[2]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 9
prof9 <- Vectorize(function(par) prof.lik(par,9))
curve(prof9,0.2, 1.2 , n = 200, lwd = 2, xlab = expression(beta[3]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 10
prof10 <- Vectorize(function(par) prof.lik(par,10))
curve(prof10,-0.2,0.8 , n = 200, lwd = 2, xlab = expression(beta[4]), ylab = "Profile Likelihood")

Alternative Models

##################################################################################
# Alternative models
##################################################################################


# Exponentiated Weibull

OPTEWGH <- GHMLE(init = c(0,0,0, rep(0, ncol(xt)), rep(0,ncol(x))), times = times, status = status, 
                 hstr = "GH", dist = "EW", des = x, des_t = xt, method = "nlminb", maxit = 10000)

# MLE EW-GH structure
MLE.EWGH <- c(exp(OPTEWGH$OPT$par[1:3]), OPTEWGH$OPT$par[-c(1:3)])
kable(MLE.EWGH, digits = 3)

x
0.006
0.219
8.957
0.886
0.847
0.407
0.969
0.071
0.656
0.302

AIC.EWGH <- 2*OPTEWGH$OPT$objective + 2*length(MLE.EWGH)
AIC.EWGH

## [1] 1533.727

#######################################################################################
# Profile likelihood
#######################################################################################


p2 <- length(MLE.EWGH) # number of parameters
ML2 <- OPTEWGH$OPT$objective # Maximum likelihood


# Profile likelihood function for parameter "ind"
prof.lik2 <- function(par1, ind){
  
  tempf <- function(par){
    tempv <- rep(0,p2)
    tempv <- replace(x = tempv, c(1:p2)[-ind] , par)
    tempv <- replace(x = tempv, ind , par1)
    out0 <- OPTEWGH$log_lik(tempv)
    return(out0)
  } 
  
  out <-  -nlminb(OPTEWGH$OPT$par[-ind],tempf, control = list(iter.max = 10000))$objective + ML2

  out2 <- ifelse(exp(out)<=1, exp(out), 0)
  return(out2)
}


# Profile likelihoods

# Profile likelihood of Parameter 1
prof21 <- Vectorize(function(par) prof.lik2(log(par),1))
curve(prof21,0.0001,0.075 , n = 200, lwd = 2, xlab = expression(sigma), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 2
prof22 <- Vectorize(function(par) prof.lik2(log(par),2))
curve(prof22,0.125,0.325 , n = 200, lwd = 2, xlab = expression(nu), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 3
prof23 <- Vectorize(function(par) prof.lik2(log(par),3))
curve(prof23,3.5,25 , n = 200, lwd = 2, xlab = expression(gamma), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 4
prof24 <- Vectorize(function(par) prof.lik2(par,4))
curve(prof24,0.2,1.6, n = 200, lwd = 2, xlab = expression(alpha[1]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 5
prof25 <- Vectorize(function(par) prof.lik2(par,5))
curve(prof25,0.25, 1.5 , n = 200, lwd = 2, xlab = expression(alpha[2]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 6
prof26 <- Vectorize(function(par) prof.lik2(par,6))
curve(prof26,-0.25, 1 , n = 200, lwd = 2, xlab = expression(alpha[3]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 7
prof27 <- Vectorize(function(par) prof.lik2(par,7))
curve(prof27,0.7,1.2, n = 200, lwd = 2, xlab = expression(beta[1]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 8
prof28 <- Vectorize(function(par) prof.lik2(par,8))
curve(prof28,-0.125, 0.3 , n = 200, lwd = 2, xlab = expression(beta[2]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 9
prof29 <- Vectorize(function(par) prof.lik2(par,9))
curve(prof29,0.4, 1.1 , n = 200, lwd = 2, xlab = expression(beta[3]), ylab = "Profile Likelihood")

# Profile likelihood of Parameter 10
prof210 <- Vectorize(function(par) prof.lik2(par,10))
curve(prof210,0,0.75 , n = 200, lwd = 2, xlab = expression(beta[4]), ylab = "Profile Likelihood")

Espindola, J. A., J. A. Montoya, and F. J. Rubio. 2023. “On Near Redundancy and Identifiability of Parametric Hazard Regression Models Under Censoring.” Biometrical Journal 65: 2300006.

Rubio, F. J., L. Remontet, N. P. Jewell, and A. Belot. 2019. “On a General Structure for Hazard-Based Regression Models: An Application to Population-Based Cancer Research.” Statistical Methods in Medical Research 28: 2404–17.

Near Redundancy and Practical non-identifiability: Leukemia data

F. Javier Rubio, J.A. Espindola, and J.A. Montoya

14 July, 2024