Merge pull request #135 from adibender/issue-134

Fixes #134

.Rbuildignore

@@ -19,3 +19,4 @@
 ^codecov\.yml$
 ^sandbox$
 ^cran-comments\.md$
+^pkgdown$

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: pammtools
 Title: Piece-Wise Exponential Additive Mixed Modeling Tools for Survival Analysis
-Version: 0.1.17
-Date: 2019-12-12
+Version: 0.1.18
+Date: 2020-01-30
 Authors@R: c(
     person("Andreas", "Bender", , "[email protected]", role = c("aut", "cre"), comment=c(ORCID = "0000-0001-5628-8611")),
     person("Fabian", "Scheipl", , "[email protected]", role = c("aut")))
@@ -45,5 +45,5 @@ License: MIT + file LICENSE
 LazyData: true
 URL: https://github.com/adibender/pammtools
 BugReports: https://github.com/adibender/pammtools/issues
-RoxygenNote: 7.0.1
+RoxygenNote: 7.0.2
 Encoding: UTF-8

R/add-functions.R

@@ -334,9 +334,8 @@ get_cumu_hazard <- function(
       if (ci_type == "sim") {
         newdata <- get_hazard(newdata, object, type = "response", ci = FALSE,
           time_var = time_var, ...)
-        ci_df <- split(newdata, group_indices(newdata)) %>%
+        newdata <- split(newdata, group_indices(newdata)) %>%
           map_dfr(get_sim_ci_cumu, object = object, nsim = nsim)
-        newdata <- newdata %>% bind_cols(ci_df)
       }
     }
   } else {
@@ -449,9 +448,8 @@ get_surv_prob <- function(
       if (ci_type == "sim") {
         newdata <- get_hazard(newdata, object, type = "response", ci = FALSE,
           time_var = time_var)
-        ci_df <- split(newdata, group_indices(newdata)) %>%
+        newdata <- split(newdata, group_indices(newdata)) %>%
           map_dfr(get_sim_ci_surv, object = object, nsim = nsim)
-        newdata <- newdata %>% bind_cols(ci_df)
       }
     }
   } else {
@@ -510,9 +508,8 @@ add_ci <- function(
             map_dfr(add_delta_ci, object = object, se_mult = se_mult)
       } else {
         if (ci_type == "sim") {
-          ci_df <- split(newdata, group_indices(newdata)) %>%
-            map_dfr(get_sim_ci, object = object)
-          newdata <- newdata %>% bind_cols(ci_df)
+          newdata <- split(newdata, group_indices(newdata)) %>%
+            map_dfr(get_sim_ci, object = object, nsim = nsim)
         }
       }
     }
@@ -576,10 +573,10 @@ get_sim_ci <- function(newdata, object, alpha = 0.05, nsim = 100L) {
   sim_coef_mat <- mvtnorm::rmvnorm(nsim, mean = coefs, sigma = V)
   sim_fit_mat <- apply(sim_coef_mat, 1, function(z) exp(X %*% z))
 
-  ci_lower <- apply(sim_fit_mat, 1, quantile, probs = alpha / 2)
-  ci_upper <- apply(sim_fit_mat, 1, quantile, probs = 1 - alpha / 2)
+  newdata$ci_lower <- apply(sim_fit_mat, 1, quantile, probs = alpha / 2)
+  newdata$ci_upper <- apply(sim_fit_mat, 1, quantile, probs = 1 - alpha / 2)
 
-  tibble(ci_lower = ci_lower, ci_upper = ci_upper)
+  newdata
 
 }
 
@@ -590,14 +587,14 @@ get_sim_ci_cumu <- function(newdata, object, alpha = 0.05, nsim = 100L) {
   V     <- object$Vp
   coefs <- coef(object)
 
-  sim_coef_mat <- mvtnorm::rmvnorm(nsim, mean = coefs, sigm = V)
+  sim_coef_mat <- mvtnorm::rmvnorm(nsim, mean = coefs, sigma = V)
   sim_fit_mat <- apply(sim_coef_mat, 1, function(z)
     cumsum(newdata$intlen * exp(X %*% z)))
 
-  cumu_lower <- apply(sim_fit_mat, 1, quantile, probs = alpha / 2)
-  cumu_upper <- apply(sim_fit_mat, 1, quantile, probs = 1 - alpha / 2)
+  newdata$cumu_lower <- apply(sim_fit_mat, 1, quantile, probs = alpha / 2)
+  newdata$cumu_upper <- apply(sim_fit_mat, 1, quantile, probs = 1 - alpha / 2)
 
-  tibble(cumu_lower = cumu_lower, cumu_upper = cumu_upper)
+  newdata
 
 }
 
@@ -611,9 +608,9 @@ get_sim_ci_surv <- function(newdata, object, alpha = 0.05, nsim = 100L) {
   sim_fit_mat <- apply(sim_coef_mat, 1, function(z)
     exp(-cumsum(newdata$intlen * exp(X %*% z))))
 
-  surv_lower <- apply(sim_fit_mat, 1, quantile, probs = alpha / 2)
-  surv_upper <- apply(sim_fit_mat, 1, quantile, probs = 1 - alpha / 2)
+  newdata$surv_lower <- apply(sim_fit_mat, 1, quantile, probs = alpha / 2)
+  newdata$surv_upper <- apply(sim_fit_mat, 1, quantile, probs = 1 - alpha / 2)
 
-  tibble(surv_lower = surv_lower, surv_upper = surv_upper)
+  newdata
 
 }

tests/testthat/test-add-functions.R

@@ -12,6 +12,9 @@ bam <- mgcv::bam(ped_status ~ s(tend, k = 5) + trt, data = ped,
 pem <- glm(ped_status ~ 0 + interval + trt, data = ped,
   family = poisson(), offset = offset)
 
+pam3 <- mgcv::gam(ped_status ~ s(tend, k=5, by = as.factor(trt)) + as.factor(trt),
+  data = ped, family = poisson(), offset = offset)
+
 test_that("hazard functions work for PAM", {
 
   expect_data_frame(haz <- add_hazard(ped_info(ped), bam), nrows = 5L,
@@ -53,6 +56,16 @@ test_that("hazard functions work for PAM", {
   hr2 <- ped_info(ped) %>%  add_hazard(pam2, reference = list(age = mean(.$age)))
   expect_equal(hr2$hazard, rep(1, 5))
 
+  ## factor group variable
+  ndf <- ped %>% make_newdata(tend = unique(tend), trt = unique(trt)) %>%
+    group_by(trt)
+  ndf1 <- ndf %>% add_cumu_hazard(pam3, ci = TRUE, ci_type = "default")
+  ndf2 <- ndf %>% add_cumu_hazard(pam3, ci = TRUE, ci_type = "delta")
+  ndf3 <- ndf %>% add_cumu_hazard(pam3, ci = TRUE, ci_type = "sim", nsim = 100L)
+  expect_true(all(ndf1$cumu_hazard > ndf1$cumu_lower & ndf1$cumu_hazard < ndf1$cumu_upper))
+  expect_true(all(ndf2$cumu_hazard > ndf2$cumu_lower & ndf2$cumu_hazard < ndf2$cumu_upper))
+  expect_true(all(ndf3$cumu_hazard > ndf3$cumu_lower & ndf3$cumu_hazard < ndf3$cumu_upper))
+
 })
 
 test_that("hazard functions work for PEM", {