Introduction to splineCox
Ren Teranishi
2025-07-11
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.
Key Features
- Computes the copula density \(c(u, v)\) or the copula
distribution function \(C(u,
v)\).
- Uses a 5×5 coefficient matrix (
R) with various built-in
presets (e.g., independence, positive/negative dependence, tail
dependence, etc.).
- Integrates with the
joint.Cox package’s
M.spline and I.spline basis functions for
computation.
Example: Visualizing Copula Densities
Below is an example illustrating how to visualize the density of the
independence copula and a positively dependent copula.
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")
Computing Dependence Measures
The spline.copula() function can also compute Kendall’s
tau and Spearman’s rho analytically from the coefficient matrix \(R\).
These are widely used rank-based measures of dependence.
# Compute Kendall's tau and Spearman's rho for preset R matrix "PE1"
out <- spline.copula(U, V, R = "PE1", Kendall = TRUE, Spearman = TRUE)
# Display the results
out$Kendall_tau
## [1] 0.5238903
## [1] 0.72