version 1.1.0

cran · Sep 18, 2023 · 251807d · 251807d
1 parent bf60f52
commit 251807d
Show file tree

Hide file tree

Showing 9 changed files with 462 additions and 214 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: WH
 Type: Package
 Title: Enhanced Implementation of Whittaker-Henderson Smoothing
-Version: 1.0.3
+Version: 1.1.0
 Authors@R: person("Guillaume", "Biessy", email = "guillaume.biessy78@gmail.com",
   role = c("aut", "cre"), comment = c(ORCID = "0000-0003-3756-7345"))
 Description: An enhanced implementation of Whittaker-Henderson smoothing for the gradation 
@@ -21,8 +21,8 @@ Config/testthat/edition: 3
 RoxygenNote: 7.2.3
 VignetteBuilder: knitr
 NeedsCompilation: no
-Packaged: 2023-06-13 17:40:05 UTC; Guillaume BIESSY
+Packaged: 2023-09-16 11:47:09 UTC; Guillaume BIESSY
 Author: Guillaume Biessy [aut, cre] (<https://orcid.org/0000-0003-3756-7345>)
 Maintainer: Guillaume Biessy <guillaume.biessy78@gmail.com>
 Repository: CRAN
-Date/Publication: 2023-06-14 08:32:20 UTC
+Date/Publication: 2023-09-18 08:20:02 UTC
diff --git a/MD5 b/MD5
@@ -1,8 +1,8 @@
-1c97db9bcfd27417ffa8389999775663 *DESCRIPTION
+2e1a1dc73c63f5072005124f7f6ee7ab *DESCRIPTION
 cda3ab64b68876145fdd6d6032379c16 *NAMESPACE
-e2693e137a7069f49ea543cb782ab69a *NEWS.md
+f9088b37832cd871673066795553b92b *NEWS.md
 6f005df720dee1e2788140dea3ac0d24 *R/WH-package.R
-661cd6119c047cab925bf37bd9ca2e15 *R/main.R
+cdbd6da5a263f78eef1b2b6e02912118 *R/main.R
 4e3a90501d0387bdb1035ebaf7c5501c *R/tools.R
 90cbe261c9ca614a2a18554cfc07b2b4 *README.md
 7c8bfb9d72b020e0c68f9873f61961de *build/vignette.rds
@@ -11,13 +11,13 @@ e2693e137a7069f49ea543cb782ab69a *NEWS.md
 13e5b5217d0ca0ae468daab576cdd2db *inst/CITATION
 d847721c40983fb15fdb62c0b9e141c8 *inst/doc/WH.R
 b2438c9365b57f8dc207cffc3b229c80 *inst/doc/WH.Rmd
-079c3f07a64175ceb26e3a4e3475977e *inst/doc/WH.html
+efba59e111c312606d010e0bd61abbf6 *inst/doc/WH.html
 2d7c65133d7934fd4a6f97fcdab7cc2e *man/WH-package.Rd
 d442149a14024490020343be20420da0 *man/WH_1d.Rd
 889d09be20652fdc6ffb4cc631ce0a5a *man/WH_1d_fixed_lambda.Rd
 6cac4e003cf20fe118935116c371f8cb *man/WH_1d_outer.Rd
 3ccec8ddf5f0423e803659ca4642ecef *man/WH_1d_perf.Rd
-706093609571493a51dd2c9a4d95c12f *man/WH_2d.Rd
+357d0d5c34cc66b4e095c0ca538f79d2 *man/WH_2d.Rd
 4ab1a2780f9ae7126464586d7d6577e7 *man/WH_2d_fixed_lambda.Rd
 59cbf4c1489a31144a94c1bac0297795 *man/WH_2d_outer.Rd
 4ab6d6eadc694865216157d2cd4d539e *man/WH_2d_perf.Rd
@@ -30,7 +30,7 @@ c9ba1845c022a3033eb0181914837eac *man/figures/README-plot-2.png
 322acdd22207c3b58d5b964ada6026d3 *man/figures/README-plot-3.png
 a97393938e4803c65f958748c3d3b49a *man/figures/README-plot-4.png
 7e48b945d48d76c2c88dd6bedd6138d5 *man/figures/README-plot-5.png
-0c196ca017663de0604926762c2e999a *man/figures/README-predict-1.png
+464d10356c7b350846c2b1a04227ce90 *man/figures/README-predict-1.png
 04fe2b0f7cfaa7d6f7d4d4a27a424c3f *man/figures/README-predict-2.png
 4e6e1418d443718c7de38fce6036de54 *man/get_diagnosis.Rd
 d5e82535ed8b3fceab4cc7509b156a4b *man/map.Rd
@@ -45,8 +45,8 @@ b82976558c28dc55d3cc2a4562b7faf9 *man/predict.WH_2d.Rd
 9a98438678dfc5b9ee8d219f97e5515e *man/print.WH_1d.Rd
 c4ed3089c4af14ec341cdcaa83e2841c *man/print.WH_2d.Rd
 828d0a66030e7522cf399e9e2c90bca4 *tests/testthat.R
-2e24e8df96029548c324f2efbbf747f2 *tests/testthat/test-1D.R
-13e9caf046004d85a5806f34d3dc0e47 *tests/testthat/test-2D.R
+98b0b2097af4a4ed1574e915654f2ed5 *tests/testthat/test-1D.R
+e430db3a04daae6211dac97fd26ca510 *tests/testthat/test-2D.R
 b2438c9365b57f8dc207cffc3b229c80 *vignettes/WH.Rmd
 ecfe2ac0176df273cadfa9e69fdb6815 *vignettes/Whittaker_Henderson_plots.R
 5f5187a8555f27efd51e5cfc43d0a846 *vignettes/biblio.bib
diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,29 @@
+# WH 1.1.0
+
+* Made significant changes to the extrapolation methods provided by be `predict.WH_1d` and `predict.WH_2d` functions. It turns out that the formula used for the extrapolation ignored the innovation error caused by the prior on the extrapolated region, resulting in smaller credibility intervals than they should have been. This now has been fixed.
+
+# WH 1.0.7
+
+* Removed tests of the form `expect_equal(f(x), f(x))` after confirmation they do not work with MKL BLAS
+
+# WH 1.0.6
+
+* Further improved test robustness
+
+* Added tests of the form `expect_equal(f(x), f(x))` for testing purposes
+
+# WH 1.0.5
+
+* Further improved tests robustness
+
+# WH 1.0.4
+
+* Simplified computation of the matrix `tUWU` by using the `crossprod` function, which should reduce memory usage (by half) and computation time (slightly)
+
+* Fixed an issue with the `WH_2d` plot and the what = "edf" argument
+
+* Improved tests robustness
+
 # WH 1.0.3
 
 * Replaced backquotes by normal quotes in the description field of the DESCRIPTION file

diff --git a/R/main.R b/R/main.R
@@ -182,7 +182,7 @@ WH_1d <- function(d, ec, lambda, criterion, method, q = 2, framework, y, wt, qui
 
 #' 2D Whittaker-Henderson Smoothing
 #'
-#' #' Main package function to apply Whittaker-Henderson smoothing in a
+#' Main package function to apply Whittaker-Henderson smoothing in a
 #' bidimensional survival analysis framework. It takes as input a matrix of
 #' observed events and a matrix of associated central exposure, both depending
 #' on two covariates, and build a smooth version of the log-hazard rate.
@@ -377,7 +377,8 @@ WH_2d <- function(d, ec, lambda, criterion, method, max_dim = 200, p,
   if (missing(p)) {
 
     max_ratio <- sqrt(max_dim / (n[[2]] * n[[1]]))
-    p <- floor(pmin(max_ratio, 1) * n)
+    p <- as.integer(pmin(max_ratio, 1) * n)
+    names(p) <- if (framework == "reg") names(dimnames(y)) else names(dimnames(d))
   } else {
     if (!is.numeric(q) || length(q) != 2 || any(q <= 0) ||
         (max(abs(q - round(q))) > .Machine$double.eps ^ 0.5)) stop(
@@ -439,32 +440,20 @@ predict.WH_1d <- function(object, newdata = NULL, ...) {
   full_data <- sort(union(data, newdata))
 
   n <- length(data)
-  n_pred <- length(full_data)
-  n_new <- n_pred - n
-
-  ind_fit <- which(full_data %in% data)
-  ind_inf <- which(full_data < min(data))
-  ind_sup <- which(full_data > max(data))
-  ind_new <- c(ind_inf, ind_sup)
-
-  D_mat_pred <- build_D_mat(n_pred, object$q)
+  n_inf <- sum(full_data < min(data))
+  n_sup <- sum(full_data > max(data))
+  n_tot <- n + n_inf + n_sup
 
-  D1_inf <- D_mat_pred[ind_inf, ind_fit, drop = FALSE]
-  D1_sup <- D_mat_pred[ind_sup - object$q, ind_fit, drop = FALSE]
-  D1 <- rbind(D1_inf, D1_sup)
+  ind_fit <- c(rep(FALSE, n_inf), rep(TRUE, n), rep(FALSE, n_sup)) |> which()
 
-  D2_inf <- D_mat_pred[ind_inf, ind_inf, drop = FALSE]
-  D2_sup <- D_mat_pred[ind_sup - object$q, ind_sup, drop = FALSE]
-  D2 <- blockdiag(D2_inf, D2_sup)
+  D_mat_pred <- build_D_mat(n_tot, object$q)
+  P_pred <- object$lambda * crossprod(D_mat_pred)
+  diag(P_pred[ind_fit, ind_fit]) <- diag(P_pred[ind_fit, ind_fit]) + object$wt
 
-  D2_inv <- if (n_new == 0) matrix(0, 0, 0) else solve(D2)
-
-  A_pred <- matrix(0, n_pred, n)
-  A_pred[ind_fit,] <- object$U
-  A_pred[ind_new,] <- - D2_inv %*% D1 %*% object$U
+  Psi_pred <- P_pred |> chol() |> chol2inv() # unconstrained variance / covariance matrix
 
-  y_pred <- c(A_pred %*% object$beta_hat)
-  std_y_pred <- sqrt(rowSums(A_pred * (A_pred %*% object$Psi)))
+  y_pred <- c(Psi_pred[,ind_fit] %*% (object$wt * object$z))
+  std_y_pred <- sqrt(diag(Psi_pred))
 
   names(y_pred) <- names(std_y_pred) <- full_data
 
@@ -509,33 +498,41 @@ predict.WH_2d <- function(object, newdata = NULL, ...) {
 
   data <- dimnames(object$y) |> lapply(as.numeric)
   full_data <- map2(data, newdata, \(x,y) sort(union(x, y)))
-  ind_fit <- map2(data, full_data, \(x,y) which(y %in% x))
 
   n <- map(data, length, "integer")
   n_inf <- map2(data, full_data, \(x,y) sum(y < min(x)), "integer")
   n_sup <- map2(data, full_data, \(x,y) sum(y > max(x)), "integer")
-  n_pred <- n + n_inf + n_sup
+  n_tot <- n + n_inf + n_sup
 
-  wt_pred <- matrix(0, n_pred[[1]], n_pred[[2]])
-  wt_pred[ind_fit[[1]], ind_fit[[2]]] <- object$wt
-
-  # W_pred <- diag(wt_pred) # extended weight matrix
-  D_mat_pred <- map2(n_pred, object$q, build_D_mat) # extended difference matrices
-  P_pred <- object$lambda[[1]] * diag(n_pred[[2]]) %x% crossprod(D_mat_pred[[1]]) +
-    object$lambda[[2]] * crossprod(D_mat_pred[[2]]) %x% diag(n_pred[[1]]) # extended penalization matrix
-  diag(P_pred) <- diag(P_pred) + c(wt_pred)
-  Psi_pred <- P_pred |> chol() |> chol2inv() # unconstrained variance / covariance matrix
+  prod_n <- prod(n)
+  prod_n_tot <- prod(n_inf + n + n_sup)
 
   ind_rows <- c(rep(FALSE, n_inf[[1]]), rep(TRUE, n[[1]]), rep(FALSE, n_sup[[1]]))
-  ind_coef_2d <- c(rep(FALSE, n_pred[[1]] * n_inf[[2]]), rep(ind_rows, n[[2]])) |> which()
+  ind_fit <- c(rep(FALSE, n_tot[[1]] * n_inf[[2]]), rep(ind_rows, n[[2]])) |> which()
+
+  D_mat_pred <- map2(n_tot, object$q, build_D_mat) # extended difference matrices
+  P_pred <- object$lambda[[1]] * diag(n_tot[[2]]) %x% crossprod(D_mat_pred[[1]]) +
+    object$lambda[[2]] * crossprod(D_mat_pred[[2]]) %x% diag(n_tot[[1]]) # extended penalization matrix
+  diag(P_pred[ind_fit, ind_fit]) <- diag(P_pred[ind_fit, ind_fit]) + c(object$wt)
+
+  Psi <- object$U %*% object$Psi %*% t(object$U)
+
+  P_12 <- P_pred[ind_fit, - ind_fit]
+  P_22_m <- P_pred[- ind_fit, - ind_fit] |> chol() |> chol2inv()
+
+  P_aux <- P_12 %*% P_22_m
+  P_aux_2 <- Psi %*% P_aux
 
-  Psi_inv <- Psi_pred[ind_coef_2d, ind_coef_2d] |> chol() |> chol2inv()
-  A_pred <- Psi_pred[,ind_coef_2d] %*% Psi_inv
+  Psi_pred <- matrix(0, prod_n_tot, prod_n_tot)
+  Psi_pred[ind_fit, ind_fit] <- Psi
+  Psi_pred[ind_fit, - ind_fit] <- - P_aux_2
+  Psi_pred[- ind_fit, ind_fit] <- t(Psi_pred[ind_fit, - ind_fit])
+  Psi_pred[- ind_fit, - ind_fit] <- t(P_aux) %*% P_aux_2 + P_22_m
 
-  y_pred <- c(A_pred %*% c(object$y_hat))
-  std_y_pred <- sqrt(rowSums(A_pred * (A_pred %*% (object$U %*% object$Psi %*% t(object$U)))))
+  y_pred <- c(Psi_pred[,ind_fit] %*% c(object$wt * object$z))
+  std_y_pred <- sqrt(diag(Psi_pred))
 
-  dim(y_pred) <- dim(std_y_pred) <- n_pred # set dimension for output matrices
+  dim(y_pred) <- dim(std_y_pred) <- n_tot # set dimension for output matrices
   dimnames(y_pred) <- dimnames(std_y_pred) <- full_data # set names for output matrices
 
   object$y_pred <- y_pred
@@ -808,7 +805,8 @@ plot.WH_2d <- function(x, what = "y_hat", trans, ...) {
   x <- unique(df$x)
   t <- unique(df$t)
 
-  data <- matrix(c(df[[what]]), length(t), length(x), byrow = TRUE)
+  data <- matrix(c(if (what == "edf") df[["edf_obs"]] else df[[what]]),
+                 length(t), length(x), byrow = TRUE)
 
   graphics::contour(
     t, x, data,
@@ -862,7 +860,7 @@ WH_1d_fixed_lambda <- function(d, ec, y, wt, lambda = 1e3, q = 2, p,
     z <- y
     wt_pos <- c(wt)[which_pos]
     z_pos <- c(z)[which_pos]
-    tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+    tUWU <- crossprod(sqrt(wt_pos) * U_pos)
     tUWz <- t(U_pos) %*% (wt_pos * z_pos)
 
   } else {
@@ -888,7 +886,7 @@ WH_1d_fixed_lambda <- function(d, ec, y, wt, lambda = 1e3, q = 2, p,
       z <- y_hat + d / wt - 1
       wt_pos <- c(wt)[which_pos]
       z_pos <- c(z)[which_pos]
-      tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+      tUWU <- crossprod(sqrt(wt_pos) * U_pos)
       tUWz <- t(U_pos) %*% (wt_pos * z_pos)
     }
 
@@ -929,7 +927,7 @@ WH_1d_fixed_lambda <- function(d, ec, y, wt, lambda = 1e3, q = 2, p,
   names(y_hat) <- names(std_y_hat) <- names(res) <- names(edf_obs) <-
     names(wt) <- names(z) <- names(y) # set names for output vectors
 
-  out <- list(y = y, wt = wt, y_hat = y_hat, std_y_hat = std_y_hat, res = res, edf_obs = edf_obs,
+  out <- list(y = y, wt = wt, z = z, y_hat = y_hat, std_y_hat = std_y_hat, res = res, edf_obs = edf_obs,
               beta_hat = beta_hat, edf_par = edf_par, diagnosis = diagnosis,
               U = U, Psi = Psi, lambda = lambda, p = p, q = q)
   class(out) <- "WH_1d"
@@ -972,7 +970,7 @@ WH_1d_outer <- function(d, ec, y, wt, q = 2, p, criterion = "REML", lambda = 1e3
     z <- y
     wt_pos <- c(wt)[which_pos]
     z_pos <- c(z)[which_pos]
-    tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+    tUWU <- crossprod(sqrt(wt_pos) * U_pos)
     tUWz <- t(U_pos) %*% (wt_pos * z_pos)
 
   } else {
@@ -1008,7 +1006,7 @@ WH_1d_outer <- function(d, ec, y, wt, q = 2, p, criterion = "REML", lambda = 1e3
         z <- y_hat + d / wt - 1
         wt_pos <- c(wt)[which_pos]
         z_pos <- c(z)[which_pos]
-        tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+        tUWU <- crossprod(sqrt(wt_pos) * U_pos)
         tUWz <- t(U_pos) %*% (wt_pos * z_pos)
       }
 
@@ -1086,7 +1084,7 @@ WH_1d_perf <- function(d, ec, y, wt, q = 2, p, criterion = "REML", lambda = 1e3,
     z <- y
     wt_pos <- c(wt)[which_pos]
     z_pos <- c(z)[which_pos]
-    tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+    tUWU <- crossprod(sqrt(wt_pos) * U_pos)
     tUWz <- t(U_pos) %*% (wt_pos * z_pos)
 
   } else {
@@ -1112,7 +1110,7 @@ WH_1d_perf <- function(d, ec, y, wt, q = 2, p, criterion = "REML", lambda = 1e3,
       z <- y_hat + d / wt - 1
       wt_pos <- c(wt)[which_pos]
       z_pos <- c(z)[which_pos]
-      tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+      tUWU <- crossprod(sqrt(wt_pos) * U_pos)
       tUWz <- t(U_pos) %*% (wt_pos * z_pos)
     }
 
@@ -1216,7 +1214,7 @@ WH_2d_fixed_lambda <- function(d, ec, y, wt, lambda = c(1e3, 1e3), q = c(2, 2),
     z <- y
     wt_pos <- c(wt)[which_pos]
     z_pos <- c(z)[which_pos]
-    tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+    tUWU <- crossprod(sqrt(wt_pos) * U_pos)
     tUWz <- t(U_pos) %*% (wt_pos * z_pos)
 
   } else {
@@ -1245,7 +1243,7 @@ WH_2d_fixed_lambda <- function(d, ec, y, wt, lambda = c(1e3, 1e3), q = c(2, 2),
       z <- y_hat + d / wt - 1
       wt_pos <- c(wt)[which_pos]
       z_pos <- c(z)[which_pos]
-      tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+      tUWU <- crossprod(sqrt(wt_pos) * U_pos)
       tUWz <- t(U_pos) %*% (wt_pos * z_pos)
     }
 
@@ -1289,7 +1287,7 @@ WH_2d_fixed_lambda <- function(d, ec, y, wt, lambda = c(1e3, 1e3), q = c(2, 2),
   dimnames(y_hat) <- dimnames(std_y_hat) <- dimnames(res) <- dimnames(edf_obs) <-
     dimnames(wt) <- dimnames(y) # set names for output matrices
 
-  out <- list(y = y, wt = wt, y_hat = y_hat, std_y_hat = std_y_hat, res = res, edf_obs = edf_obs,
+  out <- list(y = y, wt = wt, z = z, y_hat = y_hat, std_y_hat = std_y_hat, res = res, edf_obs = edf_obs,
               beta_hat = beta_hat, edf_par = edf_par, diagnosis = diagnosis,
               U = U, Psi = Psi, lambda = lambda, p = p, q = q)
   class(out) <- "WH_2d"
@@ -1327,7 +1325,7 @@ WH_2d_outer <- function(d, ec, y, wt, q = c(2, 2), p, criterion = "REML", lambda
     z <- y
     wt_pos <- c(wt)[which_pos]
     z_pos <- c(z)[which_pos]
-    tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+    tUWU <- crossprod(sqrt(wt_pos) * U_pos)
     tUWz <- t(U_pos) %*% (wt_pos * z_pos)
 
   } else {
@@ -1365,7 +1363,7 @@ WH_2d_outer <- function(d, ec, y, wt, q = c(2, 2), p, criterion = "REML", lambda
         z <- y_hat + d / wt - 1
         wt_pos <- c(wt)[which_pos]
         z_pos <- c(z)[which_pos]
-        tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+        tUWU <- crossprod(sqrt(wt_pos) * U_pos)
         tUWz <- t(U_pos) %*% (wt_pos * z_pos)
       }
 
@@ -1446,7 +1444,7 @@ WH_2d_perf <- function(d, ec, y, wt, q = c(2, 2), p, criterion = "REML", lambda
     z <- y
     wt_pos <- c(wt)[which_pos]
     z_pos <- c(z)[which_pos]
-    tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+    tUWU <- crossprod(sqrt(wt_pos) * U_pos)
     tUWz <- t(U_pos) %*% (wt_pos * z_pos)
 
   } else {
@@ -1472,7 +1470,7 @@ WH_2d_perf <- function(d, ec, y, wt, q = c(2, 2), p, criterion = "REML", lambda
       z <- y_hat + d / wt - 1
       wt_pos <- c(wt)[which_pos]
       z_pos <- c(z)[which_pos]
-      tUWU <- t(U_pos) %*% (wt_pos * U_pos)
+      tUWU <- crossprod(sqrt(wt_pos) * U_pos)
       tUWz <- t(U_pos) %*% (wt_pos * z_pos)
     }