Update example codes and descriptions

seungjae2525 · Oct 25, 2023 · 35418cd · 35418cd
1 parent f29a758
commit 35418cd
Show file tree

Hide file tree

Showing 11 changed files with 70 additions and 190 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,7 @@
 Package: RMSTSens
 Version: 1.0.0
-Title: Sensitivity analysis for unmeasured confounding in estimating the difference in restricted mean survival time
+Title: Sensitivity analysis for unmeasured confounding in estimating
+        the difference in restricted mean survival time
 Description: RMSTSens is a package aimed at providing a novel sensitivity analysis method for the RMST difference when unmeasured confounding is suspected. Given a user-specified sensitivity parameter, the sensitivity range of the difference in adjusted RMST is calculated with the percentile bootstrap confidence interval for the population sensitivity range. 
 License: MIT + file LICENSE
 URL: https://github.com/seungjae2525/RMSTSens
@@ -23,7 +24,7 @@ Authors@R:
            family = "Lee",
            role = "aut",
            comment = c(ORCID = "0000-0001-7447-7045")))
-Packaged: 2023-03-23 04:52:25 UTC; seung
+Packaged: 2023-10-25 09:59:40 UTC; seung
 Author: Seungjae Lee [aut, cre] (<https://orcid.org/0000-0001-5508-8634>),
     Woojoo Lee [aut] (<https://orcid.org/0000-0001-7447-7045>)
 Maintainer: Seungjae Lee <[email protected]>
diff --git a/R/RMSTSens.R b/R/RMSTSens.R
@@ -62,13 +62,9 @@ RMSTSens <- function(...) UseMethod("RMSTSens")
 #'
 #' @examples
 #' dat <- gbsg
-#' dat$size2 <- ifelse(dat$size <= 20, 0,
-#'                     ifelse(dat$size > 20 & dat$size <= 50, 1, 2))
-#' dat$age2 <- dat$age/100
-#' dat$er2 <- dat$er/1000
 #'
 #' ## Estimation of propensity score
-#' denom.fit <- glm(hormon~(age2)^3+(age2)^3*log(age2)+meno+factor(size2)+sqrt(nodes)+er2,
+#' denom.fit <- glm(hormon~age+meno+size+factor(grade)+nodes+pgr+er,
 #'                  data=dat, family=binomial(link="logit"))
 #' dat$Ps <- predict(denom.fit, type="response")
 #'
@@ -77,22 +73,22 @@ RMSTSens <- function(...) UseMethod("RMSTSens")
 #' results.optim <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                           level.exposed="1", ps="Ps", data=dat, methods="Optim",
 #'                           use.multicore=TRUE, n.core=2,
-#'                           lambda=1.5, tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                           lambda=1.2, tau=365.25*5, ini.par=1, verbose=FALSE)
 #' results.optim
 #'
 #' # Using approximate optimization method
 #' results.approx <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                            level.exposed="1", ps="Ps", data=dat, methods="Approx",
 #'                            use.multicore=TRUE, n.core=2,
-#'                            lambda=1.5, tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                            lambda=1.2, tau=365.25*5, ini.par=1, verbose=FALSE)
 #' results.approx
 #'
 #' ## Performing the sensitivity analysis - sensitivity range with multiple lambda
 #' # Using approximate optimization method
 #' results.approx2 <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                             level.exposed="1", ps="Ps", data=dat, methods="Approx",
 #'                             use.multicore=TRUE, n.core=2,
-#'                             lambda=c(1,1.5), tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                             lambda=c(1,1.2), tau=365.25*5, ini.par=1, verbose=FALSE)
 #' results.approx2
 #'
 #' @seealso

diff --git a/R/RMSTSens.ci.R b/R/RMSTSens.ci.R
@@ -59,13 +59,9 @@
 #'
 #' @examples
 #' dat <- gbsg
-#' dat$size2 <- ifelse(dat$size <= 20, 0,
-#'                     ifelse(dat$size > 20 & dat$size <= 50, 1, 2))
-#' dat$age2 <- dat$age/100
-#' dat$er2 <- dat$er/1000
 #'
 #' ## Estimation of propensity score
-#' denom.fit <- glm(hormon~(age2)^3+(age2)^3*log(age2)+meno+factor(size2)+sqrt(nodes)+er2,
+#' denom.fit <- glm(hormon~age+meno+size+factor(grade)+nodes+pgr+er,
 #'                  data=dat, family=binomial(link="logit"))
 #' dat$Ps <- predict(denom.fit, type="response")
 #'
@@ -74,71 +70,22 @@
 #' results.approx2 <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                             level.exposed="1", ps="Ps", data=dat, methods="Approx",
 #'                             use.multicore=TRUE, n.core=2,
-#'                             lambda=c(1,1.5), tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                             lambda=c(1,1.1), tau=365.25*5, ini.par=1, verbose=FALSE)
 #'
-#' # Additional sensitivity analysis when lambda=1.7
+#' # Additional sensitivity analysis when lambda=1.2
 #' results.approx3 <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                             level.exposed="1", ps="Ps", data=dat, methods="Approx",
 #'                             use.multicore=TRUE, n.core=2,
-#'                             lambda=c(1.7), tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                             lambda=c(1.2), tau=365.25*5, ini.par=1, verbose=FALSE)
 #'
 #' # Percentile bootstrap CI for population sensitivity range
 #' re.ap.boot <- RMSTSens.ci(x=RMSTSens.merge(x=list(results.approx2, results.approx3)),
 #'               B=50, # Set B=50 to reduce computation time for R checks
 #'               level=0.95, seed=220524,
-#'               formula=hormon~(age2)^3+(age2)^3*log(age2)+meno+factor(size2)+sqrt(nodes)+er2,
+#'               formula=hormon~age+meno+size+factor(grade)+nodes+pgr+er,
 #'               model="logistic", trunc.prop=0, use.multicore=TRUE, n.core=2)
 #' re.ap.boot
 #'
-#' ## Estimate of propensity score using random forest
-#' library(randomForest)
-#' set.seed(220528)
-#' dat$size2 <- factor(dat$size2)
-#' model.rf <- randomForest(formula=factor(hormon)~age2+meno+size2+nodes+er2, data=dat,
-#'                          ntree=1501, mtry=5)
-#' dat$Ps.rf <- as.numeric(predict(model.rf, type= "prob")[,2])
-#' # Trimming the propensity score
-#' trunc.prop <- 0.1
-#' dat$Ps.rf <- ifelse(dat$Ps.rf > quantile(dat$Ps.rf, 1-trunc.prop),
-#'                     quantile(dat$Ps.rf, 1-trunc.prop),
-#'                     ifelse(dat$Ps.rf < quantile(dat$Ps.rf, trunc.prop),
-#'                            quantile(dat$Ps.rf, trunc.prop), dat$Ps.rf))
-#' sum(dat$Ps.rf== 0) + sum(dat$Ps.rf== 1) # no 0 or 1 value
-#' # Sensitivity range for difference in adjusted RMST
-#' results.approx.rf <- RMSTSens(time="rfstime", status="status", exposure="hormon",
-#'                               level.exposed="1", ps="Ps.rf", data=dat, methods="Approx",
-#'                               use.multicore=TRUE, n.core=2,
-#'                               lambda=c(1,1.5,2), tau=365.25*5, ini.par=1, verbose=FALSE)
-#' # Percentile bootstrap CI for population sensitivity range
-#' re.rf <- RMSTSens.ci(x=results.approx.rf,
-#'                      B=50, # Set B=50 to reduce computation time for R checks
-#'                      level=0.95, seed=220528,
-#'                      formula=factor(hormon)~age2+meno+size2+nodes+er2,
-#'                      model="rf", trunc.prop=trunc.prop, use.multicore=TRUE, n.core=2,
-#'                      ntree=1501, mtry=5)
-#' re.rf
-#'
-#'
-#' ## Estimate of propensity score using generalized boosted regression modeling
-#' library(gbm)
-#' set.seed(220528)
-#' model.gbm <- gbm(formula=hormon~age2+meno+size2+nodes+er2,
-#'                  data=dat, distribution= "bernoulli", verbose= FALSE)
-#' dat$Ps.gbm <- as.numeric(predict(model.gbm, type= "response"))
-#' sum(dat$Ps.gbm== 0) + sum(dat$Ps.gbm== 1) # Check that there is 0 or 1 value
-#' # Sensitivity range for difference in adjusted RMST
-#' results.approx.gbm <- RMSTSens(time="rfstime", status="status", exposure="hormon",
-#'                                level.exposed="1", ps="Ps.gbm", data=dat, methods="Approx",
-#'                                use.multicore=TRUE, n.core=2,
-#'                                lambda=c(1,1.5,2), tau=365.25*5, ini.par=1, verbose=FALSE)
-#' # Percentile bootstrap CI for population sensitivity range
-#' re.gbm <- RMSTSens.ci(x=results.approx.gbm,
-#'                       B=50, # Set B=50 to reduce computation time for R checks
-#'                       level=0.95, seed=220528,
-#'                       formula=hormon~age2+meno+size2+nodes+er2,
-#'                       model="gbm", trunc.prop=0, use.multicore=TRUE, n.core=2)
-#' re.gbm
-#'
 #' @seealso
 #'  \code{\link[RMSTSens]{RMSTSens}}, \code{\link[RMSTSens]{print.RMSTSens}}, \code{\link[RMSTSens]{autoplot.RMSTSens}}
 #'

diff --git a/R/RMSTSens.merge.R b/R/RMSTSens.merge.R
@@ -11,13 +11,9 @@
 #'
 #' @examples
 #' dat <- gbsg
-#' dat$size2 <- ifelse(dat$size <= 20, 0,
-#'                     ifelse(dat$size > 20 & dat$size <= 50, 1, 2))
-#' dat$age2 <- dat$age/100
-#' dat$er2 <- dat$er/1000
 #'
 #' ## Estimation of propensity score
-#' denom.fit <- glm(hormon~(age2)^3+(age2)^3*log(age2)+meno+factor(size2)+sqrt(nodes)+er2,
+#' denom.fit <- glm(hormon~age+meno+size+factor(grade)+nodes+pgr+er,
 #'                  data=dat, family=binomial(link="logit"))
 #' dat$Ps <- predict(denom.fit, type="response")
 #'
@@ -26,14 +22,14 @@
 #' results.approx2 <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                             level.exposed="1", ps="Ps", data=dat, methods="Approx",
 #'                             use.multicore=TRUE, n.core=2,
-#'                             lambda=c(1,1.5,2.0), tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                             lambda=c(1,1.1,1.2), tau=365.25*5, ini.par=1, verbose=FALSE)
 #' RMSTSens.merge(x=list(results.approx2))
 #'
-#' # Additional sensitivity analysis when lambda=1.7
+#' # Additional sensitivity analysis when lambda=1.15
 #' results.approx3 <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                             level.exposed="1", ps="Ps", data=dat, methods="Approx",
 #'                             use.multicore=TRUE, n.core=2,
-#'                             lambda=c(1.7), tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                             lambda=c(1.15), tau=365.25*5, ini.par=1, verbose=FALSE)
 #' RMSTSens.merge(x=list(results.approx2, results.approx3))
 #'
 #' @export

diff --git a/R/plot.RMSTSens.R b/R/plot.RMSTSens.R
@@ -28,13 +28,9 @@
 #'
 #' @examples
 #' dat <- gbsg
-#' dat$size2 <- ifelse(dat$size <= 20, 0,
-#'                     ifelse(dat$size > 20 & dat$size <= 50, 1, 2))
-#' dat$age2 <- dat$age/100
-#' dat$er2 <- dat$er/1000
 #'
 #' ## Estimation of propensity score
-#' denom.fit <- glm(hormon~(age2)^3+(age2)^3*log(age2)+meno+factor(size2)+sqrt(nodes)+er2,
+#' denom.fit <- glm(hormon~age+meno+size+factor(grade)+nodes+pgr+er,
 #'                  data=dat, family=binomial(link="logit"))
 #' dat$Ps <- predict(denom.fit, type="response")
 #'
@@ -43,17 +39,17 @@
 #' results.approx2 <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                             level.exposed="1", ps="Ps", data=dat, methods="Approx",
 #'                             use.multicore=TRUE, n.core=2,
-#'                             lambda=c(1,1.5,2.0), tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                             lambda=c(1,1.1,1.2), tau=365.25*5, ini.par=1, verbose=FALSE)
 #' autoplot(object=results.approx2, alpha.range=0.4, alpha.ci=0.9, smooth.par=100,
 #'          ytickdiff=100, point.size=1.4, h.width=1,
 #'          axis.title.size=15, axis.text.size=12,
 #'          save.plot=FALSE, save.plot.name="Plot", save.plot.device="png",
 #'          save.plot.width=10, save.plot.height=6, save.plot.dpi=300)
-#' # Additional sensitivity analysis when lambda=1.7
+#' # Additional sensitivity analysis when lambda=1.15
 #' results.approx3 <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                             level.exposed="1", ps="Ps", data=dat, methods="Approx",
 #'                             use.multicore=TRUE, n.core=2,
-#'                             lambda=c(1.7), tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                             lambda=c(1.15), tau=365.25*5, ini.par=1, verbose=FALSE)
 #' # After merging two results, plot the analysis results.
 #' autoplot(object=RMSTSens.merge(list(results.approx2, results.approx3)), smooth.par=100,
 #'          alpha.range=0.4, alpha.ci=0.9,
@@ -66,7 +62,7 @@
 #' re.ap.boot <- RMSTSens.ci(x=RMSTSens.merge(list(results.approx2, results.approx3)),
 #'               B=50, # Set B=50 to reduce computation time for R checks
 #'               level=0.95, seed=220524,
-#'               formula=hormon~(age2)^3+(age2)^3*log(age2)+meno+factor(size2)+sqrt(nodes)+er2,
+#'               formula=hormon~age+meno+size+factor(grade)+nodes+pgr+er,
 #'               model="logistic", use.multicore=TRUE, n.core=2)
 #' autoplot(object=re.ap.boot, alpha.range=0.4, alpha.ci=0.9, smooth.par=100,
 #'          ytickdiff=100, point.size=1.4, h.width=1,
@@ -136,6 +132,23 @@ autoplot_RMSTSens <- function(xxx=NULL,
                               axis.title.size=NULL, axis.text.size=NULL,
                               save.plot=FALSE, save.plot.name=NULL, save.plot.device=NULL,
                               save.plot.width=NULL, save.plot.height=NULL, save.plot.dpi=NULL, ...) {
+  ## Function for counting the number of decimal places
+  decimalplaces <- function(x) {
+    rr <- c()
+    for(i in 1:length(x)){
+      if(i == 1){
+        rr[i] <- "1.0"
+      } else {
+        if (abs(x[i] - round(x[i])) > .Machine$double.eps^0.5) {
+          rrr <- nchar(strsplit(sub('0+$', '', as.character(x[i])), ".", fixed = TRUE)[[1]][[2]])
+        } else {
+          rrr <- 0
+        }
+        rr[i] <- sprintf(paste0("%.",rrr,"f"), x[i])
+      }
+    }
+    return(rr)
+  }
 
   ##
   xx <- xxx$result.df
@@ -190,7 +203,7 @@ autoplot_RMSTSens <- function(xxx=NULL,
         xlab(expression(bold(Lambda))) + ylab("Difference in RMST")  +
         theme_bw() +
         scale_x_continuous(breaks = xx$Lambda,
-                           labels = sprintf("%.1f", xx$Lambda), expand = c(0.005,0.005)) +
+                           labels = decimalplaces(xx$Lambda), expand = c(0.005,0.005)) +
         scale_y_continuous(breaks = ylabel, labels = ylabel) +
         theme(axis.title.x=element_text(face="bold", size=axis.title.size),
               axis.title.y=element_text(face="bold", size=axis.title.size, angle=90),
@@ -227,7 +240,7 @@ autoplot_RMSTSens <- function(xxx=NULL,
         xlab(expression(bold(Lambda))) + ylab("Difference in RMST")  +
         theme_bw() +
         scale_x_continuous(breaks = xx$Lambda,
-                           labels = sprintf("%.1f", xx$Lambda), expand = c(0.005,0.005)) +
+                           labels = decimalplaces(xx$Lambda), expand = c(0.005,0.005)) +
         scale_y_continuous(breaks = ylabel, labels = ylabel) +
         theme(axis.title.x=element_text(face="bold", size=axis.title.size),
               axis.title.y=element_text(face="bold", size=axis.title.size, angle=90),
@@ -264,7 +277,7 @@ autoplot_RMSTSens <- function(xxx=NULL,
         xlab(expression(bold(Lambda))) + ylab("Difference in RMST")  +
         theme_bw() +
         scale_x_continuous(breaks = xx$Lambda,
-                           labels = sprintf("%.1f", xx$Lambda), expand = c(0.005,0.005)) +
+                           labels = decimalplaces(xx$Lambda), expand = c(0.005,0.005)) +
         scale_y_continuous(breaks = ylabel, labels = ylabel) +
         theme(axis.title.x=element_text(face="bold", size=axis.title.size),
               axis.title.y=element_text(face="bold", size=axis.title.size, angle=90),
@@ -293,7 +306,7 @@ autoplot_RMSTSens <- function(xxx=NULL,
         xlab(expression(bold(Lambda))) + ylab("Difference in RMST")  +
         theme_bw() +
         scale_x_continuous(breaks = xx$Lambda,
-                           labels = sprintf("%.1f", xx$Lambda), expand = c(0.005,0.005)) +
+                           labels = decimalplaces(xx$Lambda), expand = c(0.005,0.005)) +
         scale_y_continuous(breaks = ylabel, labels = ylabel) +
         theme(axis.title.x=element_text(face="bold", size=axis.title.size),
               axis.title.y=element_text(face="bold", size=axis.title.size, angle=90),

diff --git a/R/print.RMSTSens.R b/R/print.RMSTSens.R
@@ -11,13 +11,9 @@
 #'
 #' @examples
 #' dat <- gbsg
-#' dat$size2 <- ifelse(dat$size <= 20, 0,
-#'                     ifelse(dat$size > 20 & dat$size <= 50, 1, 2))
-#' dat$age2 <- dat$age/100
-#' dat$er2 <- dat$er/1000
 #'
 #' ## Estimation of propensity score
-#' denom.fit <- glm(hormon~(age2)^3+(age2)^3*log(age2)+meno+factor(size2)+sqrt(nodes)+er2,
+#' denom.fit <- glm(hormon~age+meno+size+factor(grade)+nodes+pgr+er,
 #'                  data=dat, family=binomial(link="logit"))
 #' dat$Ps <- predict(denom.fit, type="response")
 #'
@@ -26,29 +22,29 @@
 #' results.optim <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                           level.exposed="1", ps="Ps", data=dat, methods="Optim",
 #'                           use.multicore=TRUE, n.core=2,
-#'                           lambda=1.5, tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                           lambda=1.2, tau=365.25*5, ini.par=1, verbose=FALSE)
 #' print(results.optim)
 #'
 #' # Using approximate optimization method
 #' results.approx <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                            level.exposed="1", ps="Ps", data=dat, methods="Approx",
 #'                            use.multicore=TRUE, n.core=2,
-#'                            lambda=1.5, tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                            lambda=1.2, tau=365.25*5, ini.par=1, verbose=FALSE)
 #' print(results.approx)
 #'
 #' ## Performing the sensitivity analysis - sensitivity range with multiple lambda
 #' # Using approximate optimization method
 #' results.approx2 <- RMSTSens(time="rfstime", status="status", exposure="hormon",
 #'                             level.exposed="1", ps="Ps", data=dat, methods="Approx",
 #'                             use.multicore=TRUE, n.core=2,
-#'                             lambda=c(1,1.5), tau=365.25*5, ini.par=1, verbose=FALSE)
+#'                             lambda=c(1,1.2), tau=365.25*5, ini.par=1, verbose=FALSE)
 #' print(results.approx2)
 #'
 #' # Percentile bootstrap CI for population sensitivity range
 #' re.ap.boot <- RMSTSens.ci(x=results.approx2,
 #'               B=50, # Set B=50 to reduce computation time for R checks
 #'               level=0.95, seed=220524,
-#'               formula=hormon~(age2)^3+(age2)^3*log(age2)+meno+factor(size2)+sqrt(nodes)+er2,
+#'               formula=hormon~age+meno+size+factor(grade)+nodes+pgr+er,
 #'               model="logistic", use.multicore=TRUE, n.core=2)
 #' print(re.ap.boot)
 #'

diff --git a/man/RMSTSens.Rd b/man/RMSTSens.Rd