Plotting of outer smooth of nested effects DONE

mfasiolo · Nov 24, 2023 · f064d52 · f064d52
1 parent 9098d6f
commit f064d52
Show file tree

Hide file tree

Showing 16 changed files with 114 additions and 63 deletions.
diff --git a/R/I_prepareInnerNested.R b/R/I_prepareInnerNested.R
@@ -0,0 +1,38 @@
+
+##########
+# Internal method
+#
+.prepareInnerNested <- function(o, n, xlim, ...){
+
+  gObj <- o$gObj
+  sm <- gObj$smooth[[ o$ism ]]
+
+  # Get single index vector
+  si <- sm$xt$si
+  dsi <- length( si$alpha )
+  prange <- (sm$first.para:sm$last.para)[1:dsi]
+
+  type <- class(o)[1]
+  if(type == "si"){
+    # Get parameters of inner transformation
+
+
+    fit <- si$B %*% si$alpha
+
+    xx <- 1:length(fit)
+
+    se <- sqrt(pmax(0, rowSums((si$B %*% gObj$Vp[prange, prange, drop = FALSE]) * si$B)))
+
+    edf   <- sum(gObj$edf[prange])
+    ylabel <- .subEDF(paste0("proj_coef(", sm$term, ")"), edf)
+    xlabel <- "Index"
+    out <- list("fit" = fit, "x" = xx, "se" = se,
+                xlab = xlabel, ylab = ylabel, main = NULL)
+    return(out)
+  }
+  if( type == "nexpsm" ){
+    raw <- exp(alpha[1]) * (expsmooth(y = si$x, Xi = si$X, beta = alpha[-1])$d0 - si$xm)
+    trnam <- "expsm"
+  }
+
+}
diff --git a/R/I_prepareInnerSI.R b/R/I_prepareInnerSI.R
diff --git a/R/I_prepareNested.R → R/I_prepareOuterNested.R b/R/I_prepareNested.R → R/I_prepareOuterNested.R
@@ -2,16 +2,28 @@
 ##########
 # Internal method that prepare non-linear effects plots (only single index at the moment)
 #
-.prepareNested<- function(o, n, xlim, ...){
+.prepareOuterNested <- function(o, n, xlim, ...){
 
   gObj <- o$gObj
   sm <- gObj$smooth[[ o$ism ]]
-
-  # Get single index vector
   si <- sm$xt$si
-  raw <- sort( si$X %*% si$alpha )
-  dsi <- length( si$alpha )
+  alpha <- si$alpha
+  dsi <- length( alpha )
 
+  type <- class(o)[1]
+  if( type == "si" ){
+   raw <- sort( si$X %*% alpha )
+   trnam <- "proj"
+  } 
+  if( type == "nexpsm" ){
+   raw <- exp(alpha[1]) * (expsmooth(y = si$x, Xi = si$X, beta = alpha[-1])$d0 - si$xm)
+   trnam <- "expsm"
+  }
+  if( type == "mgks" ){
+    raw <- exp(alpha[1]) * (mgks(y = si$x, X = si$X, X0 = si$X0, beta = alpha[-1])$d0 - si$xm)
+    trnam <- "mgks"
+  }
+
   # Get regression coeff of outer smooth
   prange <- (sm$first.para:sm$last.para)[-(1:dsi)]
   beta <- coef( gObj )[ prange ]
@@ -23,15 +35,15 @@
   xx <- seq(xlim[1], xlim[2], length = n) 
 
   # Compute outer model matrix
-  X <- sm$xt$splineDes(x = xx, deriv = 0)$X0
+  X <- sm$xt$basis$evalX(x = xx, deriv = 0)$X0
 
   fit <- X %*% beta
 
   se <- sqrt(pmax(0, rowSums((X %*% gObj$Vp[prange, prange, drop = FALSE]) * X)))
 
   edf   <- sum(gObj$edf[prange])
-  ylabel <- .subEDF(paste0("s(proj(", sm$term, "))"), edf)
-  xlabel <- paste0("proj(", sm$term, ")")
+  ylabel <- .subEDF(paste0("s(",trnam,"(", sm$term, "))"), edf)
+  xlabel <- paste0(trnam,"(", sm$term, ")")
   out <- list("fit" = fit, "x" = xx, "se" = se, "raw" = raw, "xlim" = xlim,
               xlab = xlabel, ylab = ylabel, main = NULL)
   return(out)

diff --git a/R/L_ciBar.R b/R/L_ciBar.R
@@ -27,7 +27,8 @@ l_ciBar <- function(level = 0.95, mul = NULL, ...){
 ######## Internal method for factor plots
 #' @noRd
 #'
-l_ciBar.PtermFactor <- l_ciBar.MultiPtermNumeric <- l_ciBar.MultiPtermFactor <- l_ciBar.ALE1DFactor <- function(a){
+l_ciBar.PtermFactor <- l_ciBar.MultiPtermNumeric <- l_ciBar.MultiPtermFactor <- 
+                       l_ciBar.ALE1DFactor <- l_ciBar.singleIndexInnerFactor <- function(a){
 
   xtra <- a$xtra
   a$xtra <- NULL

diff --git a/R/L_ciLine.R b/R/L_ciLine.R
@@ -27,7 +27,7 @@ l_ciLine <- function(level = 0.95, mul = NULL, ...){
 ######## Internal method 
 #' @noRd
 l_ciLine.1D <- l_ciLine.PtermNumeric <- l_ciLine.PtermMatrixNumeric <- 
-               l_ciLine.ALE1DNumeric <- l_ciLine.singleIndex1D <- function(a){
+               l_ciLine.ALE1DNumeric <- l_ciLine.nested1D <- function(a){
 
   xtra <- a$xtra
   a$xtra <- NULL

diff --git a/R/L_ciPoly.R b/R/L_ciPoly.R
@@ -26,7 +26,8 @@ l_ciPoly <- function(level = 0.95, mul = NULL, ...){
 
 ######## Internal method 
 #' @noRd
-l_ciPoly.1D <- l_ciPoly.PtermNumeric <- l_ciPoly.ALE1DNumeric <- function(a){
+l_ciPoly.1D <- l_ciPoly.PtermNumeric <- l_ciPoly.ALE1DNumeric <- 
+  l_ciPoly.nested1D <- function(a){
 
   xtra <- a$xtra
   a$xtra <- NULL

diff --git a/R/L_dens2D.R b/R/L_dens2D.R
@@ -50,7 +50,8 @@ l_dens <- l_dens2D
 
 ######## Internal method 
 #' @noRd
-l_dens2D.1D <- l_dens2D.Check1DNumeric <- l_dens2D.PtermNumeric <- function(a){
+l_dens2D.1D <- l_dens2D.Check1DNumeric <- l_dens2D.PtermNumeric <- 
+  l_dens2D.nested1D <- function(a){
 
   xtra <- a$xtra
   a$xtra <- NULL

diff --git a/R/L_fitBar.R b/R/L_fitBar.R
@@ -25,7 +25,8 @@ l_fitBar <- function(a.aes = list(), ...){
 ######## Internal method for factor plots
 #' @noRd
 #'
-l_fitBar.PtermFactor <- l_fitBar.MultiPtermNumeric <- l_fitBar.MultiPtermFactor <- l_fitBar.ALE1DFactor <- function(a){
+l_fitBar.PtermFactor <- l_fitBar.MultiPtermNumeric <- l_fitBar.MultiPtermFactor <- 
+                        l_fitBar.ALE1DFactor <- l_fitBar.singleIndexInnerFactor <- function(a){
 
   a$data <- a$data$fit
   if( is.null(a$na.rm) ){ a$na.rm <- TRUE}

diff --git a/R/L_fitLine.R b/R/L_fitLine.R
@@ -24,7 +24,7 @@ l_fitLine <- function(...){
 ######## Internal method 
 #' @noRd
 l_fitLine.1D <- l_fitLine.PtermNumeric <- l_fitLine.PtermMatrixNumeric <- 
-                l_fitLine.Multi1D <- l_fitLine.ALE1DNumeric <- l_fitLine.singleIndex1D <- function(a){
+                l_fitLine.Multi1D <- l_fitLine.ALE1DNumeric <- l_fitLine.nested1D <- function(a){
 
   a$data <- a$data$fit
   if( is.null(a$na.rm) ){ a$na.rm <- TRUE}

diff --git a/R/L_fitPoints.R b/R/L_fitPoints.R
@@ -20,7 +20,8 @@ l_fitPoints <- function(...){
 ######## Internal method for factor plots
 #' @noRd
 #'
-l_fitPoints.PtermFactor <- l_fitPoints.MultiPtermNumeric <- l_fitPoints.ALE1DFactor <- function(a){
+l_fitPoints.PtermFactor <- l_fitPoints.MultiPtermNumeric <- 
+                           l_fitPoints.ALE1DFactor <- l_fitPoints.singleIndexInnerFactor <- function(a){
 
   if( is.null(a$shape) ){ a$shape <- 19}
   if( is.null(a$size) ){ a$size <- 2}

diff --git a/R/L_points.R b/R/L_points.R
@@ -71,7 +71,8 @@ l_points.2D <- l_points.Check2DNumericNumeric <- l_points.MDslice <- function(a)
 ######## General internal method 
 #' @noRd
 l_points.1D <- l_points.sos1 <- l_points.sos0 <- 
-l_points.Check1DNumeric <- l_points.PtermNumeric <- function(a){
+l_points.Check1DNumeric <- l_points.PtermNumeric <-
+  l_points.nested1D <- function(a){
 
   a$data <- a$data$res[a$data$res$sub, ]
 

diff --git a/R/L_rug.R b/R/L_rug.R
@@ -59,7 +59,7 @@ l_rug.Check0DScalarNumeric <- function(a){
 ######## Internal method for numeric 1D plots
 #' @noRd
 l_rug.1D <- l_rug.Multi1D <- l_rug.Check1DNumeric <- 
-  l_rug.PtermNumeric <- l_rug.Check0DVectorNumeric <- l_rug.ALE1DNumeric <- l_rug.singleIndex1D <- function(a){
+  l_rug.PtermNumeric <- l_rug.Check0DVectorNumeric <- l_rug.ALE1DNumeric <- l_rug.nested1D <- function(a){
 
   if( is.null(a$mapping) ) { a$mapping <- aes(x = x) }
 

diff --git a/R/listLayers.R b/R/listLayers.R
@@ -44,7 +44,7 @@ listLayers <- function(o){
   lay <- lay[ startsWith(lay, "l_") & endsWith(lay, paste0(".", ty)) ]
 
   if( length(lay) == 0 ) {  
-    message(paste0("No layers available for objects of type", ty)) 
+    message(paste0("No layers available for objects of type ", ty)) 
     lay <- NULL
   } else {
     lay <- gsub(paste0(".", ty), '', lay)

diff --git a/R/plot_si_smooth_1D.R → R/plot_nested1D.R b/R/plot_si_smooth_1D.R → R/plot_nested1D.R
@@ -1,5 +1,5 @@
 #'
-#' Plotting one dimensional single index effects
+#' Plotting one dimensional nested effects
 #' 
 #' @description This method should be used to plot smooth effects 
 #'              of class \code{"si.smooth.1D"}.
@@ -15,23 +15,25 @@
 #'              Monotonicity is not checked. 
 #' @param ... currently unused.
 #' @return An object of class \code{c("plotSmooth", "gg")}.
-#' @name plot.si.smooth.1D
-#' @rdname plot.si.smooth.1D
-# @export plot.si.smooth.1D
-# @export
+#' @name plot.nested1D
+#' @rdname plot.nested1D
+#' @export plot.nested1D
+#' @export
 #' 
-plot.si.smooth.1D <- function(x, n = 100, xlim = NULL, maxpo = 1e4, trans = identity, inner = FALSE, ...)  {
+plot.nested1D <- function(x, n = 100, xlim = NULL, maxpo = 1e4, trans = identity, inner = FALSE, ...)  {
 
   if( inner ){
     # 1) Prepare data
-    P <- .prepareInnerSi(o = x, n = n, xlim = xlim, ...)
+    P <- .prepareInnerNested(o = x, n = n, xlim = xlim, ...)
+
+    out <- .plot.si.inner.smooth.1D(P = P, trans = trans)
 
   } else {
     # 1) Prepare data
-    P <- .prepareNested(o = x, n = n, xlim = xlim, ...)
+    P <- .prepareOuterNested(o = x, n = n, xlim = xlim, ...)
 
     # 2) Produce output object
-    out <- .plot.si.smooth.1D(x = P$smooth, P = P, trans = trans, maxpo = maxpo)
+    out <- .plot.outer.nested.1D(x = P$smooth, P = P, trans = trans, maxpo = maxpo)
   }
 
   class(out) <- c("plotSmooth", "gg")
@@ -41,14 +43,14 @@ plot.si.smooth.1D <- function(x, n = 100, xlim = NULL, maxpo = 1e4, trans = iden
 
 ########################
 #' @noRd
-.plot.si.smooth.1D <- function(x, P, trans, maxpo) {
+.plot.outer.nested.1D <- function(x, P, trans, maxpo) {
 
   .dat <- list()
 
   if ( !is.null(P$raw) ) {
     # Construct data.frame of partial residuals
     res <- data.frame("x" = as.vector(P$raw))
-
+    
     # Exclude residuals falling outside boundaries
     .dat$res <- res[res$x >= P$xlim[1] & res$x <= P$xlim[2], , drop = FALSE]
 
@@ -68,5 +70,27 @@ plot.si.smooth.1D <- function(x, n = 100, xlim = NULL, maxpo = 1e4, trans = iden
     labs(title = P$main, x = P$xlab, y = P$ylab) + theme_bw() + 
     theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
 
-  return( list("ggObj" = .pl, "data" = .dat, type = c("singleIndex", "1D")) )
+  return( list("ggObj" = .pl, "data" = .dat, type = c("nested", "1D")) )
+}
+
+########################
+#' @noRd
+.plot.si.inner.smooth.1D <- function(P, trans) {
+
+  .dat <- list()
+  .dat$fit <- data.frame("x"  = as.factor(P$x),
+                         "y"  = unname(P$fit),
+                         "ty"  = trans( unname(P$fit) ),
+                         "se" = unname(P$se) )
+  .dat$misc <- list("trans" = trans)
+
+  .pl <- ggplot(data = .dat$fit, aes("x" = x, "y" = ty)) +
+    labs(title = P$main, x = P$xlab, y = P$ylab) +
+    scale_x_discrete() +
+    theme_bw() +
+    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
+
+  return( structure(list("ggObj" = .pl, "data" = .dat, "type" = c("singleIndexInner", "Factor")), 
+                    class = c("plotSmooth",  "gg")) )
+
 }
diff --git a/R/print_plotGam.R b/R/print_plotGam.R
@@ -60,7 +60,8 @@ print.plotGam <- function(x, ask = TRUE, pages = NULL, addLay = TRUE, ...){
     } else { # [B] Smooth effect plots
 
     .l <- switch(.cl, 
-                 "singleIndex1D" = .l + l_fitLine() + l_ciLine() + l_rug(), 
+                 "nested1D" = .l + l_fitLine() + l_ciLine() + l_rug(), 
+                 "singleIndexInnerFactor" = .l + l_ciBar() + l_fitPoints(), 
                  "fs1D" = .l + l_fitLine() + theme(legend.position="none"),
                  "1D" = .l + l_fitLine() + l_ciLine() + l_rug(),
                  "2D" = .l + l_fitRaster() + l_fitContour(), 

diff --git a/R/sm.R b/R/sm.R
@@ -55,8 +55,8 @@ sm <- function(o, select){
     cl[which(cl=="fs.interaction")] <- paste("fs.interaction.", o$smooth[[select]]$dim, "D", sep='')
   }
 
-  if("si.smooth" %in% cl){
-    cl[which(cl=="si.smooth")] <- paste("si.smooth.", o$smooth[[select]]$dim, "D", sep='')
+  if("nested" %in% cl){
+    cl[which(cl=="nested")] <- paste("nested", o$smooth[[select]]$dim, "D", sep='')
   }
 
   class(out) <- cl