improvements to predict for non-gaussian models with varying u

Jouni Helske · Jouni Helske · commit 5bce041f2754 · 2025-05-19T08:55:43.000+03:00
diff --git a/ChangeLog b/ChangeLog
@@ -3,6 +3,12 @@ Changes from version 1.5.1 to 1.6.0
      users to abstract the functionality of `SSMcustom` via their own bespoke 
      function using, e.g., `SSMbespoke(myfun()) in the model formula. Thank you 
      Matteo Pelagatti for the suggestion.
+   * For non-Gaussian models, `predict` no longer checks component `u` when 
+     using n.ahead with `type=`link`.
+   * Added argument `u_new` to `predict`, which can be used together with 
+     `n.ahead` to define constant (series-specific) `u` for future predictions, 
+     thus avoiding the need for the `newdata` if model is otherwise 
+     time-invariant.
 
 Changes from version 1.5.0 to 1.5.1
    * Added explicit alias for KFAS-package due to changes in roxygen2.
diff --git a/R/predict.SSModel.R b/R/predict.SSModel.R
@@ -63,6 +63,9 @@
 #' of observations in the equations, leading to results which match with \code{glm} (where applicable).
 #' The latter case was the default behaviour of KFAS before version 1.3.8.
 #' Essentially this is the difference between observed and expected information in GLM context.
+#' @param u_new For non-Gaussian models, optional vector of length matching the 
+#' number of observation series. This defines the 'u' component to be used
+#' together with \code{n.ahead} argument.
 #' @param \dots Ignored.
 #' @return A matrix or list of matrices containing the predictions, and
 #'   optionally standard errors.
@@ -82,10 +85,12 @@
 #' pred <- predict(model,n.ahead=10,interval="prediction",level=0.9)
 #' pred
 predict.SSModel <- function(object, newdata, n.ahead,
-  interval = c("none", "confidence", "prediction"), level = 0.95,
-  type = c("response", "link"), states = NULL, se.fit = FALSE,  nsim = 0,
-  prob = TRUE, maxiter = 50, filtered = FALSE, expected = FALSE, ...) {
-
+                            interval = c("none", "confidence", "prediction"), 
+                            level = 0.95, type = c("response", "link"), 
+                            states = NULL, se.fit = FALSE, nsim = 0,
+                            prob = TRUE, maxiter = 50, filtered = FALSE, 
+                            expected = FALSE, u_new, ...) {
+  
   interval <- match.arg(interval)
   type <- match.arg(type)
   # Check that the model object is of proper form
@@ -102,8 +107,8 @@ predict.SSModel <- function(object, newdata, n.ahead,
         stop("Vector states should contain the indices or types of the states which are combined.")
     } else {
       states <- match.arg(arg = states, choices = c("all", "arima", "custom", "level","slope",
-        "cycle", "seasonal", "trend", "regression"),
-        several.ok = TRUE)
+                                                    "cycle", "seasonal", "trend", "regression"),
+                          several.ok = TRUE)
       if ("all" %in% states) {
         states <- as.integer(1:attr(object, "m"))
       } else {
@@ -131,50 +136,68 @@ predict.SSModel <- function(object, newdata, n.ahead,
     n <- attr(object, "n") <- no + nn
     timespan <- (no + 1):n
     object$y <- ts(rbind(object$y, newdata$y),
-      start = start(object$y), frequency = frequency(object$y))
+                   start = start(object$y), frequency = frequency(object$y))
     endtime <- end(object$y)
     tvo <- attr(object, "tv")
     tvn <- attr(newdata, "tv")
     same <- function(x, y) isTRUE(all.equal(x, y, tolerance = 0,
-      check.attributes = FALSE))
+                                            check.attributes = FALSE))
     if (tvo[1] || tvn[1] || !same(object$Z, newdata$Z)) {
       object$Z <- array(data = c(array(object$Z, dim = c(m, p, no)),
-        array(newdata$Z, dim = c(m, p, nn))), dim = c(p, m, n))
+                                 array(newdata$Z, dim = c(m, p, nn))), dim = c(p, m, n))
       attr(object, "tv")[1] <- 1L
     }
     if (gaussianmodel && (tvo[2] || tvn[2] || !same(object$H, newdata$H))) {
       object$H <- array(data = c(array(object$H, dim = c(p, p, no)),
-        array(newdata$H, dim = c(p, p, nn))), dim = c(p, p, n))
+                                 array(newdata$H, dim = c(p, p, nn))), dim = c(p, p, n))
       attr(object, "tv")[2] <- 1L
     } else if(!gaussianmodel) object$u <- rbind(object$u, matrix(newdata$u, nn, p))
-
+    
     if (tvo[3] || tvn[3] || !same(object$T, newdata$T)) {
       object$T <- array(data = c(array(object$T, dim = c(m, m, no)),
-        array(newdata$T, dim = c(m, m, nn))), dim = c(m, m, n))
+                                 array(newdata$T, dim = c(m, m, nn))), dim = c(m, m, n))
       attr(object, "tv")[3] <- 1L
     }
     if (tvo[4] || tvn[4] || !same(object$R, newdata$R)) {
       object$R <- array(data = c(array(object$R, dim = c(m, k, no)),
-        array(newdata$R, dim = c(m, k, nn))), dim = c(m, k, n))
+                                 array(newdata$R, dim = c(m, k, nn))), dim = c(m, k, n))
       attr(object, "tv")[4] <- 1L
     }
     if (tvo[5] || tvn[5] || !same(object$Q, newdata$Q)) {
       object$Q <- array(data = c(array(data = object$Q, dim = c(k, k, no)),
-        array(data = newdata$Q, dim = c(k, k, nn))), dim = c(k, k, n))
+                                 array(data = newdata$Q, dim = c(k, k, nn))), dim = c(k, k, n))
       attr(object, "tv")[5] <- 1L
     }
   } else {
     if (!missing(n.ahead) && !is.null(n.ahead)) {
       tv <- attr(object, "tv")
-      if(ifelse(gaussianmodel,any(tv), any(c(apply(object$u, 2, function(x) length(unique(x)) > 1))) || any(tv[-2])))
+      ok <- TRUE
+      if (gaussianmodel) {
+        ok <- !any(tv)
+      } else {
+        ok <- !any(tv[-2])
+        if (type == "response") {
+          varying_u <- !any(c(apply(object$u, 2, function(x) length(unique(x)) > 1)))
+          if (missing(u_new) && ok && !varying_u) {
+            stop("Component 'u' is time-varying. Either use 'newdata' instead of 'n.ahead', or use 'u_new' together with 'n.ahead'.")
+          }
+        }
+        if (!missing(u_new)) {
+          u_new <- rep(u_new, ncol(object$u))
+        } else {
+          u_new <- object$u[1L, ]
+        }
+      }
+      if (!ok) {
         stop("Model contains time varying system matrices, cannot use argument 'n.ahead'. Use 'newdata' instead.")
+      }
       timespan <- attr(object, "n") + 1:n.ahead
       n <- attr(object, "n") <- attr(object, "n") + as.integer(n.ahead)
       endtime<-end(object$y) + c(0, n.ahead)
       object$y <- window(object$y, end = endtime, extend = TRUE)
       if (any(object$distribution != "gaussian"))
-        object$u <- rbind(object$u, matrix(object$u[1, ], nrow = n.ahead,
-          ncol = ncol(object$u), byrow = TRUE))
+        object$u <- rbind(object$u, matrix(u_new, nrow = n.ahead,
+                                           ncol = ncol(object$u), byrow = TRUE))
     } else {
       timespan <- 1:attr(object, "n")
       endtime <- end(object$y)
@@ -193,7 +216,7 @@ predict.SSModel <- function(object, newdata, n.ahead,
         out <- KFS(model = object, filtering = "mean", smoothing = "none")
         names(out)[match(c("m", "P_mu"), names(out))]  <- c("muhat", "V_mu")
         if (out$d > 0) {
-        out$V_mu[,,1:out$d] <- Inf #diffuse phase
+          out$V_mu[,,1:out$d] <- Inf #diffuse phase
         }
       } else {
         out <- KFS(model = object, filtering = "none", smoothing = "mean")
@@ -205,25 +228,25 @@ predict.SSModel <- function(object, newdata, n.ahead,
         out <- signal(out,  states = states, filtered = TRUE)
         names(out) <- c("muhat", "V_mu")
         if (d > 0) {
-        out$V_mu[,,1:d] <- Inf #diffuse phase
+          out$V_mu[,,1:d] <- Inf #diffuse phase
         }
       } else {
         out <- signal(KFS(model = object, filtering = "none", smoothing = "state"),
-          states = states)
+                      states = states)
         names(out) <- c("muhat", "V_mu")
       }
-
+      
     }
     for (i in 1:p) {
       pred[[i]] <- cbind(fit = out$muhat[timespan, i],
-        switch(interval, none = NULL,
-          confidence = out$muhat[timespan, i] +
-            qnorm((1 - level)/2) * sqrt(out$V_mu[i, i, timespan]) %o% c(1, -1),
-          prediction = out$muhat[timespan, i] + qnorm((1 - level)/2) *
-            sqrt(out$V_mu[i, i, timespan] +
-                object$H[i, i, if (dim(object$H)[3] > 1) timespan else 1]) %o% c(1, -1)),
-        se.fit = if (se.fit)
-          sqrt(out$V_mu[i, i, timespan]))
+                         switch(interval, none = NULL,
+                                confidence = out$muhat[timespan, i] +
+                                  qnorm((1 - level)/2) * sqrt(out$V_mu[i, i, timespan]) %o% c(1, -1),
+                                prediction = out$muhat[timespan, i] + qnorm((1 - level)/2) *
+                                  sqrt(out$V_mu[i, i, timespan] +
+                                         object$H[i, i, if (dim(object$H)[3] > 1) timespan else 1]) %o% c(1, -1)),
+                         se.fit = if (se.fit)
+                           sqrt(out$V_mu[i, i, timespan]))
       if (interval != "none")
         colnames(pred[[i]])[2:3] <- c("lwr", "upr")
     }
@@ -232,41 +255,41 @@ predict.SSModel <- function(object, newdata, n.ahead,
       if (identical(states, as.integer(1:m))) {
         if (filtered) {
           out <- KFS(model = object, filtering = "signal", smoothing = "none",
-            maxiter = maxiter, expected = expected)
+                     maxiter = maxiter, expected = expected)
           names(out)[match(c("t", "P_theta"), names(out))] <- c("thetahat", "V_theta")
           if (out$d > 0) {
-          out$V_theta[,,1:out$d] <- Inf #diffuse phase
+            out$V_theta[,,1:out$d] <- Inf #diffuse phase
           }
         } else {
           out <- KFS(model = object, smoothing = "signal", maxiter = maxiter, expected = expected)
         }
       } else {
         if (filtered) {
           out <- KFS(model = object, filtering = "state", smoothing = "none",
-            maxiter = maxiter, expected = expected)
+                     maxiter = maxiter, expected = expected)
           d <- out$d
           out <- signal(out, states = states, filtered = TRUE)
           names(out) <- c("thetahat", "V_theta")
           if (d > 0) {
-          out$V_theta[,,1:d] <- Inf #diffuse phase
+            out$V_theta[,,1:d] <- Inf #diffuse phase
           }
         } else {
           out <- signal(KFS(model = object, smoothing = "state", 
-            maxiter = maxiter, expected = expected), states = states)
+                            maxiter = maxiter, expected = expected), states = states)
           names(out) <- c("thetahat", "V_theta")
         }
-
+        
       }
-
+      
       for (i in 1:p) {
         pred[[i]] <- cbind(fit = out$thetahat[timespan, i] +
-            (if (object$distribution[i] == "poisson")  log(object$u[timespan, i]) else 0),
-          switch(interval, none = NULL,
-            out$thetahat[timespan, i] +
-              (if (object$distribution[i] == "poisson") log(object$u[timespan, i]) else 0) +
-              qnorm((1 - level)/2) * sqrt(out$V_theta[i, i, timespan]) %o%
-              c(1, -1)), se.fit = if (se.fit)
-                sqrt(out$V_theta[i, i, timespan]))
+                             (if (object$distribution[i] == "poisson")  log(object$u[timespan, i]) else 0),
+                           switch(interval, none = NULL,
+                                  out$thetahat[timespan, i] +
+                                    (if (object$distribution[i] == "poisson") log(object$u[timespan, i]) else 0) +
+                                    qnorm((1 - level)/2) * sqrt(out$V_theta[i, i, timespan]) %o%
+                                    c(1, -1)), se.fit = if (se.fit)
+                                      sqrt(out$V_theta[i, i, timespan]))
         if (interval == "confidence")
           colnames(pred[[i]])[2:3] <- c("lwr", "upr")
       }
@@ -275,84 +298,84 @@ predict.SSModel <- function(object, newdata, n.ahead,
           tmp <- which(colnames(pred[[1]]) == "se.fit")
           for (i in 1:p) {
             pred[[i]][, "se.fit"] <- switch(object$distribution[i],
-              gaussian = pred[[i]][,"se.fit"],
-              poisson = pred[[i]][, "se.fit"] * exp(pred[[i]][, 1]),
-              binomial = pred[[i]][, "se.fit"] *
-                (if (!prob) object$u[timespan, i] else 1) *
-                exp(pred[[i]][, 1])/(1 + exp(pred[[i]][, 1]))^2,
-              gamma = pred[[i]][, "se.fit"] * exp(pred[[i]][, 1]),
-              `negative binomial` = pred[[i]][, "se.fit"] * exp(pred[[i]][, 1]))
+                                            gaussian = pred[[i]][,"se.fit"],
+                                            poisson = pred[[i]][, "se.fit"] * exp(pred[[i]][, 1]),
+                                            binomial = pred[[i]][, "se.fit"] *
+                                              (if (!prob) object$u[timespan, i] else 1) *
+                                              exp(pred[[i]][, 1])/(1 + exp(pred[[i]][, 1]))^2,
+                                            gamma = pred[[i]][, "se.fit"] * exp(pred[[i]][, 1]),
+                                            `negative binomial` = pred[[i]][, "se.fit"] * exp(pred[[i]][, 1]))
             pred[[i]][, -tmp] <- switch(object$distribution[i],
-              gaussian = pred[[i]][, -tmp],
-              poisson = exp(pred[[i]][, -tmp]),
-              binomial = (if (!prob) object$u[timespan, i] else 1) *
-                exp(pred[[i]][, -tmp])/(1 + exp(pred[[i]][, -tmp])),
-              gamma = exp(pred[[i]][, -tmp]),
-              `negative binomial` = exp(pred[[i]][, -tmp]))
+                                        gaussian = pred[[i]][, -tmp],
+                                        poisson = exp(pred[[i]][, -tmp]),
+                                        binomial = (if (!prob) object$u[timespan, i] else 1) *
+                                          exp(pred[[i]][, -tmp])/(1 + exp(pred[[i]][, -tmp])),
+                                        gamma = exp(pred[[i]][, -tmp]),
+                                        `negative binomial` = exp(pred[[i]][, -tmp]))
           }
         } else {
           for (i in 1:p) pred[[i]] <- switch(object$distribution[i],
-            gaussian = pred[[i]],
-            poisson = exp(pred[[i]]),
-            binomial = (if (!prob) object$u[timespan, i] else 1) *
-              exp(pred[[i]])/(1 + exp(pred[[i]])),
-            gamma = exp(pred[[i]]),
-            `negative binomial` = exp(pred[[i]]))
+                                             gaussian = pred[[i]],
+                                             poisson = exp(pred[[i]]),
+                                             binomial = (if (!prob) object$u[timespan, i] else 1) *
+                                               exp(pred[[i]])/(1 + exp(pred[[i]])),
+                                             gamma = exp(pred[[i]]),
+                                             `negative binomial` = exp(pred[[i]]))
         }
       }
-
+      
     } else {# with importance sampling
       if (filtered) {
         d <- KFS(approxSSM(object, maxiter = maxiter, expected = expected), smoothing = "none")$d
       }
-
+      
       if (interval == "none") {
         imp <- importanceSSM(object,
-          ifelse(identical(states, as.integer(1:m)), "signal", "states"),
-          nsim = nsim, antithetics = TRUE, maxiter = maxiter, filtered = filtered,
-          expected = expected)
+                             ifelse(identical(states, as.integer(1:m)), "signal", "states"),
+                             nsim = nsim, antithetics = TRUE, maxiter = maxiter, filtered = filtered,
+                             expected = expected)
         nsim <- as.integer(4 * nsim)
         if (!identical(states, as.integer(1:m))) {
           imp$samples <- .Fortran(fzalpha, NAOK = TRUE, as.integer(dim(object$Z)[3] > 1),
-            object$Z, imp$samples, signal = array(0, c(n, p, nsim)),
-            p, m, n, nsim, length(states), states)$signal
+                                  object$Z, imp$samples, signal = array(0, c(n, p, nsim)),
+                                  p, m, n, nsim, length(states), states)$signal
         }
-
+        
         w <- imp$weights/sum(imp$weights)
         if (type == "response") {
           for (i in 1:p) {
             imp$samples[timespan, i, ] <- switch(object$distribution[i],
-              gaussian = imp$samples[timespan, i, ],
-              poisson = object$u[timespan, i] * exp(imp$samples[timespan, i, ]),
-              binomial = (if (!prob) object$u[timespan, i] else 1) *
-                exp(imp$samples[timespan, i, ])/(1 + exp(imp$samples[timespan, i, ])),
-              gamma = exp(imp$samples[timespan, i, ]),
-              `negative binomial` = exp(imp$samples[timespan, i, ]))
+                                                 gaussian = imp$samples[timespan, i, ],
+                                                 poisson = object$u[timespan, i] * exp(imp$samples[timespan, i, ]),
+                                                 binomial = (if (!prob) object$u[timespan, i] else 1) *
+                                                   exp(imp$samples[timespan, i, ])/(1 + exp(imp$samples[timespan, i, ])),
+                                                 gamma = exp(imp$samples[timespan, i, ]),
+                                                 `negative binomial` = exp(imp$samples[timespan, i, ]))
           }
         } else {
           for (i in 1:p) if (object$distribution[i] == "poisson")
             imp$samples[timespan, i, ] <- imp$samples[timespan, i, ] + log(object$u[timespan,
-              i])
+                                                                                    i])
         }
         varmean <- .Fortran(fvarmeanw, NAOK = TRUE, imp$samples[timespan, , , drop = FALSE], w,
-          p, length(timespan),
-          nsim, mean = array(0, c(length(timespan), p)),
-          var = array(0, c(length(timespan), p)), as.integer(se.fit))
-
+                            p, length(timespan),
+                            nsim, mean = array(0, c(length(timespan), p)),
+                            var = array(0, c(length(timespan), p)), as.integer(se.fit))
+        
         if (se.fit) {
           if (filtered && d > 0) {
             varmean$var[1:d, ] <- Inf #diffuse phase
           }
           pred <- lapply(1:p, function(j) cbind(fit = varmean$mean[, j],
-            se.fit = sqrt(varmean$var[, j])))
-
+                                                se.fit = sqrt(varmean$var[, j])))
+          
         } else {
           pred <- lapply(1:p, function(j) varmean$mean[, j])
         }
       } else {
         pred <- interval(object, interval = interval, level = level, type = type,
-          states = states, nsim = nsim, se.fit = se.fit, timespan = timespan,
-          prob = prob, maxiter = maxiter, filtered = filtered, expected = expected)
+                         states = states, nsim = nsim, se.fit = se.fit, timespan = timespan,
+                         prob = prob, maxiter = maxiter, filtered = filtered, expected = expected)
         if (filtered && d > 0) {
           for (i in 1:p) {
             pred[[i]][1:d, "lwr"] <- -Inf
@@ -361,7 +384,7 @@ predict.SSModel <- function(object, newdata, n.ahead,
               pred[[i]][1:d, "se.fit"] <- Inf #diffuse phase
             }
           }
-
+          
         }
       }
     }
diff --git a/man/predict.SSModel.Rd b/man/predict.SSModel.Rd
