coxlps {blapsr}

Fit a Cox proportional hazards regression model with Laplace-P-splines.

Description

Fits a Cox proportional hazards regression model for right censored data by combining Bayesian P-splines and Laplace approximations.

Usage

coxlps(formula, data, K = 30, penorder = 2, tmax = NULL)

Arguments

formula

A formula object where the ~ operator separates the response from the covariates. In a Cox model, it takes the form response ~ covariates, where response is a survival object returned by the Surv function of the survival package.

data

Optional. A data frame to match the variable names provided in formula.

K

A positive integer specifying the number of cubic B-spline functions in the basis. Default is K = 30 and allowed values are 10 <= K <= 60.

penorder

The penalty order used on finite differences of the coefficients of contiguous B-splines. Can be either 2 for a second-order penalty (the default) or 3 for a third-order penalty.

tmax

A user-specified value for the upper bound of the B-spline basis. The default is NULL, so that the B-spline basis is specified in the interval [0, tup], where tup is the upper bound of the follow-up times. It is required that tmax > tup.

Details

The log-baseline hazard is modeled as a linear combination of K cubic B-splines as obtained from cubicbs. The B-spline basis is specified in the interval [0, tup], where tup is the upper bound of the follow-up times, i.e. the largest observed follow-up time. Following Jullion and Lambert (2007), a robust Gamma prior is imposed on the roughness penalty parameter. A grid-based approach is used to explore the posterior penalty space and the resulting quadrature points serve to compute the approximate (joint) marginal posterior of the latent field vector. Point and set estimates of latent field elements are obtained from a finite mixture of Gaussian densities. The routine centers the columns of the covariate matrix around their mean value for numerical stability.

Value

An object of class coxlps containing various components from the fit. Details can be found in coxlps.object. Plot of estimated smooth hazard and survival curves can be obtained using plot.coxlps. If required, estimated baseline quantities on specific time values can be obtained with coxlps.baseline.

Author(s)

Oswaldo Gressani [e-mail].

References

Cox, D.R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society: Series B (Methodological) 34(2): 187-202.

Gressani, O. and Lambert, P. (2018). Fast Bayesian inference using Laplace approximations in a flexible promotion time cure model based on P-splines. Computational Statistical & Data Analysis 124: 151-167.

Jullion, A. and Lambert, P. (2007). Robust specification of the roughness penalty prior distribution in spatially adaptive Bayesian P-splines models. Computational Statistical & Data Analysis 51(5): 2542-2558.

See Also

Surv, coxph, simsurvdata

Examples

### Example 1 (Simulated survival data)

set.seed(3)

# Simulate survival data  with simsurvdata
betas <- c(0.13, 0.52, 0.30)
simul <- simsurvdata(a = 3.8, b = 2.2, n = 250, betas = betas , censperc = 20)
simul
simdat <- simul$survdata
plot(simul) # Plot survival data

# Estimation with coxlps
fit <- coxlps(Surv(time, delta) ~ x1 + x2 + x3, data = simdat, K = 15)
# Compare coxlps and coxph
fit
summary(coxph(Surv(time, delta) ~ x1 + x2 + x3, data = simdat))

# Fitted baseline survival vs target
plot(fit, h0 = FALSE, cred.int = 0.95, overlay.km = TRUE)
domt <- seq(0, 4, length = 100)
lines(domt, simul$S0(domt), type = "l", col = "red")
legend("topright", col=c("black", "blue", "red"), lty = rep(1,3),
      c("Bayesian LPS", "Kaplan-Meier", "Target"), cex = 0.8, bty = "n")

### Example 2 (Kidney transplant data)

data(kidneytran)
Surv.obj <- Surv(kidneytran$time, kidneytran$delta)
fit <- coxlps(Surv.obj ~ age + gender + race, data = kidneytran)
coxphfit <- coxph(Surv.obj ~ age + gender + race, data = kidneytran)
## Compare coxph and coxlps results
summary(coxphfit)
fit
## Plot Kaplan-Meier curve vs Laplace-P-spline fit
plot(fit, h0 = FALSE, overlay.km = TRUE, plot.cred = FALSE)

### Example 3 (Laryngeal cancer data)

data(laryngeal)
fit <- coxlps(Surv(time, delta) ~ age + diagyr + as.factor(stage),
               data = laryngeal)
coxphfit <- coxph(Surv(time, delta) ~ age + diagyr + as.factor(stage),
                  data = laryngeal)
## Compare coxph and coxlps results
summary(coxphfit)
fit
## Plot Kaplan-Meier curve vs Laplace-P-spline fit
plot(fit, h0 = FALSE, overlay.km = TRUE, plot.cred = FALSE)