Basic visualisation of inner and outer nested smooths DONE

NAMESPACE

@@ -24,6 +24,7 @@ S3method(plot,multi.ptermInteraction)
 S3method(plot,multi.ptermLogical)
 S3method(plot,multi.ptermNumeric)
 S3method(plot,multi.random.effect)
+S3method(plot,nested1D)
 S3method(plot,ptermFactor)
 S3method(plot,ptermInteraction)
 S3method(plot,ptermLogical)
@@ -104,6 +105,7 @@ export(plot.multi.ptermFactor)
 export(plot.multi.ptermLogical)
 export(plot.multi.ptermNumeric)
 export(plot.multi.random.effect)
+export(plot.nested1D)
 export(plot.ptermFactor)
 export(plot.ptermInteraction)
 export(plot.ptermLogical)

R/I_prepareInnerNested.R

@@ -2,37 +2,94 @@
 ##########
 # Internal method
 #
-.prepareInnerNested <- function(o, n, xlim, ...){
+.prepareInnerNested <- function(o, n, xlim, ylim = NULL, ...){
 
   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]
+  alpha <- si$alpha 
+  B <- si$B
+
+  da <- length( alpha )
+  prange <- (sm$first.para:sm$last.para)[1:da]
+  Va <- gObj$Vp[prange, prange, drop = FALSE]
 
   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)))
-
+    alpha <- drop(B %*% alpha)
+    Va <- B %*% Va %*% t(B)
+    se <- sqrt(pmax(0, diag(Va)))
     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)
+    out <- list("fit" = alpha, "x" = 1:da, "se" = se,
+                xlab = xlabel, ylab = ylabel, main = NULL, type = "si")
   }
   if( type == "nexpsm" ){
-    raw <- exp(alpha[1]) * (expsmooth(y = si$x, Xi = si$X, beta = alpha[-1])$d0 - si$xm)
-    trnam <- "expsm"
+    inner <- expsmooth(y = si$x, Xi = si$X, beta = alpha[-1], deriv = 1)
+    fit <- inner$d0
+    Jac <- inner$d1
+    nobs <- length(fit)
+    se <- sqrt(pmax(0, rowSums((Jac %*% Va[-1, -1, drop = FALSE]) * Jac)))
+    edf   <- sum(gObj$edf[prange[-1]])
+    ylabel <- .subEDF(paste0("expsm(", sm$term, ")"), edf)
+    xlabel <- "Index"
+    if( !is.null(xlim) ) {
+      xlim <- sort(xlim)
+      xlim[1] <- max(xlim[1], 1)
+      xlim[2] <- min(xlim[2], nobs)
+      ii <- which(1:nobs >= xlim[1] & 1:nobs <= xlim[2])
+      nobs <- length(ii)
+    } else {
+      xlim <- c(1, nobs)
+      ii <- 1:nobs
+    }
+    out <- list("fit" = fit[ii], "x" = ii, "se" = se[ii],
+                "p.resid" = si$x[ii], "raw" = ii, 
+                "xlim" = xlim, 
+                xlab = xlabel, ylab = ylabel, main = NULL, type = "nexpsm")
   }
+  if( type == "mgks" ){
+    d <- ncol(si$X0)
+    if( d != 2 ){ # ONLY 2D case handled at the moment!!
+      return( NULL )
+    }
+    if( !is.null(xlim) ) {
+      xlim <- sort(xlim)
+    } else {
+      xlim <- range(si$X[ , 1])
+    }
+    if( !is.null(ylim) ) {
+      ylim <- sort(ylim)
+    } else {
+      ylim <- range(si$X[ , 2])
+    }
+
+    xx <- rep(seq(xlim[1], xlim[2], length.out = n), n)
+    yy <- rep(seq(ylim[1], ylim[2], length.out = n), rep(n, n))
+
+    X <- cbind(xx, yy)
+    si$x <- as.matrix(si$x)
+    if( ncol(si$x) > 1 ){
+      si$x <- colMeans(si$x)
+    }
+
+    inner <- mgks(y = si$x, X = X, X0 = si$X0, beta = alpha[-1], deriv = 1)
+    fit <- inner$d0
+    Jac <- inner$d1
+    se <- sqrt(pmax(0, rowSums((Jac %*% Va[-1, -1, drop = FALSE]) * Jac)))
+    edf   <- sum(gObj$edf[prange[-1]])
+
+    mainlab <- .subEDF(paste0("mgks(", sm$term, ")"), edf)
+    ylabel <- "X[ , 2]"
+    xlabel <- "X[ , 1]"
+    out <- list("fit" = fit, "X" = si$X, "se" = se, "x" = xx, "y" = yy,
+                "p.resid" = si$x, "X0" = si$X0, 
+                "xlim" = xlim, "ylim" = ylim,
+                "xlab" = xlabel, "ylab" = ylabel, "main" = mainlab, type = "mgks")
+  }
+
+  return(out)
 
 }

R/I_prepareOuterNested.R

@@ -10,17 +10,20 @@
   alpha <- si$alpha
   dsi <- length( alpha )
 
+  rescale <- function(x){ exp(alpha[1]) * (x - si$xm) }
+
   type <- class(o)[1]
   if( type == "si" ){
    raw <- sort( si$X %*% alpha )
+   rescale <- function(x) x # No rescaling needed!
    trnam <- "proj"
   } 
   if( type == "nexpsm" ){
-   raw <- exp(alpha[1]) * (expsmooth(y = si$x, Xi = si$X, beta = alpha[-1])$d0 - si$xm)
+   raw <- expsmooth(y = si$x, Xi = si$X, beta = alpha[-1])$d0
    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)
+    raw <- mgks(y = si$x, X = si$X, X0 = si$X0, beta = alpha[-1])$d0
     trnam <- "mgks"
   }
 
@@ -35,7 +38,7 @@
   xx <- seq(xlim[1], xlim[2], length = n) 
 
   # Compute outer model matrix
-  X <- sm$xt$basis$evalX(x = xx, deriv = 0)$X0
+  X <- sm$xt$basis$evalX(x = rescale(xx), deriv = 0)$X0
 
   fit <- X %*% beta
 

R/L_ciBar.R

@@ -28,7 +28,7 @@ l_ciBar <- function(level = 0.95, mul = NULL, ...){
 #' @noRd
 #'
 l_ciBar.PtermFactor <- l_ciBar.MultiPtermNumeric <- l_ciBar.MultiPtermFactor <- 
-                       l_ciBar.ALE1DFactor <- l_ciBar.singleIndexInnerFactor <- function(a){
+                       l_ciBar.ALE1DFactor <- l_ciBar.siFactor <- function(a){
 
   xtra <- a$xtra
   a$xtra <- NULL

R/L_fitBar.R

@@ -26,7 +26,7 @@ l_fitBar <- function(a.aes = list(), ...){
 #' @noRd
 #'
 l_fitBar.PtermFactor <- l_fitBar.MultiPtermNumeric <- l_fitBar.MultiPtermFactor <- 
-                        l_fitBar.ALE1DFactor <- l_fitBar.singleIndexInnerFactor <- function(a){
+                        l_fitBar.ALE1DFactor <- l_fitBar.siFactor <- function(a){
 
   a$data <- a$data$fit
   if( is.null(a$na.rm) ){ a$na.rm <- TRUE}

R/L_fitContour.R

@@ -20,7 +20,7 @@ l_fitContour <- function(...){
 ######## Internal method 
 #' @noRd
 l_fitContour.2D <- l_fitContour.sos1 <- 
-l_fitContour.sos0 <- l_fitContour.MDslice <- function(a){
+l_fitContour.sos0 <- l_fitContour.MDslice <- l_fitContour.mgks2D <- function(a){
 
   a$data <- a$data$fit
   if( is.null(a$mapping) ) { a$mapping <- aes(z = tz) }

R/L_fitPoints.R

@@ -21,7 +21,7 @@ l_fitPoints <- function(...){
 #' @noRd
 #'
 l_fitPoints.PtermFactor <- l_fitPoints.MultiPtermNumeric <- 
-                           l_fitPoints.ALE1DFactor <- l_fitPoints.singleIndexInnerFactor <- function(a){
+                           l_fitPoints.ALE1DFactor <- l_fitPoints.siFactor <- function(a){
 
   if( is.null(a$shape) ){ a$shape <- 19}
   if( is.null(a$size) ){ a$size <- 2}

R/L_fitRaster.R

@@ -48,7 +48,7 @@ l_fitRaster.sos0 <- function(a){
 
 ######## Internal method 
 #' @noRd
-l_fitRaster.2D <- l_fitRaster.sos1 <- l_fitRaster.MDslice <- function(a){
+l_fitRaster.2D <- l_fitRaster.sos1 <- l_fitRaster.MDslice <- l_fitRaster.mgks2D <- function(a){
 
   xtra <- a$xtra
   a$xtra <- NULL

R/L_points.R

@@ -62,7 +62,7 @@ l_points.Check2DFactorFactor <- function(a){
 
 ######## Internal method for numeric/numeric 2D plots
 #' @noRd
-l_points.2D <- l_points.Check2DNumericNumeric <- l_points.MDslice <- function(a){
+l_points.2D <- l_points.Check2DNumericNumeric <- l_points.MDslice <- l_points.mgks2D <- function(a){
 
   return( l_points.1D(a) )
 

R/L_rug.R

@@ -94,7 +94,7 @@ l_rug.Check2DFactorFactor <- function(a){
 
 ######## Internal method for numeric/numeric 2D plots
 #' @noRd
-l_rug.2D <- l_rug.sos0 <- l_rug.sos1 <- l_rug.Check2DNumericNumeric <- l_rug.MDslice <- function(a){
+l_rug.2D <- l_rug.sos0 <- l_rug.sos1 <- l_rug.Check2DNumericNumeric <- l_rug.MDslice <- l_rug.mgks2D <- function(a){
 
   if( is.null(a$mapping) ) { a$mapping <- aes(x = x, y = y) }
 

R/plot_nested1D.R

@@ -4,6 +4,8 @@
 #' @description This method should be used to plot smooth effects 
 #'              of class \code{"si.smooth.1D"}.
 #' @param x a smooth effect object.
+#' @param inner if TRUE we are doing to plot the inner transformation, rather that then
+#'              outer smooth effect. 
 #' @param n number of grid points used to compute main effect and c.i. lines. 
 #'          For a nice smooth plot this needs to be several times the estimated degrees of 
 #'          freedom for the smooth.
@@ -20,13 +22,13 @@
 #' @export plot.nested1D
 #' @export
 #' 
-plot.nested1D <- function(x, n = 100, xlim = NULL, maxpo = 1e4, trans = identity, inner = FALSE, ...)  {
+plot.nested1D <- function(x, inner = FALSE, n = 100, xlim = NULL, ylim = NULL, maxpo = 1e4, trans = identity,  ...)  {
 
   if( inner ){
     # 1) Prepare data
-    P <- .prepareInnerNested(o = x, n = n, xlim = xlim, ...)
+    P <- .prepareInnerNested(o = x, n = n, xlim = xlim, ylim = ylim, ...)
 
-    out <- .plot.si.inner.smooth.1D(P = P, trans = trans)
+    out <- .plot.inner.nested.smooth.1D(P = P, trans = trans, maxpo = maxpo)
 
   } else {
     # 1) Prepare data
@@ -70,12 +72,55 @@ plot.nested1D <- function(x, n = 100, xlim = NULL, maxpo = 1e4, trans = identity
     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("nested", "1D")) )
+  return( list("ggObj" = .pl, "data" = .dat, type = "1D") )
 }
 
 ########################
 #' @noRd
-.plot.si.inner.smooth.1D <- function(P, trans) {
+.plot.inner.nested.smooth.1D <- function(P, trans, maxpo) {
+
+  if( is.null(P) ) { return(NULL) }
+
+  if(P$type == "nexpsm"){
+    return( .plot.mgcv.smooth.1D(x = NULL, P = P, trans = trans, maxpo = maxpo) )
+  }
+
+  if(P$type == "mgks"){
+    .dat <- list()
+    # 1) Build dataset on fitted effect
+    .dat$fit <- data.frame("z" = drop( P$fit ),
+                           "tz" = drop( trans(P$fit) ),
+                           "x" = rep(P$x, length(P$fit) / length(P$x)),
+                           "y" = rep(P$y, each = length(P$fit) / length(P$x)),
+                           "se" = P$se)
+
+    # 2) Build dataset on residuals
+    P$raw <- data.frame(z = P$p.resid, x = P$X0[ , 1], y = P$X0[ , 2])
+    if( !is.null(P$raw) ){
+
+      # Exclude residuals falling outside boundaries
+      .dat$res <- P$raw[P$raw$x >= P$xlim[1] & P$raw$x <= P$xlim[2] &
+                          P$raw$y >= P$ylim[1] & P$raw$y <= P$ylim[2] , , drop = FALSE]
+
+      # Sample if too many points (> maxpo)
+      nres <- nrow( .dat$res )
+      .dat$res$sub <- if(nres > maxpo) {
+        sample( c(rep(T, maxpo), rep(F, nres-maxpo)) )
+      } else {
+        rep(T, nres)
+      }
+    }
+
+    .dat$misc <- list("trans" = trans)
+
+    .pl <- ggplot(data = .dat$fit, aes(x = x, y = y, z = z)) +
+      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("mgks", "2D")) )
+
+  }
 
   .dat <- list()
   .dat$fit <- data.frame("x"  = as.factor(P$x),
@@ -86,11 +131,10 @@ plot.nested1D <- function(x, n = 100, xlim = NULL, maxpo = 1e4, trans = identity
 
   .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")), 
+
+  return( structure(list("ggObj" = .pl, "data" = .dat, "type" = c("si", "Factor")), 
                     class = c("plotSmooth",  "gg")) )
 
 }

R/print_plotGam.R

@@ -61,7 +61,7 @@ print.plotGam <- function(x, ask = TRUE, pages = NULL, addLay = TRUE, ...){
 
     .l <- switch(.cl, 
                  "nested1D" = .l + l_fitLine() + l_ciLine() + l_rug(), 
-                 "singleIndexInnerFactor" = .l + l_ciBar() + l_fitPoints(), 
+                 "siFactor" = .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(), 
@@ -82,6 +82,7 @@ print.plotGam <- function(x, ask = TRUE, pages = NULL, addLay = TRUE, ...){
                  "MultiPtermFactor" = .l + l_ciBar() + l_fitPoints(), 
                  "ALE1DNumeric" = .l + l_fitLine() + l_ciLine() + l_rug(),
                  "ALE1DFactor" = .l + l_ciBar() + l_fitPoints() + l_rug(),
+                 "mgks2D" = .l + l_fitRaster() + l_points(),
                   .l
                  )