version 1.1.0

DESCRIPTION

@@ -2,11 +2,11 @@ Package: intccr
 Type: Package
 Title: Semiparametric Competing Risks Regression under Interval
         Censoring
-Version: 1.0.5
+Version: 1.1.0
 Author: Giorgos Bakoyannis <gbakogia@iu.edu>, Jun Park <jp84@iu.edu>
 Maintainer: Jun Park <jp84@iu.edu>
 Description: Semiparametric regression models on the cumulative incidence function with interval-censored competing risks data as described in Bakoyannis, Yu, & Yiannoutsos (2017) <doi:10.1002/sim.7350>. The main function fits the proportional subdistribution hazards model (Fine-Gray model), the proportional odds model, and other models that belong to the class of semiparametric generalized odds rate transformation models.
-Date: 2019-01-16
+Date: 2019-01-24
 Imports: alabama (>= 2015.3.1), doParallel, foreach, parallel,
         numDeriv, stats, utils
 Depends: R (>= 2.14.0)
@@ -15,6 +15,6 @@ Encoding: UTF-8
 LazyData: true
 RoxygenNote: 6.0.1
 NeedsCompilation: no
-Packaged: 2019-01-16 17:33:34 UTC; Jun
+Packaged: 2019-01-24 20:51:14 UTC; Jun
 Repository: CRAN
-Date/Publication: 2019-01-16 18:00:03 UTC
+Date/Publication: 2019-01-24 21:20:03 UTC

MD5

@@ -1,24 +1,24 @@
-c9567f3868abff995a5ab9297264e5a3 *DESCRIPTION
-1ca814efa8634030f77a2824178d39e9 *NAMESPACE
-eb22a4010601a62acd8264f17939a06b *R/bssmle.R
-0983fa325427673b475633b919cfe65d *R/bssmle_se.R
-eeca2ee1fed13d7c14489f76acd2fc46 *R/ciregic.R
-a805c7c0052ca0c72a6b5e69bb21effc *R/dataprep.R
+51acfdce57f7052512742ce8f14cf335 *DESCRIPTION
+c78b27d5b41e93ce49ba398cbac56912 *NAMESPACE
+95b520970702426bac4363f1539c7794 *R/bssmle.R
+ae63c45ab17296b46c74d970b568d371 *R/bssmle_se.R
+2427f190fa19c63e8429ebbb4cdc8af6 *R/ciregic.R
+4cfa5e18819caeca0f7080d02959f037 *R/dataprep.R
 309784dcab7b3f0386e51f1561c64341 *R/longdata.R
-1dfea5d46a5b3e7f1c5fc2e81a0b4704 *R/naive_b.R
+07f6c4ac6f675551e67b750fe31c1a67 *R/naive_b.R
 29922bfa59014246d29a0209986289f3 *R/pseudoHIVlong.R
 bbd2a95aee0d87d39bb00b1b522ad2dc *R/simdata.R
 512bacc4ce440952aab345213db32cd0 *R/surv2.R
 53cb826850deb94719d43789216bcf4d *data/longdata.rda
 4aa3fc5663a20612431874c007b62399 *data/pseudo.HIV.long.rda
 78eae9bc04d386eec1aacef483a57454 *data/simdata.rda
 ed3ceca062c6d09bc223c791b5140141 *man/Surv2.Rd
-53c074463a912e908975166cee443d29 *man/bssmle.Rd
-83ca7f49b640d2d2e8806341cff74c9d *man/bssmle_se.Rd
-2edcaf16eb1f8618826e8bb8f40a22cb *man/ciregic.Rd
+2f5aaa27cfd74de7039e432c1746326b *man/bssmle.Rd
+93167885b0fd6f20ca57820a88800d0d *man/bssmle_se.Rd
+a6bea2abadb4bea2a716e1dee41af160 *man/ciregic.Rd
 b1888fdc4ba620f037ae47b4e9aeba98 *man/dataprep.Rd
 b0c1db2255bad6d54bba7d1a1df204ad *man/longdata.Rd
-ed5fb50ba8bb0b512f5107441e5d4dd1 *man/naive_b.Rd
+74612c6e7d96c67eb6259d7a3cddcae1 *man/naive_b.Rd
 daafbe946bfbe5ab57d87268230b927a *man/predict.ciregic.Rd
 010f3620d01d8cb44a9c0c63b23249b0 *man/pseudo.HIV.long.Rd
 92b675c17a46ab186c67a12ffb99bad5 *man/simdata.Rd

NAMESPACE

@@ -14,6 +14,7 @@ import(foreach)
 import(numDeriv)
 import(parallel)
 import(stats)
+import(tcltk)
 importFrom(alabama,constrOptim.nl)
 importFrom(doParallel,registerDoParallel)
 importFrom(splines,bs)

R/bssmle.R

@@ -5,6 +5,7 @@
 #' @param formula a formula object relating survival object \code{Surv2(v, u, event)} to a set of covariates
 #' @param data a data frame that includes the variables named in the formula argument
 #' @param alpha \eqn{\alpha = (\alpha1, \alpha2)} contains parameters that define the link functions from class of generalized odds-rate transformation models. The components \eqn{\alpha1} and \eqn{\alpha2} should both be \eqn{\ge 0}. If \eqn{\alpha1 = 0}, the user assumes the proportional subdistribution hazards model or the Fine-Gray model for the cause of failure 1. If \eqn{\alpha2 = 1}, the user assumes the proportional odds model for the cause of failure 2.
+#' @param k a tuning parameter to control the number of knots. \code{k = 1} is the default, but \eqn{0.5 \le}  \code{k} \eqn{\le 1}.
 #' @keywords bssmle
 #' @import stats
 #' @importFrom alabama constrOptim.nl
@@ -25,7 +26,7 @@
 #' est.longdata <- intccr:::bssmle(Surv2(v, u, c) ~ z1 + z2, data = newdata, alpha = c(1, 1))
 
 
-bssmle <- function(formula, data, alpha){
+bssmle <- function(formula, data, alpha, k = 1){
 
   ## Create time-window, event var and design matrix
   mf <- model.frame(formula = formula, data = data)
@@ -42,7 +43,7 @@ bssmle <- function(formula, data, alpha){
 
   ## B-spline basis matrix
   t <- c(Tv, Tu[delta > 0])
-  nk <- floor(length(t)^(1/3))
+  nk <- floor(k * length(t)^(1/3))
   max <- nk + 1
   knots <- quantile(t, seq(0, 1, by = 1 / (nk + 1)))[2:max]
   Bv <- bs(Tv, knots = knots, degree = 3, intercept = TRUE, Boundary.knots = c(min(t), max(t)))
@@ -74,7 +75,7 @@ bssmle <- function(formula, data, alpha){
   colnames(comb) <- colnames(Z)
 
   ## Define the starting values for the optimizer
-  b0 <- naive_b(data, Tv, Tu, delta, q)
+  b0 <- naive_b(data, Tv, Tu, delta, q, k)
 
   ## Define the function for the negative log-likelihood
   nLL <- function(x) {

R/bssmle_se.R

@@ -5,10 +5,11 @@
 #' @param formula a formula object relating survival object \code{Surv2(v, u, event)} to a set of covariates
 #' @param data a data frame that includes the variables named in the formula argument
 #' @param alpha \eqn{\alpha = (\alpha1, \alpha2)} contains parameters that define the link functions from class of generalized odds-rate transformation models. The components \eqn{\alpha1} and \eqn{\alpha2} should both be \eqn{\ge 0}. If \eqn{\alpha1 = 0}, the user assumes the proportional subdistribution hazards model or the Fine-Gray model for the cause of failure 1. If \eqn{\alpha2 = 1}, the user assumes the proportional odds model for the cause of failure 2.
+#' @param k a tuning parameter to control the number of knots. \code{k = 1} is the default, but \eqn{0.5 \le}  \code{k} \eqn{\le 1}.
 #' @param do.par using parallel computing for bootstrap calculation. If \code{do.par = TRUE}, parallel computing will be used during the bootstrap estimation of the variance-covariance matrix for the regression parameter estimates.
 #' @param nboot a number of bootstrap samples for estimating variances and covariances of the estimated regression coefficients. If \code{nboot = 0}, the function \code{ciregic} does dot perform bootstrap estimation of the variance matrix of the regression parameter estimates and returns \code{NA} in the place of the estimated variance matrix of the regression parameter estimates.
 #' @keywords bssmle_se
-#' @import foreach parallel numDeriv
+#' @import foreach parallel numDeriv tcltk
 #' @importFrom doParallel registerDoParallel
 #' @importFrom utils txtProgressBar setTxtProgressBar
 #' @details The function \code{bssmle_se} estimates bootstrap standard errors for the estimated regression coefficients from the function \code{bssmle}.
@@ -20,7 +21,7 @@
 #' est.vcov <- intccr:::bssmle_se(formula = Surv2(v, u, c) ~ z1 + z2, data = simdata,
 #'                                alpha = c(1, 1), do.par = FALSE, nboot = 1)
 
-bssmle_se <- function(formula, data, alpha, do.par, nboot) {
+bssmle_se <- function(formula, data, alpha, k = 1, do.par, nboot) {
   tmp <- list()
   for(i in 1:nboot){
     tmp[[i]] <- data[sample(dim(data)[1], replace = TRUE), ]
@@ -35,7 +36,7 @@ bssmle_se <- function(formula, data, alpha, do.par, nboot) {
                         utils::setTxtProgressBar(pb, j)
                         par <- bssmle(formula,
                                       data = tmp[[j]],
-                                      alpha)
+                                      alpha, k)
                         q <- length(par[[2]])
                         n <- (length(par[[1]]) - 2 * q) / 2
                         pars <- par[[1]][(2 * n + 1):(2 * n + 2 * q)]
@@ -50,12 +51,12 @@ bssmle_se <- function(formula, data, alpha, do.par, nboot) {
     res.bt <- foreach(j = 1:nboot,
                       .combine = "rbind",
                       .export = c("naive_b", "bssmle"),
-                      .packages = c("splines", "stats", "alabama", "utils")) %dopar% {
-                        pb <- utils::txtProgressBar(max = nboot, style = 3)
-                        utils::setTxtProgressBar(pb, j)
+                      .packages = c("splines", "stats", "alabama", "utils", "tcltk")) %dopar% {
+                        pb <- tkProgressBar("Parallel task", min = 1, max = nboot)
+                        setTkProgressBar(pb, j)
                         par <- bssmle(formula,
                                       data = tmp[[j]],
-                                      alpha)
+                                      alpha, k)
                         q <- length(par[[2]])
                         n <- (length(par[[1]]) - 2 * q) / 2
                         pars <- par[[1]][(2 * n + 1):(2 * n + 2 * q)]

R/ciregic.R

@@ -5,6 +5,7 @@
 #' @param formula a formula object relating the survival object \code{Surv2(v, u, event)} to a set of covariates
 #' @param data a data frame that includes the variables named in the formula argument
 #' @param alpha \eqn{\alpha = (\alpha1, \alpha2)} contains parameters that define the link functions from class of generalized odds-rate transformation models. The components \eqn{\alpha1} and \eqn{\alpha2} should both be \eqn{\ge 0}. If \eqn{\alpha1 = 0}, the user assumes the proportional subdistribution hazards model or the Fine-Gray model for the cause of failure 1. If \eqn{\alpha2 = 1}, the user assumes the proportional odds model for the cause of failure 2.
+#' @param k a tuning parameter to control the number of knots. \code{k = 1} is the default, but \eqn{0.5 \le}  \code{k} \eqn{\le 1}.
 #' @param do.par an option to use parallel computing for bootstrap. If \code{do.par = TRUE}, parallel computing will be used during the bootstrap estimation of the variance-covariance matrix for the regression parameter estimates.
 #' @param nboot a number of bootstrap samples for estimating variances and covariances of the estimated regression coefficients. If \code{nboot = 0}, the function \code{ciregic} does not perform bootstrap estimation of the variance-covariance matrix of the regression parameter estimates and returns \code{NA} in the place of the estimated variance-covariance matrix of the regression parameter estimates.
 #' @return The function \code{ciregic} provides an object of class \code{ciregic} with components:
@@ -65,14 +66,14 @@
 #' points(pred$t, pred$cif2, type = "l", col = 2)
 #'
 #' @export
-ciregic <- function(formula, data, alpha, do.par, nboot) UseMethod("ciregic")
+ciregic <- function(formula, data, alpha, k = 1, do.par, nboot) UseMethod("ciregic")
 
 #' @export
-ciregic.default <- function(formula, data, alpha, do.par, nboot){
-  est <- bssmle(formula, data, alpha)
+ciregic.default <- function(formula, data, alpha, k = 1, do.par, nboot){
+  est <- bssmle(formula, data, alpha, k)
   if(min(!is.na(est$beta)) == 1) {
     if(nboot >= 1){
-      res <- bssmle_se(formula, data, alpha, do.par, nboot)
+      res <- bssmle_se(formula, data, alpha, k, do.par, nboot)
       Sigma <- res$Sigma
       numboot <- res$numboot
       q <- length(est$varnames)
@@ -94,6 +95,7 @@ ciregic.default <- function(formula, data, alpha, do.par, nboot){
               gamma = gamma,
               vcov = Sigma,
               alpha = est$alpha,
+              k = k,
               loglikelihood = est$loglikelihood,
               convergence = est$convergence,
               tms = est$tms,
@@ -115,6 +117,7 @@ ciregic.default <- function(formula, data, alpha, do.par, nboot){
               gamma = gamma,
               vcov = Sigma,
               alpha = est$alpha,
+              k = k,
               loglikelihood = est$loglikelihood,
               convergence = est$convergence,
               tms = est$tms,

R/dataprep.R

@@ -75,13 +75,23 @@ dataprep <- function(data, ID, time, event, Z) {
   }
   colnames(X) <- Z
   temp <- data.frame(id, v, u, c, X)
-  if (sum(is.na(temp)) != 0){       # detect NA in temp
-    naval <- which(is.na(v))        # detect subject id who has NA (one right-censored time record)
-    if(length(naval) == 1) {
-      warning("subject id ", naval, " is omitted because its interval is (0, Inf).")
-    } else {
-      warning("subject id ", toString(naval), " are omitted because those intervals are (0, Inf).")
+  if (sum(temp$v == 0 & temp$u == Inf) + sum(temp$v == 0 & temp$u == 0) > 0) {
+    right1 <- which(temp$v == 0 & temp$u == Inf)       # detect subject id who has one right-censored time record
+    left1 <- which(temp$v == 0 & temp$u == 0)          # detect subject id who has one left-censored time record
+    if(length(right1) == 1) {
+      warning("subject id ", right1, " is omitted because its interval is (0, Inf).")
+    }
+    if (length(left1) == 1) {
+      warning("subject id ", left1, " is omitted because its interval is (0, 0).")
+    }
+    if (length(right1) > 1) {
+      warning("subject id ", toString(right1), " are omitted because those intervals are (0, Inf).")
+    }
+    if (length(left1) > 1) {
+      warning("subject id ", toString(left1), " are omitted because those intervals are (0, 0).")
     }
+    return(temp[!(temp$id %in% c(right1, left1)),])
+  } else {
+    return(temp)
   }
-  na.omit(temp)
 }

R/naive_b.R

@@ -7,6 +7,7 @@
 #' @param u the first observation time after the failure.
 #' @param c an indicator of cause of failure. If an observation is righ-censored, \code{event = 0}; otherwise, \code{event = 1} or \code{event = 2}, where \code{1} represents the first cause of failure, and \code{2} represents the second cause of failure. The current version of package only allows for two causes of failure.
 #' @param q a dimension of design matrix.
+#' @param k a tuning parameter to control the number of knots. \code{k = 1} is the default, but \eqn{0.5 \le}  \code{k} \eqn{\le 1}.
 #' @keywords naive_b
 #' @importFrom splines bs
 #' @details The function \code{naive_b} provides initial values for the optimization procedure.
@@ -16,9 +17,9 @@
 #' attach(simdata)
 #' intccr:::naive_b(data = simdata, v = v, u = u, c = c, q = 2)
 
-naive_b <- function(data, v, u, c, q){
+naive_b <- function(data, v, u, c, q, k = 1){
   t <- c(v, u[c > 0])
-  nk <- floor(length(t)^(1/3))
+  nk <- floor(k * length(t)^(1/3))
 
   max <- nk + 1
   knots <- quantile(t, seq(0, 1, by = 1 / (nk + 1)))[2:max]