First commit

.Rbuildignore

@@ -0,0 +1,4 @@
+^.*\.Rproj$
+^\.Rproj\.user$
+^man/figures$
+^README\.Rmd$

.gitignore

@@ -0,0 +1,3 @@
+.Rproj.user
+mixcurelps.Rproj
+README.Rmd

DESCRIPTION

@@ -0,0 +1,24 @@
+Package: mixcurelps
+Type: Package
+Title: Approximate Bayesian Inference in Mixture Cure Models
+Version: 0.1.1
+Depends: R (>= 4.1)
+Authors@R: person("Oswaldo","Gressani", 
+  email="oswaldo_gressani@hotmail.fr", role=c("aut","cre"))
+Maintainer: Oswaldo Gressani <oswaldo_gressani@hotmail.fr>
+Description: This package implements Laplace approximations and P-splines for
+  fast approximate Bayesian inference in mixture cure models.
+License: file LICENSE
+Encoding: UTF-8
+Imports: 
+    Rcpp (>= 1.0.7),
+    survival  (>= 3.2-11),
+    ggplot2   (>= 3.3.5),
+    progress  (>= 1.2.2),
+    crayon    (>= 1.4.1),
+    survminer (>= 0.4.3)
+LazyData: true
+RoxygenNote: 7.1.1
+LinkingTo: 
+    RcppArmadillo,
+    Rcpp

LICENSE

@@ -0,0 +1,2 @@
+YEAR: 2021
+COPYRIGHT HOLDER: Oswaldo Gressani

NAMESPACE

@@ -0,0 +1,12 @@
+# Generated by roxygen2: do not edit by hand
+
+S3method(plot,simdatmixcure)
+S3method(print,lpsmc)
+export(curefit)
+export(lpsmc)
+export(postpendist)
+export(simdatmixcure)
+export(simlpsmc)
+export(survcurve)
+importFrom(Rcpp,sourceCpp)
+useDynLib(mixcurelps, .registration = TRUE)

NEWS.md

@@ -0,0 +1,5 @@
+# Newspaper for the mixcurelps package #
+
+### Version 0.1.1 ###
+
+* **2021-09-09:** Release of unstable version on Github. Version name: "Keeping CPU at low temperature".

R/RcppExports.R

@@ -0,0 +1,11 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+Rcpp_Laplace <- function(lat0, v, K, PDcorrect, Dloglik, D2loglik, Qv) {
+    .Call(`_mixcurelps_Rcpp_Laplace`, lat0, v, K, PDcorrect, Dloglik, D2loglik, Qv)
+}
+
+Rcpp_cubicBspline <- function(x, lower, upper, K) {
+    .Call(`_mixcurelps_Rcpp_cubicBspline`, x, lower, upper, K)
+}
+

R/S3_lpsmc_print.R

@@ -0,0 +1,65 @@
+#' Print a lpsmc object.
+#'
+#' @description Print method for a \code{lpsmc} object.
+#'
+#' @param x An object of class \code{lpsmc}.
+#' @param ... Further arguments to be passed to print routine.
+#'
+#'
+#' @author Oswaldo Gressani \email{oswaldo_gressani@hotmail.fr} .
+#'
+#' @export
+
+print.lpsmc <- function(x,...){
+
+  tabcolnames <- c("Estimate","sd","CI90%.low","CI90.up%",
+                   "CI95%.low","CI95%.up")
+  K <- x$K
+
+  # Table for the incidence part
+  tabincidrow <- x$p
+  tabincidout <- as.data.frame(matrix(0, nrow = tabincidrow, ncol = 6))
+  rownames(tabincidout) <- colnames(x$X)
+  colnames(tabincidout) <- tabcolnames
+  tabincidout[,1] <- x$betahat
+  tabincidout[,2] <- sqrt(diag(x$Covhat)[(K + 1):(K + tabincidrow)])
+  tabincidout[,3] <- as.numeric(x$CI90[,1])[1:tabincidrow]
+  tabincidout[,4] <- x$CI90[,2][1:tabincidrow]
+  tabincidout[,5] <- x$CI95[,1][1:tabincidrow]
+  tabincidout[,6] <- x$CI95[,2][1:tabincidrow]
+  iscolnum <- sapply(tabincidout, is.numeric)
+  tabincidout[iscolnum] <- lapply(tabincidout[iscolnum], round, 3)
+
+  # Table for the latency part
+  tablatencyrow <- x$q
+  dimlat <- x$K+x$p+x$q
+  tablatencyout <- as.data.frame(matrix(0, nrow = tablatencyrow, ncol = 6))
+  rownames(tablatencyout) <- colnames(x$Z)
+  colnames(tablatencyout) <- tabcolnames
+  tablatencyout[,1] <- x$gammahat
+  tablatencyout[,2] <- sqrt(diag(x$Covhat)[(K + tabincidrow + 1):dimlat])
+  tablatencyout[,3] <- as.numeric(x$CI90[,1])[(tabincidrow + 1):(dimlat-K)]
+  tablatencyout[,4] <- x$CI90[,2][(tabincidrow + 1):(dimlat-K)]
+  tablatencyout[,5] <- x$CI95[,1][(tabincidrow + 1):(dimlat-K)]
+  tablatencyout[,6] <- x$CI95[,2][(tabincidrow + 1):(dimlat-K)]
+  iscolnum <- sapply(tablatencyout, is.numeric)
+  tablatencyout[iscolnum] <- lapply(tablatencyout[iscolnum], round, 3)
+
+  # Print output table
+
+  cat("Fitting mixture cure model with Laplacian-P-splines \n")
+  cat(paste(rep("-",50),collapse = ""),"\n")
+  cat("Sample size: ", length(x$ftime), "\n")
+  cat("No. of B-splines: ", K, "\n")
+  cat(paste(rep("-",90),collapse = ""),"\n")
+  cat("                                  (Incidence)                    \n")
+  cat(paste(rep("-",90),collapse = ""),"\n")
+  print.table(as.matrix(tabincidout), digits = 3, justify = "left")
+  cat(paste(rep("-",90),collapse = ""),"\n")
+  cat("                                   (Latency)                     \n")
+  cat(paste(rep("-",90),collapse = ""),"\n")
+  print.table(as.matrix(tablatencyout), digits = 3, justify = "left")
+  cat(paste(rep("-",90),collapse = ""),"\n")
+  cat(paste0("'Real' elapsed time: ",x$timer, " seconds.\n"))
+
+}

R/S3_simdatmixcure_plot.R

@@ -0,0 +1,60 @@
+#' @method plot simdatmixcure
+#' @export
+
+# Plot method for an object of class simdatmixcure
+plot.simdatmixcure <- function(x, ...) {
+
+
+  # tdom <- seq(0, x$tup, length = 1000)
+
+  # graphics::plot(x$fitKM, mark.time = TRUE, mark = "x", xlab = "t",
+  #                ylab = expression(S[0](t)), main = "Baseline survival",
+  #                cex.main = 0.9)
+  # graphics::abline(v = x$plateau, lty = 2, lwd = 2, col = "orange")
+  # graphics::lines(tdom, sapply(tdom, x$S0), type = "l", col = "blue")
+  # graphics::legend("topright", lty = c(1,1,2),
+  #                  col = c("black", "blue", "orange"),
+  #                  c("Kaplan-Meier", "Weibull baseline survival",
+  #                    "Start of plateau"), bty = "n", cex = 0.8)
+
+
+  # With survminer
+  tobs <- x$tobs
+  status <- x$delta
+  dataKapM <- data.frame(tobs, status)
+  fitKapM <- survival::survfit(survival::Surv(tobs, status) ~ 1,
+                               data = dataKapM)
+  plotsurv <- survminer::ggsurvplot(fitKapM,
+                        data = dataKapM,
+                        censor.shape="x",
+                        censor.size = 5.5,
+                        size = 1,
+                        palette = "#0089FF",
+                        conf.int = TRUE,
+                        font.tickslab = c(14,"darkblue"),
+                        font.x =c(14,"black"),
+                        font.y = c(14,"black"),
+                        ggtheme = ggplot2::theme_light(),
+                        risk.table = "percentage",
+                        risk.table.col ="darkblue",
+                        legend="none",
+                        legend.title="",
+                        tables.theme = survminer::theme_cleantable()
+                        )
+  plotsurv$table <- plotsurv$table + ggplot2::theme(
+    axis.text.y = ggplot2::element_blank(),
+    axis.text.x = ggplot2::element_text(size=14),
+    axis.title.x = ggplot2::element_text(size=14)
+  )
+
+  plotsurv $plot<- plotsurv$plot + ggplot2::geom_vline(xintercept = x$plateau,
+                                 linetype = "dashed", size = 1,
+                                 colour = "#15BA57")
+
+  plotsurv
+
+
+
+
+
+}

R/breastcancer.R

@@ -0,0 +1,23 @@
+#' Breast cancer data.
+#'
+#' @docType data
+#'
+#' @description Breast cancer data from the \code{breastCancerVDX} package.
+#'
+#' @usage data(breastcancer)
+#'
+#' @format A data frame with 286 rows and 4 columns.
+#' \describe{
+#'  \item{\code{tobs}}{Distant-metastasis-free survival (in days).}
+#'  \item{\code{delta}}{Event indicator \code{1}=death or relapse, \code{0}=censored.}
+#'  \item{\code{AGE}}{Age of patients.}
+#'  \item{\code{ER}}{Estrogen receptor \code{0}="<=10fmol", \code{1}=">10fmol".}
+#' }
+#'
+#'
+#' @source \url{https://doi.org/doi:10.18129/B9.bioc.breastCancerVDX}
+#'
+#' @references  Schroeder M, Haibe-Kains B, Culhane A, Sotiriou C, Bontempi G,
+#'  Quackenbush J (2021). breastCancerVDX: Gene expression datasets published
+#'  by Wang et al. [2005] and Minn et al. [2007] (VDX). R package version 1.30.0.
+"breastcancer"

R/cpprelated.R

@@ -0,0 +1,5 @@
+## usethis namespace: start
+#' @useDynLib mixcurelps, .registration = TRUE
+#' @importFrom Rcpp sourceCpp
+## usethis namespace: end
+NULL

R/curefit.R

@@ -0,0 +1,68 @@
+#' Estimated cure proportion
+#'
+#' @description
+#' Computes the estimated cure proportion based on a mixture cure model fit
+#' with \code{lpsmc}. Both estimates and approximate 90% and 95% credible
+#' intervals are shown.
+#'
+#' @param x A lpsmc object.
+#' @param covarprofile The covariate profile on which to compute the
+#' cure proportion.
+#'
+#' @return A table with the estimated cure proportion.
+#'
+#' @author Oswaldo Gressani \email{oswaldo_gressani@hotmail.fr} .
+#'
+#' @examples
+#' ### Application on breast cancer data
+#' rm(list=ls())
+#' data("breastcancer")
+#' formula <- Surv(tobs, delta) ~ inci(AGE + ER) + late(AGE + ER)
+#' fitcancer <- lpsmc(formula = formula, data = breastcancer, K = 20)
+#' covarprofile <- matrix(c(1, 30, 1, 1, 40, 0), nrow = 2 , byrow = TRUE)
+#' fitcure <- curefit(fitcancer,covarprofile)
+#' fitcure$estimcure
+#'
+#' @export
+
+
+curefit <- function(x, covarprofile){
+
+  betahat <- x$betahat
+  phat <- x$px(betahat, covarprofile)
+  p <- x$p
+  K <- x$K
+  nprofiles <- nrow(covarprofile)
+
+  CIcure <- function(y, alpha){
+    gbhat     <- log(log(as.numeric(1 + exp(y %*% betahat))))
+    Sigmabhat <- x$Covhat[(K + 1):(K + p), (K + 1):(K + p)]
+    gradbhat  <- (x$px(betahat, y)/log(as.numeric(1 + exp(y %*% betahat)))) * y
+    qz_alpha  <- stats::qnorm(alpha * 0.5, lower.tail = FALSE)
+    postsd    <- sqrt(as.numeric(t(gradbhat) %*% Sigmabhat %*% gradbhat))
+    CIcure_alpha <- c(gbhat - qz_alpha * postsd, gbhat + qz_alpha * postsd)
+    CIcure_original <- rev(exp(-exp(CIcure_alpha)))
+    return(CIcure_original)
+  }
+
+  CIcuremat <- matrix(0, nrow = nprofiles , ncol = p + 5)
+  colnames(CIcuremat) <- c(as.character(colnames(x$X)),"1-p(x)",
+                           "CI90.low","CI90.up","CI95.low","CI95.up")
+  rownames(CIcuremat) <- paste0("x.profile", seq(nprofiles))
+  for(j in 1:nprofiles){
+   CIcuremat[j, 1:p] <- covarprofile[j, ]
+   CIcuremat[j, p+1] <- 1-phat[j]
+   CIcuremat[j, (p+2):(p+3)] <- CIcure(covarprofile[j,], 0.10)
+   CIcuremat[j, (p+4):(p+5)] <- CIcure(covarprofile[j,], 0.05)
+  }
+
+  # cat(paste(rep("-",90),collapse = ""),"\n")
+  # cat("Estimated cure proportion \n")
+  # cat(paste(rep("-",90),collapse = ""),"\n")
+  # print.table(CIcuremat, digits = 3, justify = "left")
+  # cat(paste(rep("-",90),collapse = ""),"\n")
+
+
+  outlist <- list(estimcure = CIcuremat)
+
+}

R/ecog1684.R

@@ -0,0 +1,33 @@
+#' Phase III Melanoma clinical trial.
+#'
+#' @docType data
+#'
+#' @description Melanoma data from the phase III Eastern Cooperative
+#'  Oncology Group (ECOG) two-arm clinical trial studied in
+#'  Kirkwood et al. (1996) and obtained from the \code{smcure} package.
+#'
+#' @usage data(ecog1684)
+#'
+#' @format A data frame with 284 rows and 5 columns.
+#' \describe{
+#'  \item{\code{tobs}}{Relapse-free survival (in years).}
+#'  \item{\code{delta}}{\code{1}=death or relapse, \code{0}=censored.}
+#'   \item{\code{TRT}}{Treatment: \code{0}=control,
+#'        \code{1}=Interferon alpha-2b (IFN).}
+#'  \item{\code{AGE}}{Age centered to the mean.}
+#'  \item{\code{SEX}}{\code{0}=Male, \code{1}=Female.}
+#' }
+#'
+#'
+#' @source \url{https://CRAN.R-project.org/package=smcure}
+#'
+#' @references  Kirkwood, J. M., Strawderman, M. H., Ernstoff, M. S.,
+#'  Smith, T. J., Borden, E. C. and Blum, R. H. (1996). Interferon alfa-2b
+#'  adjuvant therapy of high-risk resected cutaneous melanoma: the Eastern
+#'  Cooperative Oncology Group Trial EST 1684.
+#'  \emph{Journal of clinical oncology} \strong{14}(1): 7-17.
+#' @references Corbiere, F. and Joly, P. (2007). A SAS macro for parametric
+#'  and semiparametric mixture cure models. \emph{Computer methods and programs
+#'  in Biomedicine} \strong{85}(2): 173-180.
+#'  \url{https://doi.org/10.1016/j.cmpb.2006.10.008}
+"ecog1684"