Harmonic oscillator hazard function

This document provides codes and information about the Damped Harmonic Oscillator model proposed in Christen and Rubio (2024).

The damped harmonic oscillator model

The harmonic oscillator hazard model is defined through the linear second order ODE: \[\begin{equation} h''(t) + 2 \eta w_0 h'(t) + w^2_0 (h(t) - h_b) = 0; \quad h(0) = h_0, \quad h'(0) = r_0, \label{eq:harmonic_osc} \end{equation}\] where \(w_0 > 0\) is called the “natural frequency”, \(\eta >0\) is the “damping coefficient”. The shift parameter, \(h_b > 0\), represents the stability state or equilibrium level of the solution. The resulting model has closed-form hazard and cumulative hazard functions, facilitating likelihood and Bayesian inference on the parameters, as detailed in Christen and Rubio (2024).

The following R code presents an example from Christen and Rubio (2024) where we fit the damped harmonic oscillator model to the rotterdam data set from the R package survival. We present both Likelihood-based and Bayesian analyses.

Data preparation

#######################################################################################
# Data preparation
#######################################################################################

rm(list=ls())

# Required packages
library(deSolve)
library(survival)
library(ggplot2)
#library(devtools)
#install_github("FJRubio67/HazReg")
library(HazReg)
library(spBayes)
library(Rtwalk)
library(knitr)
library(bshazard)

source("routinesHO.R")

#--------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------
# Data preparation
#--------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------


head(rotterdam)

##      pid year age meno  size grade nodes pgr  er hormon chemo rtime recur dtime
## 1393   1 1992  74    1  <=20     3     0  35 291      0     0  1799     0  1799
## 1416   2 1984  79    1 20-50     3     0  36 611      0     0  2828     0  2828
## 2962   3 1983  44    0  <=20     2     0 138   0      0     0  6012     0  6012
## 1455   4 1985  70    1 20-50     3     0   0  12      0     0  2624     0  2624
## 977    5 1983  75    1  <=20     3     0 260 409      0     0  4915     0  4915
## 617    6 1983  52    0  <=20     3     0 139 303      0     0  5888     0  5888
##      death
## 1393     0
## 1416     0
## 2962     0
## 1455     0
## 977      0
## 617      0

dim(rotterdam)

## [1] 2982   15

# Kaplan-Meier estimator for the survival times
km <- survfit(Surv(rotterdam$dtime/365.24, rotterdam$death) ~ 1)

plot(km$time, km$surv, type = "l", col = "black", lwd = 2, lty = 1, 
     ylim = c(0,1), xlab = "Time", ylab = "Survival")

# New data frame: logical status, time in years, survival times sorted
df <- data.frame(time = rotterdam$dtime, status = rotterdam$death)
df$status <- as.logical(rotterdam$death)
df$time <- df$time/365.24

df <- df[order(df$time),]

# Required quantities
status <- as.logical(df$status)
t_obs <- df$time[status]
survtimes <- df$time

Maximum Likelihood Analysis

Weibull Model

#==================================================================================================
# Maximum Likelihood Analysis
#==================================================================================================


#--------------------------------------------------------------------------------------------------
# Fitting a Weibull distribution
#--------------------------------------------------------------------------------------------------

# Initial value
initW <- c(0,0)

# Optimisation step
OPTW <- GHMLE(initW, survtimes, status, hstr = "baseline", dist = "Weibull", method = "nlminb", maxit = 10000)

MLEW <- exp(OPTW$OPT$par)

# Fitted Weibull hazard
fithw <- Vectorize(function(t) hweibull( t, exp(OPTW$OPT$par[1]), exp(OPTW$OPT$par[2]) ) )

# Fitted Weibull cumulative hazard
fitchw <- Vectorize(function(t) chweibull( t, exp(OPTW$OPT$par[1]), exp(OPTW$OPT$par[2]) ) )

# Fitted Weibull survival
fitsw <- Vectorize(function(t) exp(-chweibull( t, exp(OPTW$OPT$par[1]), exp(OPTW$OPT$par[2]) ) ))

# AIC
AICW <- 2*OPTW$OPT$objective + 2*length(OPTW$OPT$par)

# BIC
BICW <- 2*OPTW$OPT$objective + length(OPTW$OPT$par)*log(length(survtimes))

The Power Generalised Weibull Model

#==================================================================================================
# Maximum Likelihood Analysis
#==================================================================================================

#--------------------------------------------------------------------------------------------------
# Fitting a PGW distribution
#--------------------------------------------------------------------------------------------------

# Initial value
initPGW <- c(0,0,0)

# Optimisation step
OPTPGW <- GHMLE(initPGW, survtimes, status, hstr = "baseline", dist = "PGW", 
                method = "nlminb", maxit = 10000)

MLEPGW <- exp(OPTPGW$OPT$par)

# Fitted Weibull hazard
fithpgw <- Vectorize(function(t) hpgw( t, MLEPGW[1], MLEPGW[2], MLEPGW[3] ) )

# Fitted Weibull cumulative hazard
fitchpgw <- Vectorize(function(t) chpgw( t, MLEPGW[1], MLEPGW[2], MLEPGW[3] ) )

# Fitted Weibull survival
fitspgw <- Vectorize(function(t) exp(-chpgw( t, MLEPGW[1], MLEPGW[2], MLEPGW[3] ) ))

# AIC
AICPGW <- 2*OPTPGW$OPT$objective + 2*length(OPTPGW$OPT$par)

# BIC
BICPGW <- 2*OPTPGW$OPT$objective + length(OPTPGW$OPT$par)*log(length(survtimes))

The Damped Harmonic Oscillator Model

#==================================================================================================
# MLE Analysis
#==================================================================================================

#--------------------------------------------------------------------------------------------------
# Harmonic Oscillator model for the hazard function: Solver solution
#--------------------------------------------------------------------------------------------------

# Initial conditions

# Survival at time 0
S0 <- 1
# Survival at 1 month
St <- 0.999
# Survival at 2 months
Stt <- 0.998
# dt = 1 month
dt <- 1/12

# Approximation of the initial conditions
Sp <- -(St-S0)/dt
Spp <- (Stt - 2*St + S0)/(dt^2)

h00 <- Sp/St
r00 <- h00^2 - Spp/St

# Initial point
initHO <- c(0.1,0,0)

# Optimisation step
OPTHO <- optim(initHO, log_likHO, control = list(maxit = 1000))
OPTHO

## $par
## [1] -0.8421869  0.2711340 -2.7299280
## 
## $value
## [1] 4778.518
## 
## $counts
## function gradient 
##      194       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

MLEHO <- exp(OPTHO$par)

# AIC
AICHO <- 2*OPTHO$value + 2*length(OPTHO$par)

# BIC
BICHO <- 2*OPTHO$value + length(OPTHO$par)*log(length(survtimes))

Comparison

# AIC comparison
AICs <- c(AICW, AICPGW, AICHO)

AICs

## [1] 9638.298 9572.033 9563.036

# Best model (Harmonic oscillator)
which.min(AICs)

## [1] 3

# BIC comparison
BICs <- c(BICW, BICPGW, BICHO)

BICs

## [1] 9650.299 9590.034 9581.037

# Best model (Harmonic oscillator)
which.min(BICs)

## [1] 3

# Fitted hazard, cumulative hazard, and survival
fithHO <- Vectorize(function(t) hHO(t,exp(OPTHO$par[1]), exp(OPTHO$par[2]), exp(OPTHO$par[3]), h00, r00))
fitchHO <- Vectorize(function(t) chHO(t,exp(OPTHO$par[1]), exp(OPTHO$par[2]), exp(OPTHO$par[3]), h00, r00))
fitsHO <- Vectorize(function(t) exp(-chHO(t,exp(OPTHO$par[1]), exp(OPTHO$par[2]), exp(OPTHO$par[3]), h00, r00)))

# Comparison: hazard functions
curve(fithHO, 0, max(survtimes), xlim= c(0,max(survtimes)), ylim = c(0,0.125), lwd = 2,
     xlab = "Time", ylab = "Hazard", main = "Harmonic Oscillator", cex.axis = 1.5, cex.lab = 1.5)
curve(fithw, 0, max(survtimes), lwd= 2, lty = 2, col = "gray", add = TRUE)
curve(fithpgw, 0, max(survtimes), lwd= 2, lty = 3, col = "gray", add = TRUE)
legend("topleft", legend = c("Harmonic Oscillator","Weibull", "PGW"), lty = c(1,2,3), 
       lwd = c(2,2,2), col = c("black","gray","gray"))

# Comparison: cumulative hazard functions
curve(fitchHO, 0, max(survtimes), xlim = c(0,max(survtimes)), ylim = c(0,2.25),  lwd = 2,
     xlab = "Time", ylab = "Cumulative Hazard", main = "Harmonic Oscillator", cex.axis = 1.5, cex.lab = 1.5)
curve(fitchw, 0, max(survtimes), lwd= 2, lty = 2, col = "gray", add = TRUE)
curve(fitchpgw, 0, max(survtimes), lwd= 2, lty = 3, col = "gray", add = TRUE)
legend("topleft", legend = c("Harmonic Oscillator","Weibull","PGW"), lty = c(1,2,3), 
       lwd = c(2,2,2), col = c("black","gray","gray"))

# Comparison: survival functions
curve(fitsHO, 0, max(survtimes), xlim = c(0,max(survtimes)), ylim = c(0,1), lwd = 2,
     xlab = "Time", ylab = "Survival", main = "Harmonic Oscillator", cex.axis = 1.5, cex.lab = 1.5)
curve(fitsw, 0, max(survtimes), lwd= 2, lty = 2, col = "gray", add = TRUE)
curve(fitspgw, 0, max(survtimes), lwd= 2, lty = 3, col = "gray", add = TRUE)
points(km$time, km$surv, type = "l", col = "green", lwd = 2, lty = 1)
legend("topright", legend = c("Harmonic Oscillator","Weibull","PGW","KM"), lty = c(1,2,3,1), 
       lwd = c(2,2,2,2), col = c("black","gray","gray","green"))

Bayesian Analysis

Weibull

#==================================================================================================
# Bayesian Analysis
#==================================================================================================

#--------------------------------------------------------------------------------------------------
# Priors
#--------------------------------------------------------------------------------------------------

par(mfrow = c(1,2))
p_sigmaW <- Vectorize(function(t) dgamma(t, shape = 0.001, scale = 1/0.001))
curve(p_sigmaW,0,10, n = 1000, xlab = ~sigma, ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2)

p_nuW <- Vectorize(function(t) dgamma(t, shape = 0.001, scale = 1/0.001))
curve(p_nuW,0,10, n = 1000, xlab = ~nu, ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2)

par(mfrow = c(1,1))

#--------------------------------------------------------------------------------------------------
# Weibull model for the hazard function
#--------------------------------------------------------------------------------------------------

# Support
SupportW <- function(x) {   TRUE }

# Random initial points
X0W <- function(x) { OPTW$OPT$par + runif(2,-0.01,0.01) }

# twalk for analytic solution
set.seed(1234)
infoW <- Runtwalk( dim=2,  Tr=110000,  Obj=log_postW, Supp=SupportW, 
                   x0=X0W(), xp0=X0W(), PlotLogPost = FALSE)

## This is the twalk for R.
## Evaluating the objective density at initial values.
## Opening 6 X 110001 matrix to save output.  Sampling (no graphics mode).

# Posterior sample after burn-in and thinning
ind=seq(10000,110000,100) 

# Summaries
summW <- apply(exp(infoW$output[ind,]),2,summary)
colnames(summW) <- c("sigma","nu")
kable(summW, digits = 3)

	sigma	nu
Min.	13.435	1.164
1st Qu.	14.290	1.233
Median	14.555	1.253
Mean	14.555	1.254
3rd Qu.	14.804	1.275
Max.	15.738	1.363

# KDEs
sigmapW <- exp(infoW$output[,1][ind])
nupW <- exp(infoW$output[,2][ind])

plot(density(sigmapW), main = "", xlab = expression(sigma), ylab = "Density",
     cex.axis = 1.5, cex.lab = 1.5, lwd = 2)
curve(p_sigmaW,13,17, n = 1000, xlab = ~lambda, ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2, add = TRUE)

plot(density(nupW), main = "", xlab = expression(nu), ylab = "Density",
     cex.axis = 1.5, cex.lab = 1.5, lwd = 2)
curve(p_nuW,1.1,1.5, n = 1000, xlab = ~kappa, ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2, add = TRUE)

The Power Generalised Weibull Model

#==================================================================================================
# Bayesian Analysis
#==================================================================================================

#--------------------------------------------------------------------------------------------------
# Priors
#--------------------------------------------------------------------------------------------------

par(mfrow = c(1,3))
p_sigmaPGW <- Vectorize(function(t) dgamma(t, shape = 0.001, scale = 1/0.001))
curve(p_sigmaPGW,0,10, n = 1000, xlab = ~sigma, ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2)

p_nuPGW <- Vectorize(function(t) dgamma(t, shape = 0.001, scale = 1/0.001))
curve(p_nuPGW,0,10, n = 1000, xlab = ~nu, ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2)

p_gammaPGW <- Vectorize(function(t) dgamma(t, shape = 0.001, scale = 1/0.001))
curve(p_gammaPGW,0,10, n = 1000, xlab = ~gamma, ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2)

par(mfrow = c(1,1))

#--------------------------------------------------------------------------------------------------
# Weibull model for the hazard function
#--------------------------------------------------------------------------------------------------

# Support
SupportPGW <- function(x) {   TRUE }

# Random initial points
X0PGW <- function(x) { OPTPGW$OPT$par + runif(3,-0.01,0.01) }


# twalk for analytic solution
set.seed(1234)
infoPGW <- Runtwalk( dim=3,  Tr=110000,  Obj=log_postPGW, Supp=SupportPGW, 
                     x0=X0PGW(), xp0=X0PGW(), PlotLogPost = FALSE)

## This is the twalk for R.
## Evaluating the objective density at initial values.
## Opening 8 X 110001 matrix to save output.  Sampling (no graphics mode).

# Posterior sample after burn-in and thinning
ind=seq(10000,110000,100) 

# Summaries
summPGW <- apply(exp(infoPGW$output[ind,]),2,summary)
colnames(summPGW) <- c("sigma","nu","gamma")
kable(summPGW, digits = 3)

	sigma	nu	gamma
Min.	2.187	1.653	3.531
1st Qu.	2.637	1.939	5.063
Median	2.808	2.018	5.469
Mean	2.832	2.025	5.509
3rd Qu.	2.990	2.110	5.982
Max.	4.156	2.506	7.697

# KDEs
sigmapPGW <- exp(infoPGW$output[,1][ind])
nupPGW <- exp(infoPGW$output[,2][ind])
gammapPGW <- exp(infoPGW$output[,3][ind])

plot(density(sigmapPGW), main = "", xlab = expression(sigma), ylab = "Density",
     cex.axis = 1.5, cex.lab = 1.5, lwd = 2)
curve(p_sigmaPGW,2,5, n = 1000, xlab = ~lambda, ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2, add = TRUE)

plot(density(nupPGW), main = "", xlab = expression(nu), ylab = "Density",
     cex.axis = 1.5, cex.lab = 1.5, lwd = 2)
curve(p_nuPGW,1.5,2.7, n = 1000, xlab = ~kappa, ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2, add = TRUE)

plot(density(gammapPGW), main = "", xlab = expression(gamma), ylab = "Density",
     cex.axis = 1.5, cex.lab = 1.5, lwd = 2)
curve(p_gammaPGW,0,15, cex.lab = 1.5, lwd = 2, lty = 2, add = TRUE)

The Harmonic Oscillator hazard Model

#==================================================================================================
# Bayesian Analysis
#==================================================================================================

#--------------------------------------------------------------------------------------------------
# Priors
#--------------------------------------------------------------------------------------------------

par(mfrow = c(2,2))
p_eta <- Vectorize(function(t) dgamma(t, shape = 0.001, scale = 1/0.001))
curve(p_eta,0,10, n = 1000, xlab = ~eta, ylab = "Prior Density",
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2)

p_w0 <- Vectorize(function(t) dgamma(t, shape = 0.001, scale = 1/0.001))
curve(p_w0,0,10, n = 1000, xlab = expression(w_0), ylab = "Prior Density",
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2)

p_hb <- Vectorize(function(t) dgamma(t, shape = 0.001, scale = 1/0.001))
curve(p_hb,0,10, n = 1000, xlab = expression(h_b), ylab = "Prior Density",
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2)


par(mfrow = c(1,1))

#--------------------------------------------------------------------------------------------------
# Posterior sampling
#--------------------------------------------------------------------------------------------------


n.batch <- 1100
batch.length <- 100

lp <- function(par) -log_postHO(par)

inits <- OPTHO$par

set.seed(1234)
infoHO <- adaptMetropGibbs(ltd=lp, starting=inits, accept.rate=0.44, batch=n.batch, 
                           batch.length=batch.length, report=100, verbose=FALSE)

chainHO <- infoHO$p.theta.samples[,1:3]

# Burning and thinning the chain
burn <- 1e4
thin <- 100
NS <- n.batch*batch.length
ind <- seq(burn,NS,thin)



# Summaries
summHO <- apply(exp(chainHO[ind,1:3]),2,summary)
colnames(summHO) <- c("eta","w_0","h_b")
kable(summHO, digits = 3)

	eta	w_0	h_b
Min.	0.253	0.977	0.059
1st Qu.	0.389	1.251	0.064
Median	0.436	1.327	0.065
Mean	0.447	1.333	0.065
3rd Qu.	0.494	1.408	0.067
Max.	0.746	1.809	0.074

# KDEs
etap <- exp(chainHO[,1][ind])
w0p <- exp(chainHO[,2][ind])
hbp <- exp(chainHO[,3][ind])


plot(density(etap), main = "", xlab = expression(eta), ylab = "Density",
     cex.axis = 1.5, cex.lab = 1.5, lwd = 2)
curve(p_eta,0,1, n = 1000, xlab = ~lambda, ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2, add = TRUE)

plot(density(w0p), main = "", xlab = expression(w_0), ylab = "Density",
     cex.axis = 1.5, cex.lab = 1.5, lwd = 2)
curve(p_w0,0.5,2.5, n = 1000, xlab = expression(w_0), ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2, add = TRUE)

plot(density(hbp), main = "", xlab = expression(h_b), ylab = "Density",
     cex.axis = 1.5, cex.lab = 1.5, lwd = 2)
curve(p_hb,0,0.1, n = 1000, xlab = expression(h_b), ylab = "Prior Density", 
      cex.axis = 1.5, cex.lab = 1.5, lwd = 2, lty = 2, add = TRUE)

Predictive Hazard functions

#---------------------------------------------------------------------------------------
# Predictive Hazard functions
#---------------------------------------------------------------------------------------

# Predictive Weibull hazard
predhW <- Vectorize(function(t){
  num <- den <- temp <- vector()
  for(i in 1:length(ind)){
    num[i] <- exp(-chweibull( t, sigmapW[i], nupW[i]))*hweibull( t, sigmapW[i], nupW[i])
    den[i] <- exp(-chweibull( t, sigmapW[i], nupW[i]))
  }
  return(mean(num)/mean(den))
})

# Predictive PGW
predhPGW <- Vectorize(function(t){
  num <- den <- temp <- vector()
  for(i in 1:length(ind)) num[i] <- exp(-chpgw( t, sigmapPGW[i], nupPGW[i], gammapPGW[i]))*
      hpgw( t, sigmapPGW[i], nupPGW[i], gammapPGW[i])
  for(i in 1:length(ind)) den[i] <- exp(-chpgw( t, sigmapPGW[i], nupPGW[i], gammapPGW[i]))
  return(mean(num)/mean(den))
})



# Creating the credible envelopes
tvec <- seq(0,20,by = 0.01)
ntvec <- length(tvec)

# Weibull
hCIW <- matrix(0, ncol = ntvec, nrow = length(ind))

for(j in 1:length(ind)){
  for(k in 1:ntvec){
    hCIW[j,k ] <- hweibull( tvec[k],sigmapW[j], nupW[j]) 
  }
} 

hW <-  predhW(tvec)


hCIWL <- apply(hCIW, 2, ql)
hCIWU <- apply(hCIW, 2, qu)

# PGW
hCIPGW <- matrix(0, ncol = ntvec, nrow = length(ind))

for(j in 1:length(ind)){
  for(k in 1:ntvec){
    hCIPGW[j,k ] <- hpgw( tvec[k],sigmapPGW[j], nupPGW[j], gammapPGW[j]) 
  }
} 

hPGWs <-  predhPGW(tvec)

hCIPGWL <- apply(hCIPGW, 2, ql)
hCIPGWU <- apply(hCIPGW, 2, qu)

# Harmonic Oscillator

# Predictive HO hazard
predhHO <- Vectorize(function(t){
  num <- den <- temp <- vector()
  for(i in 1:length(ind)){
    chi <- chHO( t, etap[i], w0p[i], hbp[i], h00, r00 )
    num[i] <- exp(-chi)*hHO( t, etap[i], w0p[i], hbp[i], h00, r00 )
    den[i] <- exp(-chi)
  }
  return(mean(num)/mean(den))
})


# Creating the credible envelopes
tvec <- seq(0,20,by = 0.01)
ntvec <- length(tvec)

hHOs <-  predhHO(tvec)

hCIHO<- matrix(0, ncol = ntvec, nrow = length(ind))

for(j in 1:length(ind)){
  for(k in 1:ntvec){
    hCIHO[j,k ] <- hHO( tvec[k],etap[j], w0p[j], hbp[j], h00, r00) 
  }
} 

hCIHOL <- apply(hCIHO, 2, ql)
hCIHOU <- apply(hCIHO, 2, qu)

# Splines model
fit <- bshazard(Surv(df$time, df$status) ~ 1, data = df, nbin = 100, degree = 3, verbose = FALSE)

## NOTE: entry.status has been set to 0 for all.

# Predictive hazard plot
curve(predhHO, 0, 20, n = 1000, xlab = "Time", ylab = "Predictive Hazard",
      main = "Harmonic Oscillator",
      cex.axis = 1.5, cex.lab = 1.5, lwd =2, lty = 1, ylim = c(0,0.09))

points(tvec,  hCIHOL, col = "gray", type = "l")
points(tvec,  hCIHOU, col = "gray", type = "l")
polygon(c(tvec, rev(tvec)), c(hCIHOL[order(tvec)], rev(hCIHOU[order(tvec)])),
        col = "gray", border = NA)
curve(predhHO, 0, 20, n = 1000, xlab = "Time", ylab = "Predictive Hazard", 
      cex.axis = 1.5, cex.lab = 1.5, lwd =2, lty = 1, ylim = c(0,0.08), add = TRUE)

# Comparison with Weibull and PGW
#pdf("predhaz.pdf", height = 6, width = 8)
curve(predhHO, 0, 20, n = 1000, xlab = "Time", ylab = "Predictive Hazard", 
      cex.axis = 1.5, cex.lab = 1.5, lwd =2, lty = 1, ylim = c(0,0.11))
curve(fithw, 0, max(survtimes), lwd= 2, lty = 2, col = "gray", add = TRUE)
curve(fithpgw, 0, max(survtimes), lwd= 2, lty = 3, col = "gray", add = TRUE)
legend("topleft", legend = c("HO","Weibull", "PGW"), lty = c(1,2,3), 
       lwd = c(2,2,2), col = c("black","gray","gray"))

#dev.off()

Predictive Survival functions

#--------------------------------------------------------------------------------------- 
# Predictive Survival functions 
#--------------------------------------------------------------------------------------- 

# Predictive Weibull survival 
predsW <- Vectorize(function(t){ 
  temp <- vector() 
  for(i in 1:length(ind)) temp[i] <- exp(-chweibull(t,sigmapW[i], nupW[i]) )  
  return(mean(temp)) 
}) 

# Predictive PGW survival 
predsPGW <- Vectorize(function(t){ 
  temp <- vector() 
  for(i in 1:length(ind)) temp[i] <- exp(-chpgw( t,sigmapPGW[i], nupPGW[i], gammapPGW[i]) )  
  return(mean(temp)) 
}) 


# Predictive HO Survival
predsHO <- Vectorize(function(t){
  temp <- vector()
  for(i in 1:length(ind)){
    temp[i] <- exp(-chHO( t, etap[i], w0p[i], hbp[i], h00, r00 ))
  }
  return(mean(temp)) 
})

# Comparison: survival functions
#pdf("predsurv.pdf", height = 6, width = 8)
curve(predsHO, 0, 20, n = 1000, xlab = "Time", ylab = "Predictive Survival", 
      cex.axis = 1.5, cex.lab = 1.5, lwd =2, lty = 1, ylim = c(0,1))
curve(predsW, 0, max(survtimes), lwd= 2, lty = 2, col = "gray", add = TRUE)
curve(predsPGW, 0, max(survtimes), lwd= 2, lty = 3, col = "gray", add = TRUE)
points(km$time, km$surv, type = "l", col = "gray", lwd = 2, lty = 1)
curve(predsHO, 0, 20, n = 1000, xlab = "Time", ylab = "Predictive Survival", 
      cex.axis = 1.5, cex.lab = 1.5, lwd =2, lty = 1, ylim = c(0,1),add=TRUE)
legend("topright", legend = c("HO","Weibull","PGW","KM"), lty = c(1,2,3,1), 
       lwd = c(2,2,2,2), col = c("black","gray","gray","gray"))

#dev.off()

Comparison with a spline-based estimator

Now, we present a comparison of the estimated hazard function with that obtained with the R command bshazard from the R package bshazard. This estimator is based on B-Splines. See bshazard: Nonparametric Smoothing of the Hazard Function.

# fitting the estimator
fit <- bshazard(Surv(df$time, df$status) ~ 1, data = df, nbin = 100, degree = 3, verbose = FALSE)

## NOTE: entry.status has been set to 0 for all.

# Hazard 
plot(tvec,  hHOs, type = "l", xlab = "Time", ylab = "Predictive Hazard", 
       cex.axis = 1.5, cex.lab = 1.5, lwd =2, lty = 1, ylim = c(0,0.125))
points(tvec,  hCIHOL, col = "gray", type = "l")
points(tvec,  hCIHOU, col = "gray", type = "l")
polygon(c(tvec, rev(tvec)), c(hCIHOL[order(tvec)], rev(hCIHOU[order(tvec)])),
        col = "gray", border = NA)
points(fit$time, fit$hazard, type='l', lwd = 2, lty = 2, ylim = c(0,0.25), col = "darkgray")
lines(fit$time, fit$lower.ci, lty = 2, lwd = 1, col = "gray")
lines(fit$time, fit$upper.ci, lty = 2, lwd = 1, col = "gray")
points(tvec,  hHOs, type = "l", xlab = "Time", ylab = "Predictive Hazard", 
       cex.axis = 1.5, cex.lab = 1.5, lwd =2, lty = 1)
legend("bottomright", legend = c("HO", "BS"), lty = c(1,2),  
       lwd = c(2,2), col = c("black","darkgray"))

Christen, J. A., and F. J. Rubio. 2024. “On Harmonic Oscillator Hazard Functions.” arXiv Preprint.

On harmonic oscillator hazard functions

F. Javier Rubio and J. Andres Christen

28 August, 2024

Harmonic oscillator hazard function

The damped harmonic oscillator model

Data preparation

Maximum Likelihood Analysis

Weibull Model

The Power Generalised Weibull Model

The Damped Harmonic Oscillator Model

Comparison

Bayesian Analysis

Weibull

The Power Generalised Weibull Model

The Harmonic Oscillator hazard Model

Predictive Hazard functions

Predictive Survival functions

Comparison with a spline-based estimator