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


#1

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

#2

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


#3

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!