Running examples in ‘vagam-Ex.R’ failed
The error most likely occurred in:

> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: vagam
> ### Title: Fitting generalized additive models (GAMs) using variational
> ###   approximations (VA).
> ### Aliases: vagam
> ### Keywords: methods
> 
> ### ** Examples
> 
> ## Example 1: Application to wage data
> data(wage_data)
> 
> south_code <- gender_code <- race_code <- union_code <- vector("numeric", nrow(wage_data))
> union_code[wage_data$union == "member"] <- 1
> south_code[wage_data$south == "yes"] <- 1
> gender_code[wage_data$gender == "female"] <- 1
> race_code[wage_data$race == "White"] <- 1
> para.X <- data.frame(south = south_code, gender = gender_code, race = race_code)
> 
> fit_va <- vagam(y = union_code, smooth.X = wage_data[,c("education", "wage", "age")],
+                        para.X = para.X,
+                        int.knots = 8, save.data = TRUE, 
+                        family = binomial(), 
+                        para.se = TRUE)
Lambda updated as part of VA estimation. Yeah baby!
Iteration: 1 	 Current VA logL: -Inf  | New VA logL: -251.9948  | Difference: Inf 
Iteration: 2 	 Current VA logL: -251.9948  | New VA logL: -249.6027  | Difference: 2.392096 
Iteration: 3 	 Current VA logL: -249.6027  | New VA logL: -248.813  | Difference: 0.789726 
Iteration: 4 	 Current VA logL: -248.813  | New VA logL: -248.2817  | Difference: 0.5312606 
Iteration: 5 	 Current VA logL: -248.2817  | New VA logL: -247.9144  | Difference: 0.3673779 
Iteration: 6 	 Current VA logL: -247.9144  | New VA logL: -247.6557  | Difference: 0.2587032 
Iteration: 7 	 Current VA logL: -247.6557  | New VA logL: -247.4697  | Difference: 0.1859226 
Iteration: 8 	 Current VA logL: -247.4697  | New VA logL: -247.3333  | Difference: 0.1364035 
Iteration: 9 	 Current VA logL: -247.3333  | New VA logL: -247.2313  | Difference: 0.1020059 
Iteration: 10 	 Current VA logL: -247.2313  | New VA logL: -247.1537  | Difference: 0.07758634 
Iteration: 11 	 Current VA logL: -247.1537  | New VA logL: -247.0939  | Difference: 0.05988457 
Iteration: 12 	 Current VA logL: -247.0939  | New VA logL: -247.047  | Difference: 0.04680679 
Iteration: 13 	 Current VA logL: -247.047  | New VA logL: -247.0101  | Difference: 0.03698225 
Iteration: 14 	 Current VA logL: -247.0101  | New VA logL: -246.9806  | Difference: 0.02949385 
Iteration: 15 	 Current VA logL: -246.9806  | New VA logL: -246.9569  | Difference: 0.02371432 
Iteration: 16 	 Current VA logL: -246.9569  | New VA logL: -246.9377  | Difference: 0.01920535 
Iteration: 17 	 Current VA logL: -246.9377  | New VA logL: -246.922  | Difference: 0.01565459 
Iteration: 18 	 Current VA logL: -246.922  | New VA logL: -246.9092  | Difference: 0.01283546 
Iteration: 19 	 Current VA logL: -246.9092  | New VA logL: -246.8986  | Difference: 0.01058095 
Iteration: 20 	 Current VA logL: -246.8986  | New VA logL: -246.8898  | Difference: 0.008766283 
Iteration: 21 	 Current VA logL: -246.8898  | New VA logL: -246.8825  | Difference: 0.007297091 
Iteration: 22 	 Current VA logL: -246.8825  | New VA logL: -246.8764  | Difference: 0.006101228 
Iteration: 23 	 Current VA logL: -246.8764  | New VA logL: -246.8713  | Difference: 0.005123037 
Iteration: 24 	 Current VA logL: -246.8713  | New VA logL: -246.867  | Difference: 0.00431922 
Iteration: 25 	 Current VA logL: -246.867  | New VA logL: -246.8633  | Difference: 0.003655846 
Iteration: 26 	 Current VA logL: -246.8633  | New VA logL: -246.8602  | Difference: 0.003106149 
Iteration: 27 	 Current VA logL: -246.8602  | New VA logL: -246.8576  | Difference: 0.002648886 
Iteration: 28 	 Current VA logL: -246.8576  | New VA logL: -246.8553  | Difference: 0.002267108 
Iteration: 29 	 Current VA logL: -246.8553  | New VA logL: -246.8533  | Difference: 0.001947219 
Iteration: 30 	 Current VA logL: -246.8533  | New VA logL: -246.8517  | Difference: 0.001678265 
Iteration: 31 	 Current VA logL: -246.8517  | New VA logL: -246.8502  | Difference: 0.00145138 
Iteration: 32 	 Current VA logL: -246.8502  | New VA logL: -246.849  | Difference: 0.001259362 
Iteration: 33 	 Current VA logL: -246.849  | New VA logL: -246.8479  | Difference: 0.001096333 
Iteration: 34 	 Current VA logL: -246.8479  | New VA logL: -246.8469  | Difference: 0.0009574839 
Calculating information matrix for model parameters...
 ----------- FAILURE REPORT -------------- 
 --- failure: the condition has length > 1 ---
 --- srcref --- 
: 
 --- package (from environment) --- 
vagam
 --- call from context --- 
vagam(y = union_code, smooth.X = wage_data[, c("education", "wage", 
    "age")], para.X = para.X, int.knots = 8, save.data = TRUE, 
    family = binomial(), para.se = TRUE)
 --- call from argument --- 
while (class(try(solve(obs_info), silent = TRUE)) == "try-error") {
    obs_info <- info.valouis(y = y, para.X = para.X, Z = Z, kappa = new_kappa, 
        phi = new_phi, lambda = new_lambda, a = new_a, A = new_A, 
        S = S, d = d, index.cov = index_cov, mc.samps = control$mc.samps, 
        family = family)
}
 --- R stacktrace ---
where 1: vagam(y = union_code, smooth.X = wage_data[, c("education", "wage", 
    "age")], para.X = para.X, int.knots = 8, save.data = TRUE, 
    family = binomial(), para.se = TRUE)

 --- value of length: 2 type: logical ---
[1] FALSE FALSE
 --- function from context --- 
function (y, smooth.X, para.X = NULL, lambda = NULL, int.knots, 
    family = gaussian(), A.struct = c("unstructured", "block"), 
    offset = NULL, save.data = FALSE, para.se = FALSE, doIC = FALSE, 
    control = list(eps = 0.001, maxit = 1000, trace = TRUE, seed.number = 123, 
        mc.samps = 4000, pois.step.size = 0.01)) 
{
    if (!(family$family %in% c("gaussian", "poisson", "binomial"))) 
        stop("Specified family currently no supported. Sorry!")
    if (family$family == "gaussian") 
        family <- gaussian(link = "identity")
    if (family$family == "poisson") 
        family <- poisson(link = "log")
    if (family$family == "binomial") 
        family <- binomial(link = "logit")
    if (family$family == "binomialp") 
        family <- binomial(link = "probit")
    A.struct <- match.arg(A.struct, c("unstructured", "block"))
    ncov_smooth_X <- ncol(smooth.X)
    if (is.null(offset)) 
        offset <- numeric(length(y))
    if (!(length(int.knots) %in% c(1, ncov_smooth_X))) 
        stop("int.knots should either be a single number of a vector with length equal to ncol(smooth.X). Thanks!")
    if (length(int.knots) == 1) 
        actual_int_knots <- rep(int.knots, ncov_smooth_X)
    if (length(int.knots) > 1) 
        actual_int_knots <- int.knots
    if (is.null(colnames(smooth.X))) 
        colnames(smooth.X) <- paste("smoothX", 1:ncov_smooth_X, 
            sep = "")
    if (!is.null(para.X)) {
        para.X <- as.matrix(para.X)
        if (is.null(colnames(para.X))) 
            colnames(para.X) <- paste("paraX", 1:ncol(para.X), 
                sep = "")
        if (any(apply(para.X, 2, function(x) all(x == 1)))) 
            stop("No intercept terms should be included in para.X, as this is included by default. Thanks!")
        para.X <- cbind(1, para.X)
        colnames(para.X)[1] <- "Intercept"
    }
    if (is.null(para.X)) {
        para.X <- matrix(1, nrow = length(y), ncol = 1)
        colnames(para.X) <- "Intercept"
    }
    if (is.null(lambda)) {
        cw_lambda <- new_lambda <- rep(2, ncov_smooth_X)
        message("Lambda updated as part of VA estimation. Yeah baby!")
    }
    if (!is.null(lambda)) {
        if (length(lambda) != ncov_smooth_X) 
            stop("lambda should a vector with length equal to ncol(smooth.X). Thanks!")
        cw_lambda <- new_lambda <- lambda
        message("Lambda given. Thanks!")
    }
    S <- basis_info <- vector("list", ncov_smooth_X)
    d <- numeric(ncov_smooth_X)
    index_cov <- NULL
    Z <- NULL
    for (k in 1:ncov_smooth_X) {
        x0 <- smooth.X[, k]
        get_basis <- smooth.construct(s(x0, bs = "ps", k = actual_int_knots[k] + 
            3, m = c(2, 1)), data = data.frame(x0), knots = NULL)
        Z_k <- get_basis$X
        transform_mat <- qr.Q(qr(as.matrix(colSums(Z_k))), complete = TRUE)[, 
            -1]
        Z_k <- Z_k %*% transform_mat
        index_cov <- c(index_cov, rep(k, ncol(Z_k)))
        colnames(Z_k) <- paste("smooth.X", k, 1:sum(index_cov == 
            k), sep = "")
        Z <- cbind(Z, Z_k)
        S[[k]] <- crossprod(transform_mat, get_basis$S[[1]]) %*% 
            transform_mat
        d[k] <- get_basis$rank
        get_basis$X <- NULL
        get_basis$transform_mat <- transform_mat
        basis_info[[k]] <- get_basis
        rm(transform_mat, get_basis, x0)
    }
    cw_kappa <- new_kappa <- numeric(ncol(para.X))
    names(new_kappa) <- names(para.X)
    cw_a <- new_a <- numeric(ncol(Z))
    if (family$family == "poisson") {
        old <- .Random.seed
        on.exit({
            .Random.seed = old
        })
        set.seed(control$seed.number)
        cw_a <- new_a <- rnorm(ncol(Z))
    }
    names(new_a) <- colnames(Z)
    cw_A <- new_A <- matrix(0, nrow = ncol(Z), ncol = ncol(Z))
    rownames(new_A) <- colnames(new_A) <- colnames(Z)
    cw_phi <- new_phi <- 1
    cw_logL <- -Inf
    new_logL <- 10
    diff_logL <- 10
    diff_lambda <- 0
    counter <- 1
    if (family$family %in% c("gaussian", "binomialp")) 
        ZZt <- crossprod(Z)
    tic <- proc.time()
    while (diff_logL > control$eps) {
        if (counter > control$maxit) 
            break
        ZA <- Z %*% cw_A
        new_offset <- Z %*% cw_a + as.numeric(family$family %in% 
            c("poisson", "binomial")) * 0.5 * rowSums(ZA * Z) + 
            offset
        fit0 <- suppressWarnings(glm.fit(x = para.X, y = y, family = family, 
            offset = new_offset, intercept = FALSE))
        new_kappa <- fit0$coefficients
        Q <- matrix(0, ncol(Z), ncol(Z))
        for (k2 in 1:ncov_smooth_X) Q[index_cov == k2, index_cov == 
            k2] <- S[[k2]] * cw_lambda[k2]
        new_offset <- para.X %*% new_kappa + as.numeric(family$family %in% 
            c("poisson", "binomial")) * 0.5 * rowSums(ZA * Z) + 
            offset
        if (family$family == "poisson") {
            fit0 <- ridge.iwls(x = Z, y = y, family = family, 
                penalty = Q, initial.beta = cw_a, offset = new_offset, 
                gamma = 0.1 + min(counter * control$pois.step.size, 
                  0.9))
        }
        if (family$family != "poisson") {
            fit0 <- ridge.iwls(x = Z, y = y, family = family, 
                penalty = cw_phi * Q, initial.beta = cw_a, offset = new_offset)
        }
        new_a <- fit0$coefficients
        rm(new_offset)
        new.eta <- para.X %*% new_kappa + Z %*% new_a + offset
        if (family$family == "gaussian") {
            new_phi <- (sum((y - new.eta)^2) + sum(ZA * Z))/length(y)
        }
        if (A.struct == "unstructured") {
            if (family$family %in% c("gaussian", "binomialp")) 
                new_A <- chol2inv(chol(Q + (1/new_phi) * ZZt))
            if (family$family == "poisson") {
                err2 <- 10
                cw_A <- new_A
                while (err2 > 0.01) {
                  denom <- Z * matrix(sqrt(exp(new.eta + 0.5 * 
                    rowSums(ZA * Z))), nrow = length(y), ncol = ncol(Z), 
                    byrow = FALSE)
                  denom <- Q + crossprod(denom)
                  new_A <- chol2inv(chol(denom))
                  err2 <- sum((cw_A - new_A)^2)/2
                  cw_A <- new_A
                }
            }
            if (family$family == "binomial") {
                err2 <- 10
                cw_A <- new_A
                while (err2 > 0.01) {
                  denom <- Z * matrix(sqrt(binomial()$linkinv(new.eta + 
                    0.5 * rowSums(ZA * Z))), nrow = length(y), 
                    ncol = ncol(Z), byrow = FALSE)
                  denom <- Q + crossprod(denom)
                  new_A <- chol2inv(chol(denom))
                  err2 <- sum((cw_A - new_A)^2)/2
                  cw_A <- new_A
                }
            }
        }
        if (A.struct == "block") {
            for (k2 in 1:ncov_smooth_X) {
                if (family$family %in% c("gaussian", "binomialp")) {
                  new_A[index_cov == k2, index_cov == k2] <- chol2inv(chol(Q[index_cov == 
                    k2, index_cov == k2] + (1/new_phi) * crossprod(Z[, 
                    index_cov == k2])))
                }
                if (family$family == "poisson") {
                  err2 <- 10
                  while (err2 > 0.01) {
                    denom <- Z[, index_cov == k2] * matrix(sqrt(exp(new.eta + 
                      0.5 * rowSums(ZA * Z))), nrow = length(y), 
                      ncol = sum(index_cov == k2), byrow = FALSE)
                    denom <- Q[index_cov == k2, index_cov == 
                      k2] + crossprod(denom)
                    new_A[index_cov == k2, index_cov == k2] <- chol2inv(chol(denom))
                    err2 <- sum((cw_A - new_A)^2)/2
                    cw_A <- new_A
                  }
                }
                if (family$family == "binomial") {
                  err2 <- 10
                  while (err2 > 0.01) {
                    denom <- Z[, index_cov == k2] * matrix(sqrt(binomial()$linkinv(new.eta + 
                      0.5 * rowSums(ZA * Z))), nrow = length(y), 
                      ncol = sum(index_cov == k2), byrow = FALSE)
                    denom <- Q[index_cov == k2, index_cov == 
                      k2] + crossprod(denom)
                    new_A[index_cov == k2, index_cov == k2] <- chol2inv(chol(denom))
                    err2 <- sum((cw_A - new_A)^2)/2
                    cw_A <- new_A
                  }
                }
            }
        }
        if (is.null(lambda)) {
            for (k2 in 1:ncov_smooth_X) {
                new_lambda[k2] <- d[k2]/(sum(diag(crossprod(S[[k2]], 
                  new_A[index_cov == k2, index_cov == k2]))) + 
                  tcrossprod(new_a[index_cov == k2], S[[k2]]) %*% 
                    new_a[index_cov == k2])
            }
        }
        new_logL <- calc.VAlogL(y = y, Z = Z, para.X = para.X, 
            family = family, lambda = new_lambda, kappa = new_kappa, 
            a = new_a, A = new_A, phi = new_phi, S = S, d = d, 
            index.cov = index_cov, eta = new.eta)
        diff_logL <- new_logL - cw_logL
        diff_lambda <- sum((new_lambda - cw_lambda)^2)
        if (control$trace) 
            cat("Iteration:", counter, "\t Current VA logL:", 
                cw_logL, " | New VA logL:", new_logL, " | Difference:", 
                new_logL - cw_logL, "\n")
        counter <- counter + 1
        cw_logL <- new_logL
        cw_lambda <- new_lambda
        cw_kappa <- new_kappa
        cw_a <- new_a
        cw_A <- new_A
        cw_phi <- new_phi
    }
    toc <- proc.time()
    ic_out <- rep(NA, 2)
    names(ic_out) <- c("AIC", "BIC")
    if (doIC) {
        ic_out <- -2 * new_logL - new_phi + c(2 * fit0$tr.Inf * 
            new_phi, log(length(y)) * fit0$tr.Inf * new_phi)
    }
    out <- list(kappa = new_kappa, a = new_a, A = new_A, lambda = new_lambda, 
        IC = ic_out, phi = new_phi, linear.predictors = new.eta, 
        offset = offset, logL = new_logL, no.knots = actual_int_knots, 
        index.cov = index_cov, basis.info = basis_info, family = family, 
        time.taken = toc - tic)
    if (save.data) {
        out$y
        out$para.X <- para.X
        out$smooth.X <- smooth.X
        out$Z <- Z
    }
    smooth.X.waldstat <- smooth.X.waldpval <- numeric(ncov_smooth_X)
    new_A.inv <- chol2inv(chol(new_A))
    for (k2 in 1:ncov_smooth_X) {
        smooth.X.waldstat[k2] <- crossprod(new_a[index_cov == 
            k2], new_A.inv[index_cov == k2, index_cov == k2]) %*% 
            new_a[index_cov == k2]
        smooth.X.waldpval[k2] <- pchisq(smooth.X.waldstat[k2], 
            df = sum(index_cov == k2), lower.tail = FALSE)
    }
    smooth.X.wald <- round(rbind(smooth.X.waldstat, smooth.X.waldpval), 
        5)
    rownames(smooth.X.wald) <- c("Wald Statistic", "p-value")
    colnames(smooth.X.wald) <- colnames(smooth.X)
    rm(smooth.X.waldstat, smooth.X.waldpval)
    out$smooth.stat <- smooth.X.wald
    if (para.se) {
        message("Calculating information matrix for model parameters...")
        obs_info <- info.valouis(y = y, para.X = para.X, Z = Z, 
            kappa = new_kappa, phi = new_phi, lambda = new_lambda, 
            a = new_a, A = new_A, S = S, d = d, index.cov = index_cov, 
            mc.samps = control$mc.samps, family = family)
        if (any(diag(obs_info) < 0) & family$family == "poisson") {
            message("Redoing...")
            err2 <- 10
            cw_A <- new_A
            while (err2 > 0.01) {
                denom <- Z * matrix(sqrt(exp(new.eta + 0.5 * 
                  rowSums((Z %*% cw_A) * Z))), nrow = length(y), 
                  ncol = ncol(Z), byrow = FALSE)
                denom <- Q + crossprod(denom)
                new_A <- chol2inv(chol(denom))
                err2 <- sum((cw_A - new_A)^2)/2
                cw_A <- new_A
            }
            obs_info <- info.valouis(y = y, para.X = para.X, 
                Z = Z, kappa = new_kappa, phi = new_phi, lambda = new_lambda, 
                a = new_a, A = cw_A, S = S, d = d, index.cov = index_cov, 
                mc.samps = control$mc.samps, family = family)
        }
        while (class(try(solve(obs_info), silent = TRUE)) == 
            "try-error") {
            obs_info <- info.valouis(y = y, para.X = para.X, 
                Z = Z, kappa = new_kappa, phi = new_phi, lambda = new_lambda, 
                a = new_a, A = new_A, S = S, d = d, index.cov = index_cov, 
                mc.samps = control$mc.samps, family = family)
        }
        while (!all(diag(solve(obs_info)) > 0)) {
            obs_info <- info.valouis(y = y, para.X = para.X, 
                Z = Z, kappa = new_kappa, phi = new_phi, lambda = new_lambda, 
                a = new_a, A = new_A, S = S, d = d, index.cov = index_cov, 
                mc.samps = control$mc.samps, family = family)
            while (class(try(solve(obs_info), silent = TRUE)) == 
                "try-error") {
                obs_info <- info.valouis(y = y, para.X = para.X, 
                  Z = Z, kappa = new_kappa, phi = new_phi, lambda = new_lambda, 
                  a = new_a, A = new_A, S = S, d = d, index.cov = index_cov, 
                  mc.samps = control$mc.samps, family = family)
            }
        }
        para_sderr <- sqrt(diag(solve(obs_info)))[1:ncol(para.X)]
        para_X_stat <- round(rbind(new_kappa, para_sderr, new_kappa/para_sderr, 
            2 * pnorm(abs(new_kappa/para_sderr), lower.tail = FALSE)), 
            5)
        rownames(para_X_stat) <- c("Estimate", "Std. Error", 
            "Wald Statistic", "p-value")
        colnames(para_X_stat) <- colnames(para.X)
        out$para.stat <- para_X_stat
        out$obs.info <- obs_info
    }
    if (!para.se) {
        out$para.stat <- NULL
    }
    class(out) <- "vagam"
    out$call <- match.call()
    return(out)
}
<bytecode: 0x564ae35585a8>
<environment: namespace:vagam>
 --- function search by body ---
Function vagam in namespace vagam has this body.
 ----------- END OF FAILURE REPORT -------------- 
Fatal error: the condition has length > 1