* using log directory ‘/data/gannet/ripley/R/packages/tests-LENGTH1/multinbmod.Rcheck’ * using R Under development (unstable) (2022-04-03 r82074) * using platform: x86_64-pc-linux-gnu (64-bit) * using session charset: UTF-8 * using option ‘--no-stop-on-test-error’ * checking for file ‘multinbmod/DESCRIPTION’ ... OK * this is package ‘multinbmod’ version ‘1.0’ * checking package namespace information ... OK * checking package dependencies ... OK * checking if this is a source package ... OK * checking if there is a namespace ... OK * checking for executable files ... OK * checking for hidden files and directories ... OK * checking for portable file names ... OK * checking for sufficient/correct file permissions ... OK * checking whether package ‘multinbmod’ can be installed ... OK * checking package directory ... OK * checking DESCRIPTION meta-information ... OK * checking top-level files ... OK * checking for left-over files ... OK * checking index information ... OK * checking package subdirectories ... OK * checking R files for non-ASCII characters ... OK * checking R files for syntax errors ... OK * checking whether the package can be loaded ... OK * checking whether the package can be loaded with stated dependencies ... OK * checking whether the package can be unloaded cleanly ... OK * checking whether the namespace can be loaded with stated dependencies ... OK * checking whether the namespace can be unloaded cleanly ... OK * checking loading without being on the library search path ... OK * checking use of S3 registration ... OK * checking dependencies in R code ... OK * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK * checking R code for possible problems ... NOTE multinb.fit: no visible global function definition for ‘nlminb’ multinbmod: no visible global function definition for ‘is.empty.model’ multinbmod: no visible global function definition for ‘model.matrix’ multinbmod: no visible binding for global variable ‘contrasts’ multinbmod: no visible global function definition for ‘model.response’ multinbmod: no visible global function definition for ‘model.offset’ Undefined global functions or variables: contrasts is.empty.model model.matrix model.offset model.response nlminb Consider adding importFrom("stats", "contrasts", "is.empty.model", "model.matrix", "model.offset", "model.response", "nlminb") to your NAMESPACE file. * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd line widths ... OK * checking Rd cross-references ... OK * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK * checking Rd \usage sections ... OK * checking Rd contents ... OK * checking for unstated dependencies in examples ... OK * checking examples ... ERROR Running examples in ‘multinbmod-Ex.R’ failed The error most likely occurred in: > ### Name: multinbmod-package > ### Title: Regression analysis of overdispersed correlated count data > ### Aliases: multinbmod-package > ### Keywords: multivariate negative binomial model, overdispersed > ### correlated count data > > ### ** Examples > > id <- factor(rep(1:20, rep(5, 20))) > y <- rnbinom(100, mu = rexp(100,1)+rep(rexp(20,.3),rep(5,20)),size=2.5) > x<-rbinom(100,1,.5) > dat <- data.frame(y = y, x = x, id = id) > multinbmod(y~x,data=dat,id=id) Warning in model.matrix.default(mt, mf, contrasts) : non-list contrasts argument ignored ----------- FAILURE REPORT -------------- --- failure: length > 1 in coercion to logical --- --- srcref --- : --- package (from environment) --- multinbmod --- call from context --- multinb.fit(y = Y, x = X, offset, id = id, start.par = c(start.coef, start.phi), control) --- call from argument --- offset != 1 && length(offset) != NROW(y) --- R stacktrace --- where 1: multinb.fit(y = Y, x = X, offset, id = id, start.par = c(start.coef, start.phi), control) where 2: multinbmod(y ~ x, data = dat, id = id) --- value of length: 100 type: logical --- [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [49] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [97] FALSE FALSE FALSE FALSE --- function from context --- function (y, x, offset = 1, id, start.par, control = list()) { if (!is.list(control)) stop("control must be a list") np <- NCOL(x) + 1 if (missing(start.par)) { start.par <- numeric(np) start.par[1] <- log(mean(y + 0.5)) start.par[np] <- 0.5 } else { if (length(start.par) != np) stop("start.par has wrong length") if (start.par[np] <= 0) stop("starting value of phi must be positive") } no.id <- (missing(id) || is.null(id)) if (no.id) { warning("id variable is missing!") return(NULL) } if (offset != 1 && length(offset) != NROW(y)) stop(paste("Number of offsets is", length(offset), ", should equal", NROW(y), "(number of observations)")) tapplys.fun <- function(x, y) { res <- tapply(x, y, sum) return(res[order(unique(y))]) } ydot <- tapplys.fun(y, id) mnbloglik <- function(param, y, x, id, ydot) { p <- length(param) beta <- param[1:(p - 1)] a <- param[p] mu <- offset * exp(x %*% beta) mudot <- tapplys.fun(mu, id) minlog <- -1 * sum(lgamma(a^(-1) + ydot) - lgamma(a^(-1)) - (a^(-1) + ydot) * log(1 + mudot * a) + ydot * log(a) + ydot * log(mudot)) - sum(tapplys.fun(y * log(mu), id) - log(mudot) * ydot) return(minlog) } grr <- function(param, y, x, id, ydot) { p <- length(param) beta <- param[1:(p - 1)] a <- param[p] mu <- offset * exp(x %*% beta) mudot <- tapplys.fun(mu, id) derivbeta <- rep(0, p) for (i in 1:(p - 1)) { derivbeta[i] <- -1 * sum(tapplys.fun(y * x[, i], id) - tapplys.fun(x[, i] * mu, id) * (1 + ydot * a)/(1 + mudot * a)) } theta <- 1/a derivbeta[p] <- (1/a^2) * sum(digamma(theta + ydot) - log(1 + mudot/theta) + (theta + ydot) * mudot/(theta * (theta + mudot)) - ydot/theta - digamma(theta)) return(derivbeta) } hes <- function(param, y, x, id, ydot) { p <- length(param) beta <- param[1:(p - 1)] a <- param[p] mu <- offset * exp(x %*% beta) mudot <- tapplys.fun(mu, id) hesmat <- matrix(numeric(p^2), p, p) for (j in (1:(p - 1))) { for (i in (1:(p - 1))) { hesmat[i, j] <- -1 * sum(-1 * tapplys.fun(mu * x[, i] * x[, j], id) * (1 + a * ydot)/(1 + a * mudot) + tapplys.fun(mu * x[, i], id) * tapplys.fun(mu * x[, j], id) * a * (1 + a * ydot)/(1 + a * mudot)^2) } } for (i in (1:(p - 1))) { hesmat[i, p] <- hesmat[p, i] <- sum(tapplys.fun(mu * x[, i], id) * (ydot - mudot)/(1 + mudot * a)^2) } hesmat[p, p] <- -1 * sum(-ydot/a^2 + trigamma(a^(-1) + ydot)/a^4 + 2 * digamma(a^(-1) + ydot)/a^3 - trigamma(a^(-1))/a^4 - 2 * digamma(a^(-1))/a^3 + 2 * mudot/(a^2 * (1 + mudot * a)) - 2 * log(1 + mudot * a)/a^3 + (a^(-1) + ydot) * mudot^2/(1 + mudot * a)^2) return(hesmat) } res <- nlminb(start.par, mnbloglik, grr, lower = c(rep(-Inf, (np - 1)), 0.01), y = y, x = x, id = id, ydot = ydot, control = control) betamle <- res$par[1:(np - 1)] names(betamle) <- names(as.data.frame(x)) amle <- as.numeric(res$par[np]) muest <- offset * exp(x %*% betamle) mudotest <- tapplys.fun(muest, id) svec <- (y - muest) Nid <- length(unique(id)) var0 <- var1 <- matrix(numeric((np - 1)^2), (np - 1)) for (i in 1:Nid) { eks <- x[(pick <- id == unique(id)[i]), , drop = F] eks <- as.matrix(eks) mus <- muest[pick] ders <- eks * mus if (length(mus) > 1) { visinv <- diag(1/mus) - (amle/(1 + amle * sum(mus))) * rep(1, length(mus)) %*% t(rep(1, length(mus))) var0 <- var0 + t(ders) %*% visinv %*% ders vsan <- svec[pick] %*% t(svec[pick]) var1 <- var1 + t(ders) %*% visinv %*% vsan %*% visinv %*% ders } else { visinv <- 1/(mus * (1 + amle * mus)) var0 <- var0 + visinv * (t(ders) %*% ders) var1 <- var1 + (t(ders) %*% ders) * (visinv^2 * (svec[pick])^2) } } modcov <- solve(var0) robcovbetamle <- solve(var0) %*% var1 %*% solve(var0) robsebetamle <- sqrt(diag(robcovbetamle)) names(robsebetamle) <- names(betamle) tvalues <- betamle/robsebetamle modsebemle <- sqrt(diag(modcov)) names(modsebemle) <- names(betamle) amle <- res$par[np] varesta <- 1/(-1 * sum(-ydot/amle^2 + trigamma(amle^(-1) + ydot)/amle^4 + 2 * digamma(amle^(-1) + ydot)/amle^3 - trigamma(amle^(-1))/amle^4 - 2 * digamma(amle^(-1))/amle^3 + 2 * mudotest/(amle^2 * (1 + mudotest * amle)) - 2 * log(1 + mudotest * amle)/amle^3 + (amle^(-1) + ydot) * mudotest^2/(1 + mudotest * amle)^2)) myres <- list(betamle, modsebemle, robsebetamle, tvalues, modcov, robcovbetamle, amle, sqrt(varesta), 2 * (res$objective + sum(lgamma(y + 1))), ifelse(res$convergence == 0, TRUE, FALSE), res$iterations) names(myres) <- c("estimated regression coefficients", "se from model", "robust se", "robust t-values", "covariance of beta estimates from model", "robust covariance of beta estimates", "estimated phi", "se(phi)", "-2 x loglikelihood", "converged?", "iterations") myres } --- function search by body --- Function multinb.fit in namespace multinbmod has this body. ----------- END OF FAILURE REPORT -------------- Fatal error: length > 1 in coercion to logical * checking PDF version of manual ... OK * checking for non-standard things in the check directory ... OK * checking for detritus in the temp directory ... OK * DONE Status: 1 ERROR, 1 NOTE See ‘/data/gannet/ripley/R/packages/tests-LENGTH1/multinbmod.Rcheck/00check.log’ for details. Command exited with non-zero status 1 Time 0:44.97, 25.18 + 5.16