rversiongam

predict.gam behaviour gives error in R4 but works in R3. full reprex included


In R version 3.6.3, my predict.gam function gives a result. When I use R 4.x, I get errors for some types of predict, but not others.

My goal is to be able to use this code in R4.

library(gam)


gm1 <- predCol ~ s(VCF_nontree) + s(VCF_nonvegt) + s(VCF_treecov) + 
  s(CCI_URBAN) + s(GPW_POP)

gm2 <- predCol ~ offset(offset) + s(ACC) + s(EARS) + s(MAT) + s(PTA) + 
  s(SLP) + s(SOC) + s(TWI)


set.seed(42)
c2 <- data.frame(
  predCol = abs(rnorm(677072, 0.022, 0.03)),
  VCF_nontree = rnorm(677072, 6533, 16.3),
  VCF_nonvegt =rnorm(677072, 15.9, 12.03),
  VCF_treecov = rnorm(677072, 18.1, 18.3),
  CCI_URBAN = rnorm(677072, 0.00412, 0.013), 
  GPW_POP = rnorm(677072, 2.23, 1.28), 
  offset = rep(NA, 677072), 
  ACC = rnorm(677072, 612.6, 456.5), 
  EARS = rnorm(677072, 731.88, 291.4), 
  MAT = rnorm(677072, 22.69, 2.56), 
  PTA = rnorm(677072, 124.13, 307.9), 
  SLP = rnorm(677072, 2.41, 3.71), 
  SOC = rnorm(677072, 0.609, 0.33), 
  TWI = rnorm(677072, 118.2, 15.35)
)

fitModel1.1       <- gam::gam(gm1,
                              data = c2,
                              family = quasibinomial)
c2$offset  <- as.numeric(predict(fitModel1.1,
                                        type = "link"))

fitModel1         <- gam::gam(gm2,
                              data = c2,
                              family=quasibinomial)
fm1               <- predict(fitModel1,
                             se.fit = TRUE,
                             type = "response")

in R3 I get a result, in R4 I get this error:

Error in NextMethod("predict") : no method to invoke

R3, gam 1.20, R4 gam 1.22-2

in R4, this call to predict works:

fm2               <- predict(fitModel1,
                             se.fit = TRUE,
                             type = "link")

Solution

  • update the package maintainer has fixed the issue, and new version of gam on cran no longer need a workaround, and work as designed!

    Following the advice of @GavinSimpson, I made a workaround: a custom predict.gam() function that can be called in an Rscript using source("custom_predict_Gam.R") I save the below script in my working directory, then in R 4, use source(), then call the function like this:

    fm1               <- custom_predict_Gam(fitModel1,
                                            se.fit = TRUE,
                                            type = "response")
    

    and get results identical to those from R 3.

    Thanks for everyone's input. This has made me happy. gam::happy...

    On the off chance that this is useful to anyone else, here is the function that I save as custom_predict_Gam.R:

    
    # custom_predict_Gam.R
    
    # this file replaces the predict.gam on line 18 with predict. The reason is because predict.gam is broken in R 4.
    
    
    
    custom_predict_Gam <- function (object, newdata, type = c("link", "response", "terms"), 
                                    dispersion = NULL, se.fit = FALSE, na.action = na.pass, terms = labels(object), 
                                    ...) 
    {
      type <- match.arg(type)
      if (missing(newdata)) {
        if (inherits(object, "Gam") && !is.null(object$smooth)) {
          if (se.fit) 
            switch(type, 
                   response = {
                     out <- predict(object, type = "link", se.fit = TRUE, ...)
                     famob <- family(object)
                     out$se.fit <- drop(out$se.fit * abs(famob$mu.eta(out$fit)))
                     out$fit <- fitted(object)
                     out
                   }, 
                   link = {
                     out <- NextMethod("predict")
                     out$fit <- object$additive.predictors
                     TS <- out$residual.scale^2
                     TT <- ncol(object$var)
                     out$se.fit <- sqrt(out$se.fit^2 + TS * object$var %*% 
                                          rep(1, TT))
                     out
                   }, 
                   terms = {
                     out <- NextMethod("predict")
                     TT <- dimnames(s <- object$smooth)[[2]]
                     TT = intersect(terms, TT)
                     out$fit[, TT] <- out$fit[, TT] + s[, TT]
                     TS <- out$residual.scale^2
                     out$se.fit[, TT] <- sqrt(out$se.fit[, TT]^2 + 
                                                TS * object$var[, TT])
                     out
                   })
          else switch(type, 
                      terms = {
                        out <- NextMethod("predict")
                        TT <- dimnames(s <- object$smooth)[[2]]
                        TT = intersect(terms, TT)
                        out[, TT] <- out[, TT] + s[, TT]
                        out
                      }, 
                      link = object$additive.predictors, 
                      response = object$fitted)
        }
        else {
          if (inherits(object, "Gam")) {
            if (type == "link" && !se.fit) 
              object$additive.predictors
            else NextMethod("predict")
          }
          else UseMethod("predict")
        }
      }
      else newdata.predict.Gam(object, newdata, type, dispersion, 
                               se.fit, na.action, terms, ...)
    }