Merge pull request #131 from adibender/issue-129

Issue 129

.gitignore

@@ -18,3 +18,4 @@ vignettes/*_files
 *.dropbox.attr
 
 docs/*
+pkgdown/

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: pammtools
 Title: Piece-Wise Exponential Additive Mixed Modeling Tools for Survival Analysis
-Version: 0.1.18
+Version: 0.2.0
 Date: 2020-01-30
 Authors@R: c(
     person("Andreas", "Bender", , "[email protected]", role = c("aut", "cre"), comment=c(ORCID = "0000-0001-5628-8611")),
@@ -20,15 +20,15 @@ Imports:
     checkmate,
     magrittr,
     rlang,
-    tidyr (>= 0.8.3),
+    tidyr (>= 1.0.0),
     ggplot2,
     dplyr (>= 0.7.0),
     purrr (>= 0.2.3),
     tibble,
-    msm,
     lazyeval,
     Formula,
-    mvtnorm
+    mvtnorm,
+    pec
 Suggests:
     RColorBrewer,
     scam,

NAMESPACE

@@ -2,9 +2,12 @@
 
 S3method(arrange,nested_fdf)
 S3method(arrange,ped)
+S3method(as.data.frame,crps)
 S3method(as_ped,data.frame)
 S3method(as_ped,list)
 S3method(as_ped,nested_fdf)
+S3method(as_ped,pamm)
+S3method(as_ped,ped)
 S3method(fill,nested_fdf)
 S3method(fill,ped)
 S3method(filter,nested_fdf)
@@ -37,6 +40,9 @@ S3method(mutate,nested_fdf)
 S3method(mutate,ped)
 S3method(nest_tdc,default)
 S3method(nest_tdc,list)
+S3method(plot,pamm)
+S3method(predictSurvProb,pamm)
+S3method(print,pamm)
 S3method(rename,nested_fdf)
 S3method(rename,ped)
 S3method(right_join,nested_fdf)
@@ -56,6 +62,7 @@ S3method(summarise,nested_fdf)
 S3method(summarise,ped)
 S3method(summarize,nested_fdf)
 S3method(summarize,ped)
+S3method(summary,pamm)
 S3method(tidy_fixed,coxph)
 S3method(tidy_fixed,gam)
 S3method(transmute,nested_fdf)
@@ -103,11 +110,13 @@ export(gg_tensor)
 export(group_by)
 export(inner_join)
 export(int_info)
+export(is.pamm)
 export(is.ped)
 export(left_join)
 export(make_newdata)
 export(mutate)
 export(nest_tdc)
+export(pamm)
 export(ped_info)
 export(rename)
 export(right_join)
@@ -160,8 +169,8 @@ importFrom(lazyeval,f_eval)
 importFrom(magrittr,"%>%")
 importFrom(mgcv,predict.bam)
 importFrom(mgcv,predict.gam)
-importFrom(msm,rpexp)
 importFrom(mvtnorm,rmvnorm)
+importFrom(pec,predictSurvProb)
 importFrom(purrr,compose)
 importFrom(purrr,cross)
 importFrom(purrr,cross_df)
@@ -190,10 +199,12 @@ importFrom(stats,coef)
 importFrom(stats,median)
 importFrom(stats,model.frame)
 importFrom(stats,model.matrix)
+importFrom(stats,poisson)
 importFrom(stats,ppoints)
 importFrom(stats,predict)
 importFrom(stats,qnorm)
 importFrom(stats,quantile)
+importFrom(stats,rexp)
 importFrom(stats,terms)
 importFrom(stats,update)
 importFrom(stats,vcov)
@@ -204,5 +215,6 @@ importFrom(tidyr,fill)
 importFrom(tidyr,fill_)
 importFrom(tidyr,gather)
 importFrom(tidyr,nest)
+importFrom(tidyr,pivot_longer)
 importFrom(tidyr,replace_na)
 importFrom(tidyr,unnest)

NEWS.md

@@ -1,3 +1,10 @@
+# pammtools 0.2.0
+* Adds a new interface for model estimation called `pamm`, which is a thin wrapper
+around `mgcv::gam` with some arguments pre-set.
+* Adds S3 methods `predict.pamm` and `predictSurvProb.pamm`
+* Adds support for model evaluation using package **`pec`**
+* Fixed bug when CIs were calculated simulation based and model contained factor variables
+
 # pammtools 0.1.15
 * Interface for specification of data transformation in `as_ped` changed. The vertical bar `|` is no longer necessary to indicate concurrent or cumulative effects
 

R/as-ped.R

@@ -41,7 +41,7 @@
 #' tumor[1:3, ] %>% as_ped(Surv(days, status)~ age + sex)
 #' @return A data frame class \code{ped} in piece-wise exponential data format.
 #' @export
-as_ped <- function(data, formula, ...) {
+as_ped <- function(data, ...) {
   UseMethod("as_ped", data)
 }
 
@@ -163,3 +163,35 @@ as_ped.list <- function(
 #' @param x any R object.
 #' @export
 is.ped <- function(x) inherits(x, "ped")
+
+
+#' @rdname as_ped
+#' @param newdata A new data set (\code{data.frame}) that contains the same
+#' variables that were used to create the PED object (code{data}).
+#' @export
+as_ped.ped <- function(data, newdata, ...) {
+
+  if (is.ped(newdata)) {
+    stop("newdata already in ped format.")
+  }
+
+  trafo_args <- attr(data, "trafo_args")
+  trafo_args[["data"]] <- newdata
+  do.call(split_data,  trafo_args)
+
+}
+
+
+
+#' @rdname as_ped
+#' @export
+as_ped.pamm <- function(data, newdata, ...) {
+
+  if (is.ped(newdata)) {
+    stop("newdata already in ped format.")
+  }
+  trafo_args      <- data[["trafo_args"]]
+  trafo_args$data <- newdata
+  do.call(split_data, trafo_args)
+
+}

R/convenience-plots.R

@@ -4,7 +4,7 @@
 #' smooth terms contained in the model. If more than one smooth is present, the
 #' different smooth are faceted.
 #'
-#' @param x A data frame or object of class \code{pamm}.
+#' @param x A data frame or object of class \code{ped}.
 #' @param ... Further arguments passed to \code{\link{get_terms}}
 #' @import ggplot2
 #' @return A \code{\link[ggplot2]{ggplot}} object.

R/formula-specials.R

@@ -268,11 +268,7 @@ add_concurrent <- function(ped, data, id_var) {
     ccr_vars_i <- c(tz_var_i, tdc_vars_i)
     ccr_i_df   <- data %>%
       select(one_of(c(id_var, ccr_vars_i)))
-    if(tidyr_new_interface()) {
-      ccr_i_df <- ccr_i_df %>% unnest(cols = -one_of(id_var))
-    } else {
-      ccr_i_df <- ccr_i_df %>% unnest()
-    }
+    ccr_i_df <- ccr_i_df %>% unnest(cols = -one_of(id_var))
 
     li <- map2(ped_split, split(ccr_i_df, f = ccr_i_df[[id_var]]),
       function(.x, .y) {

R/interval-information.R

@@ -81,7 +81,7 @@ int_info.ped <- function(x, ...) {
 #' @keywords internal
 int_info.pamm <- function(x, ...) {
 
-  int_info(x$breaks)
+  int_info(x[["trafo_args"]][["cut"]])
 
 }
 
@@ -97,24 +97,22 @@ int_info.pamm <- function(x, ...) {
 #' @import dplyr
 #' @return A \code{data.frame} containing information on intervals in which
 #' values of \code{times} fall.
+#' @seealso \code{\link[base]{findInterval}} \code{\link{int_info}}
+#' @rdname get_intervals
+#' @export
 #' @examples
 #' set.seed(111018)
 #' brks <- c(0, 4.5, 5, 10, 30)
 #' int_info(brks)
 #' x <- runif (3, 0, 30)
 #' x
 #' get_intervals(brks, x)
-#'
-#' @seealso \code{\link[base]{findInterval}} \code{\link{int_info}}
-#' @rdname get_intervals
-#' @export
 get_intervals <- function(x, times, ...) {
   UseMethod("get_intervals", x)
 }
 
-#' @inherit get_intervals
-#' @inheritParams base::findInterval
 #' @rdname get_intervals
+#' @inheritParams base::findInterval
 #' @export
 get_intervals.default <- function(
   x,

R/model-evaluation.R

@@ -0,0 +1,22 @@
+#' Transform crps object to data.frame
+#'
+#' A\code{as.data.frame} S3 method for objects of class \code{\link[pec]{crps}}.
+#'
+#' @inheritParams base::as.data.frame
+#' @param x An object of class \code{crps}. See \code{\link[pec]{crps}}.
+#' @importFrom tidyr pivot_longer
+#'
+#' @export
+as.data.frame.crps <- function(x, row.names = NULL, optional = FALSE, ...) {
+
+  m <- matrix(x, nrow = dim(x)[1], ncol = dim(x)[2])
+  colnames(m) <- attr(x, "dimnames")[[2]]
+
+  m <- as.data.frame(m)
+  m$method <- attr(x, "dimnames")[[1]]
+
+  m <- m %>%
+    pivot_longer(cols = -.data$method, values_to = "IBS") %>%
+    dplyr::rename(time = .data$name)
+
+}

R/nest-utils.R

@@ -45,18 +45,12 @@ nest_tdc.default <- function(data, formula, ...) {
     vars_to_exclude <- intersect(colnames(data), tdc_vars)
     return(data %>% select(-one_of(vars_to_exclude)))
   } else {
-    df_list <-
-        if (tidyr_new_interface()) {
-          map(
-            tdc_vars,
-            ~ tidyr::nest(.data = data[, c(id, .x)], {{.x}} := one_of(.x)))
-        } else {
-          map(
-            tdc_vars,
-            ~tidyr::nest(data = data[, c(id, .)], -one_of(id), .key = !!.))
-      }
-      suppressMessages(nested_df <- df_list %>% reduce(left_join)) # better: numeric vectors in each list element
-      class(nested_df) <- c("nested_fdf", class(nested_df))
+    df_list <- map(
+      tdc_vars,
+      ~ tidyr::nest(.data = data[, c(id, .x)], {{.x}} := one_of(.x)))
+
+    suppressMessages(nested_df <- df_list %>% reduce(left_join)) # better: numeric vectors in each list element
+    class(nested_df) <- c("nested_fdf", class(nested_df))
   }
 
   nested_df %>% as_tibble()

R/pammfit.R

@@ -0,0 +1,113 @@
+unpam <- function(pamm) {
+  class(pamm) <- class(pamm)[-1]
+  pamm
+}
+repam <- function(x) {
+  class(x) <- c("pamm", class(x))
+  x
+}
+
+append_ped_attr <- function(pamm, ped) {
+
+  attr_ped <- ped_attr(ped)
+  pamm[["attr_ped"]] <- attr_ped
+
+  pamm
+
+}
+
+
+#' Fit a piece-wise exponential additive model
+#'
+#' A thin wrapper around \code{\link[mgcv]{gam}}, however, some arguments are
+#' prespecified:
+#' \code{family=poisson()}, \code{offset=data$offset} and \code{method="REML"}.
+#' The first two can not be overriden. The \code{method} argument
+#' can be specified as usual, but defaults to \code{GCV.cp} in \code{\link[mgcv]{gam}}.
+#'
+#' @inheritParams mgcv::gam
+#' @param ... Further arguments passed to \code{engine}.
+#' @param trafo_args A named list. If data is not in PED format, \code{as_ped}
+#' will be called internally with arguments provided in \code{trafo_args}.
+#' @param engine Character name of the function that will be called to fit the
+#' model. The intended entries are either \code{"gam"} or \code{"bam"}
+#' (both from package \code{mgcv}).
+#' @import mgcv
+#' @importFrom stats poisson
+#' @rdname pamm
+#' @seealso \code{\link[mgcv]{gam}}
+#' @examples
+#' ped <- tumor %>%
+#'  as_ped(Surv(days, status) ~ complications, cut = seq(0, 3000, by = 50))
+#' pam <- pamm(ped_status ~ s(tend) + complications, data = ped)
+#' summary(pam)
+#' ## Alternatively
+#' pamm(
+#'  ped_status ~ s(tend) + complications,
+#'  data = tumor,
+#' trafo_args = list(formula = Surv(days, status)~complications))
+#' @export
+pamm <- function(
+  formula,
+  data       = list(),
+  method     = "REML",
+  ...,
+  trafo_args = NULL,
+  engine     = "gam") {
+
+  dots <- list(...)
+  dots$formula <- formula
+  dots$family  <- poisson()
+  if (!is.null(trafo_args)) {
+    trafo_args$data <- data
+    data <- do.call(split_data, trafo_args)
+  }
+  dots$data   <- data
+  dots$offset <- data$offset
+
+  pamm_fit        <- do.call(engine, dots)
+  class(pamm_fit) <- c("pamm", class(pamm_fit))
+  pamm_fit        <- append_ped_attr(pamm_fit, data)
+  pamm_fit[["trafo_args"]] <- attr(data, "trafo_args")
+
+  pamm_fit
+
+}
+
+
+#' Check if object is of class pamm
+#'
+#' @param x Any R object.
+#' @rdname pamm
+#' @keywords internal
+#' @export
+is.pamm <- function(x) inherits(x, "pamm")
+
+
+#' @rdname pamm
+#' @keywords internal
+#' @export
+print.pamm <- function(x, ...) {
+
+  print(unpam(x), ...)
+
+}
+
+#' @rdname pamm
+#' @param object An object of class \code{pamm} as returned by \code{\link{pamm}}.
+#' @keywords internal
+#' @export
+summary.pamm <- function(object, ...) {
+
+  summary(unpam(object), ...)
+
+}
+
+#' @rdname pamm
+#' @keywords internal
+#' @export
+plot.pamm <- function(x, ...) {
+
+  plot(unpam(x), ...)
+
+}