Introduction

The splineCox package provides functions for fitting spline-based Cox regression models. These models allow for flexible baseline hazard shapes and efficient model selection based on log-likelihood. The package supports predefined baseline hazard shapes as well as user-defined numeric vectors, which are normalized to have an L1 norm of 1.

Loading the Package

library(splineCox)
library(joint.Cox)  # Required for example data

## Warning: package 'joint.Cox' was built under R version 4.2.3

## Loading required package: survival

Example Dataset

The dataOvarian dataset from the joint.Cox package contains time-to-event data, event indicators, and covariates for ovarian cancer patients.

# Load the dataset
data(dataOvarian)

# Display the first few rows
head(dataOvarian)

##           t.event event t.death death group     CXCL12
## GSM432220    1650     0    1650     0     4  1.3059416
## GSM432221      30     1      30     1     4  1.2862164
## GSM432222     720     0     720     0     4 -1.3690315
## GSM432223     450     1     780     0     4  1.6132696
## GSM432224     510     1     990     1     4  0.6115144
## GSM432225    1110     0    1110     0     4 -0.5214953

Fitting the Model with Predefined Shapes

We fit a spline-based Cox regression model using three predefined baseline hazard shapes: “constant”, “increase”, and “decrease”.

# Define variables
t.event <- dataOvarian$t.event
event <- dataOvarian$event
Z <- dataOvarian$CXCL12
M <- c("constant", "increase", "decrease")

# Fit the model
reg2 <- splineCox.reg2(t.event, event, Z, model = M, plot = TRUE)

# Display the results
print(reg2)

## $model
## [1] "constant"
## 
## $parameter
## [1] 0.125 0.250 0.250 0.250 0.125
## 
## $beta
##   estimate         SE      Lower      Upper 
## 0.21342140 0.04250986 0.13010207 0.29674073 
## 
## $gamma
##  estimate        SE     Lower     Upper 
## 4.8603503 0.2037146 4.4610697 5.2596309 
## 
## $loglik
## LogLikelihood           AIC           BIC 
##     -4603.751      9211.501      9221.323 
## 
## $other_loglik
##          Loglikelihodd
## increase     -4629.807
## decrease     -4611.546

Fitting the Model with Custom Numeric Vectors

The package also allows users to specify custom numeric vectors to define the baseline hazard shape. These vectors will be normalized to have an L1 norm of 1.

# Define custom numeric vectors for baseline hazard shapes
custom_models <- list(c(0.1, 0.2, 0.3, 0.2, 0.2), c(0.2, 0.3, 0.3, 0.1, 0.1))

# Fit the model
reg2_custom <- splineCox.reg2(t.event, event, Z, model = custom_models, plot = TRUE)

# Display the results
print(reg2_custom)

## $model
## [1] 0.2 0.3 0.3 0.1 0.1
## 
## $parameter
## [1] 0.2 0.3 0.3 0.1 0.1
## 
## $beta
##   estimate         SE      Lower      Upper 
## 0.20680947 0.04245562 0.12359645 0.29002249 
## 
## $gamma
##  estimate        SE     Lower     Upper 
## 3.4358301 0.1440085 3.1535735 3.7180866 
## 
## $loglik
## LogLikelihood           AIC           BIC 
##     -4601.307      9206.615      9216.436 
## 
## $other_loglik
##                            Loglikelihodd
## c(0.1, 0.2, 0.3, 0.2, 0.2)     -4611.873

Interpreting Results

The output of the model includes: - The best-fitting baseline hazard shape or normalized custom vector. - Estimates for the regression coefficients (beta) and the baseline hazard scale parameter (gamma). - Log-likelihood for model selection. - A plot of the estimated baseline hazard function with 95% confidence intervals (if plot = TRUE).

Below are the results from the predefined shapes example:

# Print a summary of the results
print(reg2)

## $model
## [1] "constant"
## 
## $parameter
## [1] 0.125 0.250 0.250 0.250 0.125
## 
## $beta
##   estimate         SE      Lower      Upper 
## 0.21342140 0.04250986 0.13010207 0.29674073 
## 
## $gamma
##  estimate        SE     Lower     Upper 
## 4.8603503 0.2037146 4.4610697 5.2596309 
## 
## $loglik
## LogLikelihood           AIC           BIC 
##     -4603.751      9211.501      9221.323 
## 
## $other_loglik
##          Loglikelihodd
## increase     -4629.807
## decrease     -4611.546

And here are the results from the custom numeric vectors example:

# Print a summary of the results
print(reg2_custom)

## $model
## [1] 0.2 0.3 0.3 0.1 0.1
## 
## $parameter
## [1] 0.2 0.3 0.3 0.1 0.1
## 
## $beta
##   estimate         SE      Lower      Upper 
## 0.20680947 0.04245562 0.12359645 0.29002249 
## 
## $gamma
##  estimate        SE     Lower     Upper 
## 3.4358301 0.1440085 3.1535735 3.7180866 
## 
## $loglik
## LogLikelihood           AIC           BIC 
##     -4601.307      9206.615      9216.436 
## 
## $other_loglik
##                            Loglikelihodd
## c(0.1, 0.2, 0.3, 0.2, 0.2)     -4611.873

Modeling Dependence with B-spline Copulas

The splineCox package (version 0.0.4 and later) provides the spline.copula function, which implements a flexible B-spline copula model based on the five-parameter M-spline basis. This allows users to model dependence structures between two variables on the unit square \([0, 1]^2\), supporting both the copula density and distribution function.

library(ggplot2)

N <- 50
u <- v <- seq(from = 0, to = 1, length.out = N)
U <- rep(u, N)
V <- rep(v, each = N)

# Positive Exchangeable
c.data <- data.frame(
  U = U, V = V,
  C = spline.copula(U, V, R = "PE1", density = TRUE, mat = FALSE)
)

ggplot(aes(x=U, y=V), data=c.data) +
  geom_contour(aes(x=U,y=V,z=C,colour=after_stat(level)),
               data=c.data,bins=25)+xlab("u")+ylab("v")

# Negative Exchangeable
c.data <- data.frame(
  U = U, V = V,
  C = spline.copula(U, V, R = "NE3", density = TRUE, mat = FALSE)
)

ggplot(aes(x=U, y=V), data=c.data) +
  geom_contour(aes(x=U,y=V,z=C,colour=after_stat(level)),
               data=c.data,bins=25)+xlab("u")+ylab("v")