Very mysterious debug error when running rstanarm/rstan chains — “Error in x$.self$finalize() : attempt to apply non-function”?

I am currently on:

  • Operating System: Linux, CentOS 7
  • rstanarm Version: rstanarm_2.17.4

I have been using rstanarm to fit bayesian models, and every now and then, I have two errors stack as shown below:

Error in (function (x)  : attempt to apply non-function
Error in (function (x)  : attempt to apply non-function

Doing a cursory check, more information I got was:

Error in x$.self$finalize() : attempt to apply non-function
Error in x$.self$finalize() : attempt to apply non-function
Error in (function (x)  : attempt to apply non-function

I found the exact error occurring in this link:

Does anyone have any ideas what might be causing this problem? The strangest thing is that there is some stochasticity to when the error appears. I can run the same model 5 times and it will output perhaps 2-3 times. I have also included a traceback below:

> traceback()
19: .External(list(name = "CppMethod__invoke_notvoid", address = <pointer: 0x5e53a10>,
        dll = list(name = "Rcpp", path = "/n/home07/user_admin/apps/R-3.5.0/Rcpp/libs/Rcpp.so",
            dynamicLookup = TRUE, handle = <pointer: 0x8d22c10>,
            info = <pointer: 0x1abcf80>), numParameters = -1L), <pointer: 0x1fb90a0>,
        <pointer: 0x8bde1a0>, .pointer, ...)
18: sampler$call_sampler(args_list[[i]])
17: doTryCatch(return(expr), name, parentenv, handler)
16: tryCatchOne(expr, names, parentenv, handlers[[1L]])
15: tryCatchList(expr, classes, parentenv, handlers)
14: tryCatch(expr, error = function(e) {
        call <- conditionCall(e)
        if (!is.null(call)) {
            if (identical(call[[1L]], quote(doTryCatch)))
                call <- sys.call(-4L)
            dcall <- deparse(call)[1L]
            prefix <- paste("Error in", dcall, ": ")
            LONG <- 75L
            sm <- strsplit(conditionMessage(e), "\n")[[1L]]
            w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w")
            if (is.na(w))
                w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L],
                    type = "b")
            if (w > LONG)
                prefix <- paste0(prefix, "\n  ")
        }
        else prefix <- "Error : "
        msg <- paste0(prefix, conditionMessage(e), "\n")
        .Internal(seterrmessage(msg[1L]))
        if (!silent && isTRUE(getOption("show.error.messages"))) {
            cat(msg, file = outFile)
            .Internal(printDeferredWarnings())
        }
        invisible(structure(msg, class = "try-error", condition = e))
    })
13: try(sampler$call_sampler(args_list[[i]]))
12: .local(object, ...)
11: (new("nonstandardGenericFunction", .Data = function (object,
        ...)
    {
        standardGeneric("sampling")
    }, generic = "sampling", package = "rstan", group = list(), valueClass = character(0),
        signature = "object", default = NULL, skeleton = (function (object,
            ...)
        stop("invalid call in method dispatch to 'sampling' (no default method)",
            domain = NA))(object, ...)))(object = new("stanmodel",
        model_name = "continuous", model_code = "#include /pre/Columbia_copyright.stan\n#include /pre/license.stan\n\n// GLM for a Gaussian, Gamma, inverse Gaussian, or Beta outcome\nfunctions {\n#include /functions/common_functions.stan\n#include /functions/        0
        "mean_PPD"), show_messages = FALSE, iter = 4000, warmup = 100,
        thin = 10, chains = 1, control = list(adapt_delta = 0.95,
            max_treedepth = 15L), save_warmup = FALSE)
10: (new("nonstandardGenericFunction", .Data = function (object,
        ...)
    {
        pars = c("alpha", "beta", "beta_smooth", "aux", "smooth_sd",
        "mean_PPD"), show_messages = FALSE, iter = 4000, warmup = 100,
        thin = 10, chains = 1, control = list(adapt_delta = 0.95,
            max_treedepth = 15L), save_warmup = FALSE)
9: do.call(sampling, sampling_args)
8: stan_glm.fit(x = X, y = y, weights = weights, offset = offset,
       family = family, prior = prior, prior_intercept = prior_intercept,
       prior_aux = prior_aux, prior_smooth = prior_smooth, prior_PD = prior_PD,
       algorithm = algorithm, adapt_delta = adapt_delta, group = group,
       QR = QR, ...)
7: (function (formula, random = NULL, family = gaussian(), data,
       weights = NULL, subset = NULL, na.action, knots = NULL, drop.unused.levels = TRUE,
       ..., prior = normal(), prior_intercept = normal(), prior_smooth = exponential(autoscale = FALSE),
       prior_aux = exponential(), prior_covariance = decov(), prior_PD = FALSE,
       algorithm = c("sampling", "meanfield", "fullrank"), adapt_delta = NULL,
       QR = FALSE, sparse = FALSE)
   {
       data <- validate_data(data, if_missing = list())
       family <- validate_family(family)
       if (!is.null(random)) {
           fake.formula <- as.character(mgcv::interpret.gam(formula)$fake.formula)
           form <- paste(fake.formula[2], fake.formula[1], fake.formula[3],
               "+", random[2], collapse = " ")
           glmod <- lme4::glFormula(as.formula(form), data, family = gaussian,
               subset, weights, na.action, control = make_glmerControl())
           data <- glmod$fr
           weights <- validate_weights(glmod$fr$weights)
       }
       else {
           weights <- validate_weights(weights)
           glmod <- NULL
       }
       if (family$family == "binomial") {
           data$temp_y <- rep(1, NROW(data))
           temp_formula <- update(formula, temp_y ~ .)
           jd <- mgcv::jagam(formula = temp_formula, family = gaussian(),
               data = data, file = tempfile(fileext = ".jags"),
               weights = NULL, na.action = na.action, offset = NULL,
               knots = knots, drop.unused.levels = drop.unused.levels,
               diagonalize = TRUE)
           if (!is.null(random)) {
               y <- data[, as.character(formula[2L])]
           }
           else {
               y <- eval(formula[[2L]], data)
           }
           if (binom_y_prop(y, family, weights)) {
               y1 <- as.integer(as.vector(y) * weights)
               y <- cbind(y1, y0 = weights - y1)
               weights <- double(0)
           }
       }
       else {
           jd <- mgcv::jagam(formula = formula, family = gaussian(),
               data = data, file = tempfile(fileext = ".jags"),
               weights = NULL, na.action = na.action, offset = NULL,
               knots = knots, drop.unused.levels = drop.unused.levels,
               diagonalize = TRUE)
           y <- jd$jags.data$y
       }
       offset <- validate_offset(as.vector(model.offset(jd$pregam$model)),
           y = y)
       X <- jd$jags.data$X
       mark <- which(colnames(X) != "")
       colnames(X) <- colnames(jd$pregam$X) <- jd$pregam$term.names
       S <- lapply(jd$pregam$smooth, FUN = function(s) {
           ranks <- s$rank
           start <- s$first.para
           out <- list()
           for (r in seq_along(ranks)) {
               end <- start + ranks[r] - 1L
               out[[r]] <- X[, start:end, drop = FALSE]
               start <- end + 1L
           }
           return(out)
       })
       if (any(sapply(S, length) > 1))
           S <- unlist(S, recursive = FALSE)
       names(S) <- names(jd$pregam$sp)
       X <- X[, mark, drop = FALSE]
       X <- c(list(X), S)
       if (is.null(prior))
           prior <- list()
       if (is.null(prior_intercept))
           prior_intercept <- list()
       if (is.null(prior_aux))
           prior_aux <- list()
       if (is.null(prior_smooth))
           prior_smooth <- list()
       if (is.null(random)) {
           group <- list()
           prior_covariance <- list()
       }
       else {
           group <- glmod$reTrms
           group$decov <- prior_covariance
       }
       algorithm <- match.arg(algorithm)
       stanfit <- stan_glm.fit(x = X, y = y, weights = weights,
           offset = offset, family = family, prior = prior, prior_intercept = prior_intercept,
           prior_aux = prior_aux, prior_smooth = prior_smooth, prior_PD = prior_PD,
           algorithm = algorithm, adapt_delta = adapt_delta, group = group,
           QR = QR, ...)
       if (family$family == "Beta regression")
           family$family <- "beta"
       X <- do.call(cbind, args = X)
       if (is.null(random))
           Z <- Matrix::Matrix(nrow = NROW(y), ncol = 0, sparse = TRUE)
       else {
           Z <- pad_reTrms(Ztlist = group$Ztlist, cnms = group$cnms,
               flist = group$flist)$Z
           colnames(Z) <- b_names(names(stanfit), value = TRUE)
       }
       XZ <- cbind(X, Z)
       mat <- as.matrix(stanfit)
       mark <- 1:ncol(X)
       jd$pregam$Vp <- cov(mat[, mark, drop = FALSE])
       jd$pregam$coefficients <- colMeans(mat[, mark, drop = FALSE])
       jd$pregam$sig2 <- if ("sigma" %in% colnames(mat))
           mean(mat[, "sigma"])
       else 1
       eta <- X %*% t(mat[, mark, drop = FALSE])
       mu <- rowMeans(family$linkinv(eta))
       eta <- rowMeans(eta)
       w <- as.numeric(jd$pregam$w * family$mu.eta(eta)^2/family$variance(mu))
       XWX <- t(X) %*% (w * X)
       jd$pregam$edf <- rowSums(jd$pregam$Vp * t(XWX))/jd$pregam$sig2
       class(jd$pregam) <- c("jam", "gam")
       fit <- nlist(stanfit, family, formula, offset, weights, x = XZ,
           y = y, data, terms = jd$pregam$terms, model = if (is.null(random))
               jd$pregam$model
           else glmod$fr, call = match.call(expand.dots = TRUE),
           algorithm, glmod = glmod, stan_function = "stan_gamm4")
       out <- stanreg(fit)
       out$jam <- jd$pregam
       class(out) <- c(class(out), "gamm4", if (!is.null(glmod)) "lmerMod")
       return(out)
   })(Y ~ Z + s(GG, k = 15) + s(gpsG, k = 30), random = NULL, data = list(0.5, 0.5, 0.5, 0.5, 1), warmup = 100, thin = 10, chains = 1)
6: do.call(stan_gamm4, args = list(formula, random = random, data = data,
       iter = iter, weights = weights, warmup = burn, thin = thin,
       chains = 1)) at Functions.R#144

Sounds like the stanfit or the stanmodel object is getting garbage collected but I have never seen that error before.

Thanks for your reply. This error has been very strange for me because it doesn’t shut down the process, nor does it seem to alter the results I have been getting in any theoretical way. When you mention garbage being collected, may I ask what part of the code indicates that? thanks!

Someone suggested that the “recompiling in order not to crash R” message was because the model object had been GC-ed. Is that not right?

GC-ed or otherwise not found via lexical search