CRAN Package Check Results for Package aster

Last updated on 2024-03-28 09:51:36 CET.

Flavor Version Tinstall Tcheck Ttotal Status Flags
r-devel-linux-x86_64-debian-clang 1.1-3 11.14 161.21 172.35 OK
r-devel-linux-x86_64-debian-gcc 1.1-3 9.50 127.04 136.54 OK
r-devel-linux-x86_64-fedora-clang 1.1-3 213.01 OK
r-devel-linux-x86_64-fedora-gcc 1.1-3 225.33 OK
r-devel-windows-x86_64 1.1-3 18.00 2098.00 2116.00 ERROR
r-patched-linux-x86_64 1.1-3 14.01 160.81 174.82 OK
r-release-linux-x86_64 1.1-3 11.72 160.10 171.82 OK
r-release-macos-arm64 1.1-3 109.00 OK
r-release-macos-x86_64 1.1-3 223.00 OK
r-release-windows-x86_64 1.1-3 26.00 206.00 232.00 OK
r-oldrel-macos-arm64 1.1-3 89.00 OK
r-oldrel-windows-x86_64 1.1-3 22.00 208.00 230.00 OK

Check Details

Version: 1.1-3
Check: tests
Result: ERROR Running 'aster.R' [1s] Comparing 'aster.Rout' to 'aster.Rout.save' ... OK Running 'bogus.R' [1s] Comparing 'bogus.Rout' to 'bogus.Rout.save' ... OK Running 'cache.R' [3s] Comparing 'cache.Rout' to 'cache.Rout.save' ... OK Running 'chk.R' [0s] Comparing 'chk.Rout' to 'chk.Rout.save' ... OK Running 'famber.R' [159s] Running 'famnegbin.R' [158s] Running 'famnorloc.R' [169s] Running 'famnzp.R' [0s] Comparing 'famnzp.Rout' to 'famnzp.Rout.save' ... OK Running 'fampoi.R' [0s] Comparing 'fampoi.Rout' to 'fampoi.Rout.save' ... OK Running 'famtnb.R' [0s] Comparing 'famtnb.Rout' to 'famtnb.Rout.save' ... OK Running 'foobar.R' [0s] Comparing 'foobar.Rout' to 'foobar.Rout.save' ... OK Running 'formula.R' [0s] Comparing 'formula.Rout' to 'formula.Rout.save' ... OK Running 'gradmat.R' [158s] Running 'ktnb.R' [164s] Running 'ktp.R' [0s] Comparing 'ktp.Rout' to 'ktp.Rout.save' ... OK Running 'matops.R' [0s] Comparing 'matops.Rout' to 'matops.Rout.save' ... OK Running 'mlogl-c-export.R' [161s] Running 'mlogl-cond.R' [169s] Running 'mlogl-unco.R' [153s] Running 'newfam.R' [0s] Comparing 'newfam.Rout' to 'newfam.Rout.save' ... OK Running 'newpickle.R' [147s] Running 'nzp.R' [0s] Comparing 'nzp.Rout' to 'nzp.Rout.save' ... OK Running 'parm.R' [158s] Running 'penmlogl.R' [1s] Comparing 'penmlogl.Rout' to 'penmlogl.Rout.save' ... OK Running 'pickle.R' [3s] Comparing 'pickle.Rout' to 'pickle.Rout.save' ... OK Running 'predict.R' [1s] Comparing 'predict.Rout' to 'predict.Rout.save' ... OK Running 'quickle.R' [164s] Running 'raster.R' [163s] Running 'reaster-warn.R' [2s] Comparing 'reaster-warn.Rout' to 'reaster-warn.Rout.save' ... OK Running 'reaster.R' [7s] Comparing 'reaster.Rout' to 'reaster.Rout.save' ... OK Running 'reaster1.R' [1s] Comparing 'reaster1.Rout' to 'reaster1.Rout.save' ... OK Running 'trans.R' [0s] Comparing 'trans.Rout' to 'trans.Rout.save' ... OK Running the tests in 'tests/famber.R' failed. Complete output: > > library(aster) > > ifam <- fam.bernoulli() > > p <- seq(0.1, 0.9, 0.1) > theta <- log(p) - log(1 - p) > > zeroth <- double(length(p)) > first <- double(length(p)) > second <- double(length(p)) > > for (i in seq(along = p)) { + zeroth[i] <- famfun(ifam, 0, theta[i]) + first[i] <- famfun(ifam, 1, theta[i]) + second[i] <- famfun(ifam, 2, theta[i]) + } > > all.equal(zeroth, - log(1 - p)) [1] TRUE > all.equal(first, p) [1] TRUE > all.equal(second, p * (1 - p)) [1] TRUE > > > proc.time() user system elapsed 0.25 0.06 0.29 Running the tests in 'tests/famnegbin.R' failed. Complete output: > > library(aster) > > alpha <- 2.222 > ifam <- fam.negative.binomial(alpha) > > ##### usual size theta ##### > > p <- seq(0.9, 0.1, -0.1) > theta <- log(1 - p) > qq <- exp(theta) > pp <- (- expm1(theta)) > all.equal(p, pp) [1] TRUE > all.equal(pp, 1 - qq) [1] TRUE > > zeroth <- double(length(theta)) > first <- double(length(theta)) > second <- double(length(theta)) > > for (i in seq(along = theta)) { + zeroth[i] <- famfun(ifam, 0, theta[i]) + first[i] <- famfun(ifam, 1, theta[i]) + second[i] <- famfun(ifam, 2, theta[i]) + } > > all.equal(zeroth, alpha * (- log(pp))) [1] TRUE > all.equal(first, alpha * qq / pp) [1] TRUE > all.equal(second, alpha * qq / pp^2) [1] TRUE > > ##### very large negative theta ##### > > rm(p) > > theta <- seq(-100, -10, 10) > qq <- exp(theta) > pp <- (- expm1(theta)) > > zeroth <- double(length(theta)) > first <- double(length(theta)) > second <- double(length(theta)) > > for (i in seq(along = theta)) { + zeroth[i] <- famfun(ifam, 0, theta[i]) + first[i] <- famfun(ifam, 1, theta[i]) + second[i] <- famfun(ifam, 2, theta[i]) + } > > all.equal(zeroth, alpha * (- log(pp))) [1] TRUE > all.equal(first, alpha * qq / pp) [1] TRUE > all.equal(second, alpha * qq / pp^2) [1] TRUE > > ##### very small negative theta ##### > > theta <- (- 10^(- c(1:9, seq(10, 100, 10)))) > qq <- exp(theta) > pp <- (- expm1(theta)) > > zeroth <- double(length(theta)) > first <- double(length(theta)) > second <- double(length(theta)) > > for (i in seq(along = theta)) { + zeroth[i] <- famfun(ifam, 0, theta[i]) + first[i] <- famfun(ifam, 1, theta[i]) + second[i] <- famfun(ifam, 2, theta[i]) + } > > all.equal(zeroth, alpha * (- log(pp))) [1] TRUE > all.equal(first, alpha * qq / pp) [1] TRUE > all.equal(second, alpha * qq / pp^2) [1] TRUE > > ##### random ##### > > nind <- 50 > theta <- rep(- 1.75, nind) > pred <- 0 > fam <- 1 > root <- seq(1, by = 0.5, length = nind) > theta <- cbind(theta) > root <- cbind(root) > > set.seed(42) > rout <- raster(theta, pred, fam, root, famlist = list(ifam)) > > set.seed(42) > rout.too <- rnbinom(nind, size = alpha * as.numeric(root), + prob = as.numeric(1 - exp(theta))) > > all.equal(as.numeric(rout), rout.too) [1] TRUE > > > > > proc.time() user system elapsed 0.23 0.06 0.28 Running the tests in 'tests/famnorloc.R' failed. Complete output: > > library(aster) > > sigma <- 2.222 > ifam <- fam.normal.location(sigma) > > ##### usual size theta ##### > > theta <- seq(-3.0, 3.0, 0.1) > > zeroth <- double(length(theta)) > first <- double(length(theta)) > second <- double(length(theta)) > > for (i in seq(along = theta)) { + zeroth[i] <- famfun(ifam, 0, theta[i]) + first[i] <- famfun(ifam, 1, theta[i]) + second[i] <- famfun(ifam, 2, theta[i]) + } > > all.equal(zeroth, sigma^2 * theta^2 / 2) [1] TRUE > all.equal(first, sigma^2 * theta) [1] TRUE > all.equal(second, sigma^2 * theta^0) [1] TRUE > > ##### random ##### > > theta <- seq(-3.5, 3.5, 0.1) > nind <- length(theta) > pred <- 0 > fam <- 1 > root <- seq(1, 5, length = nind) > theta <- cbind(theta) > root <- cbind(root) > > set.seed(42) > rout <- raster(theta, pred, fam, root, famlist = list(ifam)) > > set.seed(42) > moo <- sigma^2 * theta * as.numeric(root) > cow <- sigma * sqrt(as.numeric(root)) > rout.too <- rnorm(nind, mean = moo, sd = cow) > > all.equal(as.numeric(rout), rout.too) [1] TRUE > > > proc.time() user system elapsed 0.23 0.06 0.28 Running the tests in 'tests/gradmat.R' failed. Complete output: > > library(aster) > > # needed because of the change in R function "sample" in R-devel > suppressWarnings(RNGversion("3.5.2")) > > set.seed(42) > > nind <- 25 > > vars <- c("l2", "l3", "f2", "f3", "h2", "h3") > pred <- c(0, 1, 1, 2, 3, 4) > fam <- c(1, 1, 1, 1, 3, 3) > length(pred) == length(fam) [1] TRUE > nnode <- length(pred) > > theta <- matrix(0, nind, nnode) > root <- matrix(1, nind, nnode) > x <- raster(theta, pred, fam, root) > dimnames(x) <- list(NULL, vars) > > data <- as.data.frame(x) > site <- factor(sample(LETTERS[1:4], nind, replace = TRUE)) > foo <- rnorm(nind) > data <- data.frame(x, site = site, foo = foo, root = 1) > > redata <- reshape(data, varying = list(vars), + direction = "long", timevar = "varb", times = as.factor(vars), + v.names = "resp") > > out <- aster(resp ~ foo + site + varb, pred, fam, varb, id, root, + data = redata) > sout1 <- summary(out, show.graph = TRUE) > > out2 <- aster(x, root, pred, fam, modmat = out$modmat) > sout2 <- summary(out2) > > out3 <- aster(x, root, pred, fam, modmat = out$modmat, type = "cond") > sout3 <- summary(out3) > > foo <- new.env(parent = emptyenv()) > bar <- suppressWarnings(try(load("gradmat.rda", foo), silent = TRUE)) > if (inherits(bar, "try-error")) { + save(sout1, sout2, sout3, file = "gradmat.rda") + } else { + print(all.equal(sout1, foo$sout1)) + print(all.equal(sout2, foo$sout2)) + print(all.equal(sout3, foo$sout3)) + } [1] TRUE [1] TRUE [1] TRUE > > foo <- match(sort(unique(site)), site) > modmat.pred <- out$modmat[foo, , ] > > ##### case 1: eta = theta, zeta = phi > > beta.hat <- out3$coef > > modmat.pred.mat <- matrix(modmat.pred, ncol = length(beta.hat)) > > theta.hat <- matrix(modmat.pred.mat %*% beta.hat, nrow = dim(modmat.pred)[1]) > > nind <- dim(modmat.pred)[1] > nnode <- dim(modmat.pred)[2] > ncoef <- dim(modmat.pred)[3] > > aster:::setfam(fam.default()) > > phi.hat <- .C(aster:::C_aster_theta2phi, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + theta = as.double(theta.hat), + phi = matrix(as.double(0), nind, nnode))$phi > > my.phi.hat <- theta.hat > my.phi.hat[ , 4] <- my.phi.hat[ , 4] - log(exp(exp(theta.hat[ , 6])) - 1) > my.phi.hat[ , 3] <- my.phi.hat[ , 3] - log(exp(exp(theta.hat[ , 5])) - 1) > my.phi.hat[ , 2] <- my.phi.hat[ , 2] - log(1 + exp(theta.hat[ , 4])) > my.phi.hat[ , 1] <- my.phi.hat[ , 1] - log(1 + exp(theta.hat[ , 3])) > my.phi.hat[ , 1] <- my.phi.hat[ , 1] - log(1 + exp(theta.hat[ , 2])) > all.equal(my.phi.hat, phi.hat) [1] TRUE > > gradmat <- 0 * modmat.pred > storage.mode(gradmat) <- "double" > > gradmat <- .C(aster:::C_aster_D_beta2theta2phi, + nind = as.integer(nind), + nnode = as.integer(nnode), + ncoef = as.integer(ncoef), + pred = as.integer(pred), + fam = as.integer(fam), + theta = as.double(theta.hat), + modmat = as.double(modmat.pred), + gradmat = gradmat)$gradmat > > my.gradmat <- 0 * gradmat > epsilon <- 1e-9 > for (k in 1:ncoef) { + beta.epsilon <- beta.hat + beta.epsilon[k] <- beta.hat[k] + epsilon + theta.epsilon <- matrix(modmat.pred.mat %*% beta.epsilon, nrow = nind) + phi.epsilon <- .C(aster:::C_aster_theta2phi, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + theta = as.double(theta.epsilon), + phi = matrix(as.double(0), nind, nnode))$phi + my.gradmat[ , , k] <- (phi.epsilon - phi.hat) / epsilon + } > > all.equal(gradmat, my.gradmat, tolerance = sqrt(epsilon)) [1] TRUE > > ##### case 2: eta = phi, zeta = theta > > beta.hat <- out2$coef > > phi.hat <- matrix(modmat.pred.mat %*% beta.hat, nrow = nind) > > theta.hat <- .C(aster:::C_aster_phi2theta, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + phi = as.double(phi.hat), + theta = matrix(as.double(0), nind, nnode))$theta > > gradmat <- .C(aster:::C_aster_D_beta2phi2theta, + nind = as.integer(nind), + nnode = as.integer(nnode), + ncoef = as.integer(ncoef), + pred = as.integer(pred), + fam = as.integer(fam), + theta = as.double(theta.hat), + modmat = as.double(modmat.pred), + gradmat = gradmat)$gradmat > > my.gradmat <- 0 * gradmat > epsilon <- 1e-9 > for (k in 1:ncoef) { + beta.epsilon <- beta.hat + beta.epsilon[k] <- beta.hat[k] + epsilon + phi.epsilon <- matrix(modmat.pred.mat %*% beta.epsilon, nrow = nind) + theta.epsilon <- .C(aster:::C_aster_phi2theta, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + phi = as.double(phi.epsilon), + theta = matrix(as.double(0), nind, nnode))$theta + my.gradmat[ , , k] <- (theta.epsilon - theta.hat) / epsilon + } > > all.equal(gradmat, my.gradmat, tolerance = sqrt(epsilon)) [1] TRUE > > ##### case 3: eta = phi, zeta = tau > > root.pred <- matrix(1, nind, nnode) > > beta.hat <- out2$coef > > beta2tau <- function(beta) { + + phi <- matrix(modmat.pred.mat %*% beta, nrow = nind) + + theta <- .C(aster:::C_aster_phi2theta, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + phi = as.double(phi), + theta = matrix(as.double(0), nind, nnode))$theta + + ctau <- .C(aster:::C_aster_theta2ctau, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + theta = as.double(theta), + ctau = double(nind * nnode))$ctau + + tau <- .C(aster:::C_aster_ctau2tau, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + root = as.double(root.pred), + ctau = as.double(ctau), + tau = double(nind * nnode))$tau + + return(tau) + } > > tau.hat <- beta2tau(beta.hat) > > gradmat <- .C(aster:::C_aster_D_beta2phi2tau, + nind = as.integer(nind), + nnode = as.integer(nnode), + ncoef = as.integer(ncoef), + pred = as.integer(pred), + fam = as.integer(fam), + beta = as.double(beta.hat), + root = as.double(root.pred), + origin = rep(as.double(0), nind * nnode), + modmat = as.double(modmat.pred), + gradmat = gradmat)$gradmat > > my.gradmat <- 0 * gradmat > epsilon <- 1e-9 > for (k in 1:ncoef) { + beta.epsilon <- beta.hat + beta.epsilon[k] <- beta.hat[k] + epsilon + tau.epsilon <- beta2tau(beta.epsilon) + my.gradmat[ , , k] <- (tau.epsilon - tau.hat) / epsilon + } > > all.equal(gradmat, my.gradmat, tolerance = sqrt(epsilon)) [1] TRUE > > ##### case 4: eta = theta, zeta = tau > > beta.hat <- out3$coef > > beta2tau <- function(beta) { + + theta <- matrix(modmat.pred.mat %*% beta, nrow = nind) + + ctau <- .C(aster:::C_aster_theta2ctau, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + theta = as.double(theta), + ctau = double(nind * nnode))$ctau + + tau <- .C(aster:::C_aster_ctau2tau, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + root = as.double(root.pred), + ctau = as.double(ctau), + tau = double(nind * nnode))$tau + + return(tau) + } > > tau.hat <- beta2tau(beta.hat) > > gradmat <- .C(aster:::C_aster_D_beta2theta2tau, + nind = as.integer(nind), + nnode = as.integer(nnode), + ncoef = as.integer(ncoef), + pred = as.integer(pred), + fam = as.integer(fam), + beta = as.double(beta.hat), + root = as.double(root.pred), + modmat = as.double(modmat.pred), + gradmat = gradmat)$gradmat > > my.gradmat <- 0 * gradmat > for (k in 1:ncoef) { + beta.epsilon <- beta.hat + beta.epsilon[k] <- beta.hat[k] + epsilon + tau.epsilon <- beta2tau(beta.epsilon) + my.gradmat[ , , k] <- (tau.epsilon - tau.hat) / epsilon + } > > all.equal(gradmat, my.gradmat, tolerance = sqrt(epsilon)) [1] TRUE > > > proc.time() user system elapsed 0.34 0.03 0.36 Running the tests in 'tests/ktnb.R' failed. Complete output: > > library(aster) > > do.chisq.test <- function(x, alpha, k, mu, max.bin) { + stopifnot(all(x > k)) + stopifnot(k + 1 < max.bin) + xx <- seq(k + 1, max.bin) + yy <- dnbinom(xx, size = alpha, mu = mu) + yy[length(yy)] <- pnbinom(max.bin - 1, size = alpha, mu = mu, + lower.tail = FALSE) + pp <- yy / sum(yy) + ecc <- length(x) * pp + if (any(ecc < 5.0)) + warning("violates rule of thumb about > 5 expected in each cell") + cc <- tabulate(x, max.bin) + cc <- cc[xx] + cc[length(cc)] <- nsim - sum(cc[- length(cc)]) + chisqstat <- sum((cc - ecc)^2 / ecc) + pval <- pchisq(chisqstat, length(ecc) - 1, lower.tail = FALSE) + if (exists("save.min.pval")) { + save.min.pval <<- min(pval, save.min.pval) + save.ntests <<- save.ntests + 1 + } else { + save.min.pval <<- pval + save.ntests <<- 1 + } + list(chisqstat = chisqstat, df = length(ecc) - 1, pval = pval, + observed = cc, expected = ecc, x = xx) + } > > set.seed(42) > nsim <- 1e4 > > alpha <- 2.222 > mu <- 10 > k <- 2 > x <- rktnb(nsim, alpha, k, mu) > chisqout1 <- do.chisq.test(x, alpha, k, mu, 40) > > alpha <- 2.222 > mu <- 3.5 > k <- 2 > x <- rktnb(nsim, alpha, k, mu) > chisqout2 <- do.chisq.test(x, alpha, k, mu, 20) > > alpha <- 2.222 > mu <- 2.5 > k <- 2 > x <- rktnb(nsim, alpha, k, mu) > chisqout3 <- do.chisq.test(x, alpha, k, mu, 16) > > alpha <- 2.222 > mu <- 1.5 > k <- 2 > x <- rktnb(nsim, alpha, k, mu) > chisqout4 <- do.chisq.test(x, alpha, k, mu, 12) > > alpha <- 2.222 > mu <- 0.5 > k <- 2 > x <- rktnb(nsim, alpha, k, mu) > chisqout5 <- do.chisq.test(x, alpha, k, mu, 8) > > alpha <- 2.222 > mu <- 0.1 > k <- 2 > x <- rktnb(nsim, alpha, k, mu) > chisqout6 <- do.chisq.test(x, alpha, k, mu, 5) > > nsim <- 2e5 > alpha <- 2.222 > mu <- 0.01 > k <- 2 > x <- rktnb(nsim, alpha, k, mu) > chisqout7 <- do.chisq.test(x, alpha, k, mu, 5) > > alpha <- 2.222 > mu <- 1.5 > xpred <- 0:10 > save.seed <- .Random.seed > x <- rktp(xpred, k, mu, xpred) > .Random.seed <- save.seed > my.x <- rep(0, length(xpred)) > for (i in seq(along = xpred)) + if (xpred[i] > 0) + for (j in 1:xpred[i]) + my.x[i] <- my.x[i] + rktp(1, k, mu) > all.equal(x, my.x) [1] TRUE > > nsim <- 1e4 > alpha <- 5.55 > k <- 5 > mu <- pi > x <- rktnb(nsim, alpha, k, mu) > chisqout8 <- do.chisq.test(x, alpha, k, mu, 16) > > alpha <- 5.55 > k <- 10 > mu <- exp(2) > x <- rktnb(nsim, alpha, k, mu) > chisqout9 <- do.chisq.test(x, alpha, k, mu, 29) > > cat("number of tests:", save.ntests, "\n") number of tests: 9 > save.ntests * save.min.pval > 0.05 [1] TRUE > > ##### > > set.seed(42) > nind <- 25 > nnode <- 1 > ncoef <- 1 > alpha <- 3.333 > k <- 2 > > pred <- 0 > fam <- 1 > ifam <- fam.truncated.negative.binomial(size = alpha, trunc = k) > aster:::setfam(list(ifam)) > theta.origin <- aster:::getfam()[[1]]$origin > > theta <- (- 4 / 3) > p <- 1 - exp(theta) > x <- rnbinom(1000, size = alpha, prob = p) > x <- x[x > k] > x <- x[1:nind] > modmat <- matrix(1, nrow = nind, ncol = 1) > > out <- mlogl(theta - theta.origin, pred, fam, x, modmat, modmat, + deriv = 2, type = "conditional", famlist = list(ifam)) > > xxx <- seq(0, 100) > ppp <- dnbinom(xxx, size = alpha, prob = p) > ppp[xxx <= k] <- 0 > ppp <- ppp / sum(ppp) > tau <- sum(xxx * ppp) > > my.grad.logl <- sum(x - tau) > all.equal(- out$gradient, my.grad.logl) [1] TRUE > > my.fish.info <- length(x) * sum((xxx - tau)^2 * ppp) > all.equal(as.numeric(out$hessian), my.fish.info) [1] TRUE > > foo <- new.env(parent = emptyenv()) > bar <- suppressWarnings(try(load("ktnb.rda", foo), silent = TRUE)) > if (inherits(bar, "try-error")) { + save(list = c(paste("chisqout", 1:9, sep = ""), "out"), file = "ktnb.rda") + } else { + print(all.equal(chisqout1, foo$chisqout1)) + print(all.equal(chisqout2, foo$chisqout2)) + print(all.equal(chisqout3, foo$chisqout3)) + print(all.equal(chisqout4, foo$chisqout4)) + print(all.equal(chisqout5, foo$chisqout5)) + print(all.equal(chisqout6, foo$chisqout6)) + print(all.equal(chisqout7, foo$chisqout7)) + print(all.equal(chisqout8, foo$chisqout8)) + print(all.equal(chisqout9, foo$chisqout9)) + print(all.equal(out, foo$out)) + } [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE > > > proc.time() user system elapsed 0.37 0.09 0.42 Running the tests in 'tests/mlogl-c-export.R' failed. Complete output: > > # copied from ../man/aster.Rd examples section > > library(aster) > data(echinacea) > vars <- c("ld02", "ld03", "ld04", "fl02", "fl03", "fl04", + "hdct02", "hdct03", "hdct04") > redata <- reshape(echinacea, varying = list(vars), direction = "long", + timevar = "varb", times = as.factor(vars), v.names = "resp") > redata <- data.frame(redata, root = 1) > pred <- c(0, 1, 2, 1, 2, 3, 4, 5, 6) > fam <- c(1, 1, 1, 1, 1, 1, 3, 3, 3) > hdct <- grepl("hdct", as.character(redata$varb)) > redata <- data.frame(redata, hdct = as.integer(hdct)) > level <- gsub("[0-9]", "", as.character(redata$varb)) > redata <- data.frame(redata, level = as.factor(level)) > aout <- aster(resp ~ varb + level : (nsloc + ewloc) + hdct : pop, + pred, fam, varb, id, root, data = redata) > > # redo with new function > > theta <- predict(aout, parm.type = "canonical", model.type = "conditional") > phi <- predict(aout, parm.type = "canonical", model.type = "unconditional") > > nind <- nrow(echinacea) > nnode <- length(pred) > > aster:::setfam(fam.default()) > > val.theta <- .C(aster:::C_aster_export_exerciser, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + parm = as.double(theta), + root = as.double(redata$root), + response = as.double(redata$resp), + is.unco = FALSE, + value = double(1))$value > > val.phi <- .C(aster:::C_aster_export_exerciser, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + parm = as.double(phi), + root = as.double(redata$root), + response = as.double(redata$resp), + is.unco = TRUE, + value = double(1))$value > > aster:::clearfam() > > all.equal(val.theta, aout$deviance / 2) [1] TRUE > all.equal(val.phi, aout$deviance / 2) [1] TRUE > > > proc.time() user system elapsed 0.68 0.04 0.71 Running the tests in 'tests/mlogl-cond.R' failed. Complete output: > > library(aster) > library(numDeriv) > > # needed because of the change in R function "sample" in R-devel > suppressWarnings(RNGversion("3.5.2")) > > set.seed(42) > > nind <- 25 > nnode <- 5 > ncoef <- nnode + 1 > > famlist <- fam.default() > fam <- c(1, 1, 2, 3, 3) > pred <- c(0, 1, 1, 2, 3) > > modmat <- array(0, c(nind, nnode, ncoef)) > modmat[ , , 1] <- 1 > for (i in 2:nnode) + modmat[ , i, i] <- 1 > modmat[ , , ncoef] <- rnorm(nind * nnode) > > beta <- rnorm(ncoef) / 10 > > theta <- matrix(modmat, ncol = ncoef) %*% beta > theta <- matrix(theta, ncol = nnode) > > root <- sample(1:3, nind * nnode, replace = TRUE) > root <- matrix(root, nind, nnode) > > x <- raster(theta, pred, fam, root) > > out <- mlogl(beta, pred, fam, x, root, modmat, deriv = 2, + type = "conditional") > > my.value <- 0 > for (j in 1:nnode) { + ifam <- fam[j] + k <- pred[j] + if (k > 0) + xpred <- x[ , k] + else + xpred <- root[ , j] + for (i in 1:nind) + my.value <- my.value - + sum(x[i, j] * theta[i, j] - + xpred[i] * famfun(famlist[[ifam]], 0, theta[i, j])) + } > all.equal(out$value, my.value) [1] TRUE > > foo <- function(beta) { + stopifnot(is.numeric(beta)) + stopifnot(is.finite(beta)) + mlogl(beta, pred, fam, x, root, modmat, type = "conditional")$value + } > > g <- grad(foo, beta) > all.equal(g, out$gradient) [1] TRUE > > h <- hessian(foo, beta) > all.equal(h, out$hessian) [1] TRUE > > ########## > > objfun <- function(beta) { + out <- mlogl(beta, pred, fam, x, root, modmat, deriv = 1, + type = "conditional") + result <- out$value + attr(result, "gradient") <- out$gradient + return(result) + } > out1 <- nlm(objfun, beta, fscale = nind) > > ########## > > objfun <- function(beta) { + out <- mlogl(beta, pred, fam, x, root, modmat, deriv = 2, + type = "conditional") + result <- out$value + attr(result, "gradient") <- out$gradient + attr(result, "hessian") <- out$hessian + return(result) + } > out <- nlm(objfun, beta, fscale = nind) > all.equal(out1$minimum, out$minimum) [1] TRUE > all.equal(out1$estimate, out$estimate, tolerance = 1e-4) [1] TRUE > > ########## expected Fisher information ########## > > aster:::setfam(fam.default()) > > fout <- .C(aster:::C_aster_fish_cond, + nind = as.integer(nind), + nnode = as.integer(nnode), + ncoef = as.integer(ncoef), + pred = as.integer(pred), + fam = as.integer(fam), + beta = as.double(out$estimate), + root = as.double(root), + x = as.double(x), + modmat = as.double(modmat), + fish = matrix(as.double(0), ncoef, ncoef)) > > mout <- mlogl(out$estimate, pred, fam, x, root, modmat, + deriv = 2, type = "conditional") > > aster:::setfam(fam.default()) > > theta <- .C(aster:::C_aster_mat_vec_mult, + nrow = as.integer(nind * nnode), + ncol = as.integer(ncoef), + a = as.double(modmat), + b = as.double(out$estimate), + c = matrix(as.double(0), nind, nnode))$c > ctau <- .C(aster:::C_aster_theta2ctau, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + theta = as.double(theta), + ctau = matrix(as.double(0), nind, nnode))$ctau > tau <- .C(aster:::C_aster_ctau2tau, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + root = as.double(root), + ctau = as.double(ctau), + tau = matrix(as.double(0), nind, nnode))$tau > psiDoublePrime <- .C(aster:::C_aster_theta2whatsis, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + deriv = as.integer(2), + theta = as.double(theta), + result = matrix(as.double(0), nind, nnode))$result > > my.hessian <- theta * NaN > my.fish <- theta * NaN > > for (i in 1:nind) { + for (j in 1:nnode) { + k <- pred[j] + if (k > 0) { + my.hessian[i, j] <- x[i, k] * psiDoublePrime[i, j] + my.fish[i, j] <- tau[i, k] * psiDoublePrime[i, j] + } else { + my.hessian[i, j] <- root[i, j] * psiDoublePrime[i, j] + my.fish[i, j] <- root[i, j] * psiDoublePrime[i, j] + } + } + } > > modmat <- matrix(as.double(modmat), ncol = ncoef) > my.hessian <- as.numeric(my.hessian) > my.fish <- as.numeric(my.fish) > my.hessian <- t(modmat) %*% diag(my.hessian) %*% modmat > my.fish <- t(modmat) %*% diag(my.fish) %*% modmat > > all.equal(my.hessian, mout$hessian) [1] TRUE > all.equal(my.fish, fout$fish) [1] TRUE > > > proc.time() user system elapsed 0.48 0.07 0.50 Running the tests in 'tests/mlogl-unco.R' failed. Complete output: > > library(aster) > library(numDeriv) > > set.seed(42) > > # needed because of the change in R function "sample" in R-devel > suppressWarnings(RNGversion("3.5.2")) > > nind <- 25 > nnode <- 5 > ncoef <- nnode + 1 > > famlist <- fam.default() > fam <- c(1, 1, 2, 3, 3) > pred <- c(0, 1, 1, 2, 3) > > modmat <- array(0, c(nind, nnode, ncoef)) > modmat[ , , 1] <- 1 > for (i in 2:nnode) + modmat[ , i, i] <- 1 > modmat[ , , ncoef] <- rnorm(nind * nnode) > > beta <- rnorm(ncoef) / 10 > > phi <- matrix(modmat, ncol = ncoef) %*% beta > phi <- matrix(phi, ncol = nnode) > > aster:::setfam(fam.default()) > > theta <- .C(aster:::C_aster_phi2theta, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + phi = as.double(phi), + theta = matrix(as.double(0), nind, nnode))$theta > > root <- sample(1:3, nind * nnode, replace = TRUE) > root <- matrix(root, nind, nnode) > > x <- raster(theta, pred, fam, root) > > zip <- rep(0, nind * nnode) > > out <- mlogl(beta, pred, fam, x, root, modmat, deriv = 2, + type = "unco", origin = zip) > > aster:::setfam(fam.default()) > > a <- .C(aster:::C_aster_theta2phi, + nind = as.integer(nind), + nnode = as.integer(nnode), + pred = as.integer(pred), + fam = as.integer(fam), + theta = as.double(zip), + phi = matrix(as.double(0), nind, nnode))$phi > > M <- matrix(modmat, ncol = ncoef) > > alpha <- as.numeric(lm(as.numeric(a) ~ 0 + M)$coefficients) > > out.too <- mlogl(beta - alpha, pred, fam, x, root, modmat, deriv = 2, + type = "unco") > all.equal(out, out.too) [1] TRUE > > beta.old <- beta > beta <- beta - alpha > > my.value <- 0 > for (j in 1:nnode) { + ifam <- fam[j] + k <- pred[j] + if (k > 0) + xpred <- x[ , k] + else + xpred <- root[ , j] + for (i in 1:nind) + my.value <- my.value - + sum(x[i, j] * theta[i, j] - + xpred[i] * famfun(famlist[[ifam]], 0, theta[i, j])) + } > all.equal(out$value, my.value) [1] TRUE > > foo <- function(beta) { + stopifnot(is.numeric(beta)) + stopifnot(is.finite(beta)) + mlogl(beta, pred, fam, x, root, modmat, type = "unco")$value + } > > g <- grad(foo, beta) > all.equal(g, out$gradient) [1] TRUE > > h <- hessian(foo, beta) > all.equal(h, out$hessian) [1] TRUE > > ########## > > objfun <- function(beta) { + out <- mlogl(beta, pred, fam, x, root, modmat, deriv = 1, + type = "unco") + result <- out$value + attr(result, "gradient") <- out$gradient + return(result) + } > nout1 <- nlm(objfun, beta, fscale = nind) Warning message: In nlm(objfun, beta, fscale = nind) : NA/Inf replaced by maximum positive value > nout <- nlm(objfun, nout1$estimate, fscale = nind) > all.equal(nout1$minimum, nout$minimum) [1] TRUE > all.equal(nout1$estimate, nout$estimate) [1] TRUE > > beta.mle.new <- nout$estimate > beta.mle.old <- beta.mle.new + alpha > mout.new <- mlogl(beta.mle.new, pred, fam, x, root, modmat, deriv = 1, + type = "unco") > mout.old <- mlogl(beta.mle.old, pred, fam, x, root, modmat, deriv = 1, + type = "unco", origin = zip) > all.equal(mout.new, mout.old, tol = 1e-7) [1] TRUE > > ########## > > objfun <- function(beta) { + out <- mlogl(beta, pred, fam, x, root, modmat, deriv = 2, + type = "unco") + result <- out$value + attr(result, "gradient") <- out$gradient + attr(result, "hessian") <- out$hessian + return(result) + } > nout1 <- nlm(objfun, beta, fscale = nind) > nout <- nlm(objfun, nout1$estimate, fscale = nind, iterlim = 1000) > all.equal(nout1$minimum, nout$minimum, tol = 1e-4) [1] TRUE > all.equal(nout1$estimate, nout$estimate, tol = 2e-2) [1] TRUE > > objfun.old <- function(beta) { + out <- mlogl(beta, pred, fam, x, root, modmat, deriv = 2, + type = "unco", origin = zip) + result <- out$value + attr(result, "gradient") <- out$gradient + attr(result, "hessian") <- out$hessian + return(result) + } > nout.old <- nlm(objfun.old, beta.mle.old, fscale = nind, iterlim = 1000) > all.equal(nout$minimum, nout.old$minimum) [1] TRUE > all.equal(nout$estimate, nout.old$estimate - alpha, tol = 1e-4) [1] TRUE > > > proc.time() user system elapsed 0.40 0.10 0.46 Running the tests in 'tests/newpickle.R' failed. Complete output: > > library(aster) > > data(radish) > > pred <- c(0,1,2) > fam <- c(1,3,2) > > ### need object of type aster to supply to penmlogl and pickle > > aout <- aster(resp ~ varb + fit : (Site * Region + Block + Pop), + pred, fam, varb, id, root, data = radish) > > ### model matrices for fixed and random effects > > modmat.fix <- model.matrix(resp ~ varb + fit : (Site * Region), + data = radish) > modmat.blk <- model.matrix(resp ~ 0 + fit:Block, data = radish) > modmat.pop <- model.matrix(resp ~ 0 + fit:Pop, data = radish) > > rownames(modmat.fix) <- NULL > rownames(modmat.blk) <- NULL > rownames(modmat.pop) <- NULL > > idrop <- match(aout$dropped, colnames(modmat.fix)) > idrop <- idrop[! is.na(idrop)] > modmat.fix <- modmat.fix[ , - idrop] > > nfix <- ncol(modmat.fix) > nblk <- ncol(modmat.blk) > npop <- ncol(modmat.pop) > > ### try newpickle > > sigma.start <- c(1, 1) > alpha.start <- aout$coefficients[match(colnames(modmat.fix), + names(aout$coefficients))] > cee.start <- rep(0, nblk + npop) > alphaceesigma.start <- c(alpha.start, cee.start, sigma.start) > > foo <- newpickle(alphaceesigma.start, fixed = modmat.fix, + random = list(modmat.blk, modmat.pop), obj = aout) > names(foo) [1] "value" > > > proc.time() user system elapsed 0.37 0.14 0.50 Running the tests in 'tests/parm.R' failed. Complete output: > > library(aster) > > data(echinacea) > > # cut and paste from help(predict.aster) > vars <- c("ld02", "ld03", "ld04", "fl02", "fl03", "fl04", + "hdct02", "hdct03", "hdct04") > redata <- reshape(echinacea, varying = list(vars), direction = "long", + timevar = "varb", times = as.factor(vars), v.names = "resp") > redata <- data.frame(redata, root = 1) > pred <- c(0, 1, 2, 1, 2, 3, 4, 5, 6) > fam <- c(1, 1, 1, 1, 1, 1, 3, 3, 3) > hdct <- grepl("hdct", as.character(redata$varb)) > redata <- data.frame(redata, hdct = as.integer(hdct)) > level <- gsub("[0-9]", "", as.character(redata$varb)) > redata <- data.frame(redata, level = as.factor(level)) > aout <- aster(resp ~ varb + level : (nsloc + ewloc) + hdct : pop, + pred, fam, varb, id, root, data = redata) > # end of cut and paste from help(predict.aster) > > length(aout$dropped) [1] 1 > # so would have had a problem with restart in version 1.0-3 and before > aout <- aster(resp ~ varb + level : (nsloc + ewloc) + hdct : pop, + pred, fam, varb, id, root, data = redata, parm = aout$coefficients, + optout = TRUE) > aout$optout$converged [1] TRUE > aout$optout$iterations == 1 [1] TRUE > > > proc.time() user system elapsed 0.95 0.12 1.06 Running the tests in 'tests/quickle.R' failed. Complete output: > > library(aster) > > data(radish) > > pred <- c(0,1,2) > fam <- c(1,3,2) > > ### try reaster > > rout <- reaster(resp ~ varb + fit : (Site * Region), + list(block = ~ 0 + fit : Block, pop = ~ 0 + fit : Pop), + pred, fam, varb, id, root, data = radish) > srout <- summary(rout) > names(rout) [1] "obj" "fixed" "random" "dropped" "sigma" [6] "nu" "c" "b" "alpha" "zwz" [11] "response" "origin" "iterations" "counts" "deviance" [16] "formula" "call" > > foo <- new.env(parent = emptyenv()) > bar <- suppressWarnings(try(load("quickle.rda", foo), silent = TRUE)) > if (inherits(bar, "try-error")) { + save(srout, file = "quickle.rda") + } else { + srout$object$iterations <- NULL + foo$srout$object$iterations <- NULL + print(all.equal(srout, foo$srout, tol = 1e-4)) + } [1] TRUE > > alpha.mle <- rout$alpha > bee.mle <- rout$b > nu.mle <- rout$sigma^2 > zwz.mle <- rout$zwz > obj <- rout$obj > fixed <- rout$fixed > random <- rout$random > alphanu.mle <- c(alpha.mle, nu.mle) > > qout <- quickle(alphanu.mle, bee.mle, fixed, random, obj, + zwz = zwz.mle, deriv = 2) > > objfun <- function(alphanu) quickle(alphanu, bee.mle, fixed, random, + obj, zwz = zwz.mle)$value > gradfun <- function(alphanu) quickle(alphanu, bee.mle, fixed, random, + obj, zwz = zwz.mle, deriv = 1)$gradient > oout <- optim(alphanu.mle, objfun, gradfun, method = "BFGS", hessian = TRUE) > all.equal(qout$hessian, oout$hessian, check.attributes = FALSE, + tolerance = 0.002) [1] TRUE > > > proc.time() user system elapsed 3.35 0.40 3.79 Running the tests in 'tests/raster.R' failed. Complete output: > > library(aster) > > # needed because of the change in R function "sample" in R-devel > suppressWarnings(RNGversion("3.5.2")) > > set.seed(42) > nind <- 25 > nnode <- 5 > > theta <- rnorm(nind * nnode) / 10 > theta <- matrix(theta, nind, nnode) > > root <- sample(1:3, nind * nnode, replace = TRUE) > root <- matrix(root, nind, nnode) > > fam <- c(1, 1, 2, 3, 3) > pred <- c(0, 1, 1, 2, 3) > famnam <- sapply(fam.default(), as.character) > > save.seed <- .Random.seed > x <- raster(theta, pred, fam, root) > > .Random.seed <- save.seed > my.x <- NaN * x > > for (i in 1:nnode) { + ipred <- pred[i] + if (ipred == 0) { + xpred <- root[ , i] + } else { + xpred <- my.x[ , ipred] + } + if (famnam[fam[i]] == "bernoulli") { + p <- 1 / (1 + exp(- theta[ , i])) + xi <- rep(0, nind) + ii <- xpred > 0 + xi[ii] <- rbinom(sum(ii), xpred[ii], p[ii]) + my.x[ , i] <- xi + } + if (famnam[fam[i]] == "poisson") { + mu <- exp(theta[ , i]) + xi <- rep(0, nind) + ii <- xpred > 0 + xi[ii] <- rpois(sum(ii), xpred[ii] * mu[ii]) + my.x[ , i] <- xi + } + if (famnam[fam[i]] == "truncated.poisson(truncation = 0)") { + mu <- exp(theta[ , i]) + my.x[ , i] <- rnzp(nind, mu, xpred) + } + } > > all.equal(x, my.x) [1] TRUE > > > proc.time() user system elapsed 0.21 0.09 0.29 Flavor: r-devel-windows-x86_64