Last updated on 2020-02-10 07:51:44 CET.
| Flavor | Version | Tinstall | Tcheck | Ttotal | Status | Flags |
|---|---|---|---|---|---|---|
| r-devel-linux-x86_64-debian-clang | 2.0.2 | 11.76 | 37.04 | 48.80 | ERROR | |
| r-devel-linux-x86_64-debian-gcc | 2.0.2 | 10.09 | 29.28 | 39.37 | ERROR | |
| r-devel-linux-x86_64-fedora-clang | 2.0.2 | 63.99 | ERROR | |||
| r-devel-linux-x86_64-fedora-gcc | 2.0.2 | 62.71 | ERROR | |||
| r-devel-windows-ix86+x86_64 | 2.0.2 | 40.00 | 64.00 | 104.00 | OK | |
| r-devel-windows-ix86+x86_64-gcc8 | 2.0.2 | 42.00 | 91.00 | 133.00 | OK | |
| r-patched-linux-x86_64 | 2.0.2 | 10.83 | 33.38 | 44.21 | OK | |
| r-patched-solaris-x86 | 2.0.2 | 84.40 | ERROR | |||
| r-release-linux-x86_64 | 2.0.2 | 10.32 | 34.19 | 44.51 | OK | |
| r-release-windows-ix86+x86_64 | 2.0.2 | 32.00 | 87.00 | 119.00 | OK | |
| r-release-osx-x86_64 | 2.0.2 | OK | ||||
| r-oldrel-windows-ix86+x86_64 | 2.0.2 | 31.00 | 53.00 | 84.00 | OK | |
| r-oldrel-osx-x86_64 | 2.0.2 | OK |
Version: 2.0.2
Check: examples
Result: ERROR
Running examples in 'norm2-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: impNorm
> ### Title: Imputation and prediction for incomplete multivariate normal
> ### data
> ### Aliases: impNorm impNorm.default impNorm.formula impNorm.norm
> ### Keywords: multivariate NA
>
> ### ** Examples
>
>
> ## run EM for marijuana data with ridge prior
> data(marijuana)
> emResult <- emNorm(marijuana, prior="ridge", prior.df=0.5)
>
> ## generate 25 multiple imputations by running 25 chains
> ## of 100 iterations each, starting each chain at the
> ## posterior mode
> set.seed(456)
> imp.list <- as.list(NULL)
> for(m in 1:25){
+ mcmcResult <- mcmcNorm(emResult, iter=100)
+ imp.list[[m]] <- impNorm(mcmcResult)}
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
where 2: mcmcNorm.norm(emResult, iter = 100)
where 3: mcmcNorm(emResult, iter = 100)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, starting.values, iter = 1000,
multicycle = NULL, seeds = NULL, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, save.all.series = TRUE, save.worst.series = FALSE,
worst.linear.coef = NULL, impute.every = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len <- as.integer(255)
msg <- format("", width = msg.len)
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
length(prior.df) <- 1L
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1L
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter)) {
length(iter) <- 1L
}
storage.mode(iter) <- "integer"
if (iter < 1L)
stop("Argument iter must be positive.")
if (!is.null(multicycle)) {
length(multicycle) <- 1L
storage.mode(multicycle) <- "integer"
}
else {
multicycle <- 1L
}
if (multicycle < 1)
stop("Argument multicycle must be positive.")
if (!is.null(seeds)) {
if (length(seeds) != 2L)
stop("Two integer seeds were expected.")
}
else {
seeds = ceiling(runif(2) * .Machine$integer.max)
}
storage.mode(seeds) <- "integer"
if (is.null(impute.every)) {
impute.every <- as.integer(0)
}
else {
impute.every <- as.integer(impute.every)
length(impute.every) <- 1L
if (impute.every < 0L)
stop("Argument impute.every cannot be negative.")
if (impute.every > iter)
stop("Argument impute.every cannot exceed iter.")
}
if (impute.every == 0) {
nimps <- as.integer(0)
}
else {
nimps <- as.integer(floor(iter/impute.every))
}
if (mode(save.all.series) != "logical")
stop("Argument save.all.series should be of mode \"logical\".")
length(save.all.series) <- 1L
if (mode(save.worst.series) != "logical")
stop("Argument save.worst.series should be of mode \"logical\".")
length(save.worst.series) <- 1L
if (save.worst.series) {
if (is.null(worst.linear.coef)) {
stop("If save.worst.series, worst.linear.coef must be provided.")
}
}
if (is.null(worst.linear.coef))
worst.linear.coef <- rep(0, ncol(x) * ncol(y) + (ncol(y) *
(ncol(y) + 1L))/2L)
storage.mode(worst.linear.coef) <- "double"
y.imp <- y
if (save.all.series) {
series.length <- as.integer(iter)
series.beta <- array(0, c(ncol(x), ncol(y), iter))
dimnames(series.beta) <- list(rownames(beta), colnames(beta),
NULL)
series.sigma <- array(0, c(ncol(y), ncol(y), iter))
dimnames(series.sigma) <- list(rownames(sigma), colnames(sigma),
NULL)
}
else {
series.length <- as.integer(0)
series.beta <- array(0, c(ncol(x), ncol(y), 0))
series.sigma <- array(0, c(ncol(y), ncol(y), 0))
}
storage.mode(series.beta) <- "double"
storage.mode(series.sigma) <- "double"
imp.list <- array(0, c(nrow(y), ncol(y), nimps))
storage.mode(imp.list) <- "double"
series.worst <- numeric(iter)
storage.mode(series.worst) <- "double"
tmp <- .Fortran("norm_mcmc", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, impute.every = impute.every,
nimps = nimps, save.all.series = save.all.series, save.worst.series = save.worst.series,
worst.linear.coef = worst.linear.coef, series.length = series.length,
beta = beta, sigma = sigma, y.imp = y.imp, iter.actual = integer(1),
series.beta = series.beta, series.sigma = series.sigma,
loglik = numeric(iter), logpost = numeric(iter), series.worst = series.worst,
npatt = integer(1), mis = matrix(logical(1), nrow(y),
ncol(y)), n.in.patt = integer(nrow(y)), n.obs = integer(ncol(y)),
which.patt = integer(nrow(y)), ybar = numeric(ncol(y)),
ysdv = numeric(ncol(y)), imp.list = imp.list, nimps.actual = integer(1),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
iter <- tmp$iter.actual
if (iter >= 1) {
loglik <- tmp$loglik[1:iter]
logpost <- tmp$logpost[1:iter]
}
else {
loglik <- NULL
logpost <- NULL
}
nimps <- tmp$nimps.actual
if ((impute.every == 0) | (nimps == 0)) {
imp.list <- NULL
impute.every <- NULL
}
else {
imp.list <- as.list(1:nimps)
for (i in 1:nimps) {
imp.list[[i]] <- tmp$imp.list[, , i]
dimnames(imp.list[[i]]) <- list(rownames(y),
colnames(y))
}
}
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (save.all.series & (iter > 0)) {
beta.array <- tmp$series.beta[, , 1:iter, drop = FALSE]
series.beta <- matrix(numeric(1), iter, ncol(x) *
ncol(y))
beta.names <- character(ncol(x) * ncol(y))
posn <- 0
for (j in 1:ncol(x)) {
for (k in 1:ncol(y)) {
posn <- posn + 1
series.beta[, posn] <- beta.array[j, k, ]
beta.names[posn] <- paste(colnames(x)[j], ",",
colnames(y)[k], sep = "")
}
}
colnames(series.beta) <- beta.names
series.beta <- as.ts(series.beta)
sigma.array <- tmp$series.sigma[, , 1:iter, drop = FALSE]
series.sigma <- matrix(numeric(1), iter, ncol(y) *
(ncol(y) + 1)/2)
sigma.names <- character(ncol(y) * (ncol(y) + 1)/2)
posn <- 0
for (k in 1:ncol(y)) {
for (j in k:ncol(y)) {
posn <- posn + 1
series.sigma[, posn] <- sigma.array[j, k, ]
sigma.names[posn] <- paste(colnames(y)[j],
",", colnames(y)[k], sep = "")
}
}
colnames(series.sigma) <- sigma.names
series.sigma <- as.ts(series.sigma)
}
else {
series.beta <- NULL
series.sigma <- NULL
}
if (save.worst.series & (iter > 0)) {
series.worst <- as.ts(tmp$series.worst[1:iter])
}
else {
series.worst <- NULL
worst.linear.coef <- NULL
}
starting.values <- list(beta = beta, sigma = sigma)
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
result <- list(y = y, x = x, method = "MCMC", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, starting.values = starting.values,
param = list(beta = tmp$beta, sigma = tmp$sigma),
loglik = loglik, logpost = logpost, series.worst = series.worst,
series.beta = series.beta, series.sigma = series.sigma,
y.imp = tmp$y.imp, impute.every = impute.every, imp.list = imp.list,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, worst.linear.coef = worst.linear.coef,
ybar = ybar, ysdv = ysdv, msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x28e2708>
<environment: namespace:norm2>
--- function search by body ---
Function mcmcNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: mcmcNorm -> mcmcNorm.norm -> mcmcNorm.default
Execution halted
Flavor: r-devel-linux-x86_64-debian-clang
Version: 2.0.2
Check: tests
Result: ERROR
Running 'emTests.R' [1s/1s]
Running 'impTests.R' [1s/1s]
Running 'logpostTests.R' [1s/1s]
Running 'mcmcTests.R' [1s/1s]
Running 'miInferenceTests.R' [1s/2s]
Running the tests in 'tests/emTests.R' failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with no missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(6), Y2=rnorm(6), Y3=rnorm(6), X1=rnorm(6) )
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x2d87c50>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in 'tests/impTests.R' failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x30911b0>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in 'tests/logpostTests.R' failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x2f611b0>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in 'tests/mcmcTests.R' failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x44631b0>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in 'tests/miInferenceTests.R' failed.
Complete output:
> library(norm2)
>
> # generate M=25 imputations for cholesterol data
> data(cholesterol)
> emResult <- emNorm(cholesterol)
> set.seed(999)
> mcmcResult <- mcmcNorm(emResult, iter=2500, impute.every=100)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
where 2: mcmcNorm.norm(emResult, iter = 2500, impute.every = 100)
where 3: mcmcNorm(emResult, iter = 2500, impute.every = 100)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, starting.values, iter = 1000,
multicycle = NULL, seeds = NULL, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, save.all.series = TRUE, save.worst.series = FALSE,
worst.linear.coef = NULL, impute.every = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len <- as.integer(255)
msg <- format("", width = msg.len)
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
length(prior.df) <- 1L
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1L
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter)) {
length(iter) <- 1L
}
storage.mode(iter) <- "integer"
if (iter < 1L)
stop("Argument iter must be positive.")
if (!is.null(multicycle)) {
length(multicycle) <- 1L
storage.mode(multicycle) <- "integer"
}
else {
multicycle <- 1L
}
if (multicycle < 1)
stop("Argument multicycle must be positive.")
if (!is.null(seeds)) {
if (length(seeds) != 2L)
stop("Two integer seeds were expected.")
}
else {
seeds = ceiling(runif(2) * .Machine$integer.max)
}
storage.mode(seeds) <- "integer"
if (is.null(impute.every)) {
impute.every <- as.integer(0)
}
else {
impute.every <- as.integer(impute.every)
length(impute.every) <- 1L
if (impute.every < 0L)
stop("Argument impute.every cannot be negative.")
if (impute.every > iter)
stop("Argument impute.every cannot exceed iter.")
}
if (impute.every == 0) {
nimps <- as.integer(0)
}
else {
nimps <- as.integer(floor(iter/impute.every))
}
if (mode(save.all.series) != "logical")
stop("Argument save.all.series should be of mode \"logical\".")
length(save.all.series) <- 1L
if (mode(save.worst.series) != "logical")
stop("Argument save.worst.series should be of mode \"logical\".")
length(save.worst.series) <- 1L
if (save.worst.series) {
if (is.null(worst.linear.coef)) {
stop("If save.worst.series, worst.linear.coef must be provided.")
}
}
if (is.null(worst.linear.coef))
worst.linear.coef <- rep(0, ncol(x) * ncol(y) + (ncol(y) *
(ncol(y) + 1L))/2L)
storage.mode(worst.linear.coef) <- "double"
y.imp <- y
if (save.all.series) {
series.length <- as.integer(iter)
series.beta <- array(0, c(ncol(x), ncol(y), iter))
dimnames(series.beta) <- list(rownames(beta), colnames(beta),
NULL)
series.sigma <- array(0, c(ncol(y), ncol(y), iter))
dimnames(series.sigma) <- list(rownames(sigma), colnames(sigma),
NULL)
}
else {
series.length <- as.integer(0)
series.beta <- array(0, c(ncol(x), ncol(y), 0))
series.sigma <- array(0, c(ncol(y), ncol(y), 0))
}
storage.mode(series.beta) <- "double"
storage.mode(series.sigma) <- "double"
imp.list <- array(0, c(nrow(y), ncol(y), nimps))
storage.mode(imp.list) <- "double"
series.worst <- numeric(iter)
storage.mode(series.worst) <- "double"
tmp <- .Fortran("norm_mcmc", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, impute.every = impute.every,
nimps = nimps, save.all.series = save.all.series, save.worst.series = save.worst.series,
worst.linear.coef = worst.linear.coef, series.length = series.length,
beta = beta, sigma = sigma, y.imp = y.imp, iter.actual = integer(1),
series.beta = series.beta, series.sigma = series.sigma,
loglik = numeric(iter), logpost = numeric(iter), series.worst = series.worst,
npatt = integer(1), mis = matrix(logical(1), nrow(y),
ncol(y)), n.in.patt = integer(nrow(y)), n.obs = integer(ncol(y)),
which.patt = integer(nrow(y)), ybar = numeric(ncol(y)),
ysdv = numeric(ncol(y)), imp.list = imp.list, nimps.actual = integer(1),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
iter <- tmp$iter.actual
if (iter >= 1) {
loglik <- tmp$loglik[1:iter]
logpost <- tmp$logpost[1:iter]
}
else {
loglik <- NULL
logpost <- NULL
}
nimps <- tmp$nimps.actual
if ((impute.every == 0) | (nimps == 0)) {
imp.list <- NULL
impute.every <- NULL
}
else {
imp.list <- as.list(1:nimps)
for (i in 1:nimps) {
imp.list[[i]] <- tmp$imp.list[, , i]
dimnames(imp.list[[i]]) <- list(rownames(y),
colnames(y))
}
}
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (save.all.series & (iter > 0)) {
beta.array <- tmp$series.beta[, , 1:iter, drop = FALSE]
series.beta <- matrix(numeric(1), iter, ncol(x) *
ncol(y))
beta.names <- character(ncol(x) * ncol(y))
posn <- 0
for (j in 1:ncol(x)) {
for (k in 1:ncol(y)) {
posn <- posn + 1
series.beta[, posn] <- beta.array[j, k, ]
beta.names[posn] <- paste(colnames(x)[j], ",",
colnames(y)[k], sep = "")
}
}
colnames(series.beta) <- beta.names
series.beta <- as.ts(series.beta)
sigma.array <- tmp$series.sigma[, , 1:iter, drop = FALSE]
series.sigma <- matrix(numeric(1), iter, ncol(y) *
(ncol(y) + 1)/2)
sigma.names <- character(ncol(y) * (ncol(y) + 1)/2)
posn <- 0
for (k in 1:ncol(y)) {
for (j in k:ncol(y)) {
posn <- posn + 1
series.sigma[, posn] <- sigma.array[j, k, ]
sigma.names[posn] <- paste(colnames(y)[j],
",", colnames(y)[k], sep = "")
}
}
colnames(series.sigma) <- sigma.names
series.sigma <- as.ts(series.sigma)
}
else {
series.beta <- NULL
series.sigma <- NULL
}
if (save.worst.series & (iter > 0)) {
series.worst <- as.ts(tmp$series.worst[1:iter])
}
else {
series.worst <- NULL
worst.linear.coef <- NULL
}
starting.values <- list(beta = beta, sigma = sigma)
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
result <- list(y = y, x = x, method = "MCMC", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, starting.values = starting.values,
param = list(beta = tmp$beta, sigma = tmp$sigma),
loglik = loglik, logpost = logpost, series.worst = series.worst,
series.beta = series.beta, series.sigma = series.sigma,
y.imp = tmp$y.imp, impute.every = impute.every, imp.list = imp.list,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, worst.linear.coef = worst.linear.coef,
ybar = ybar, ysdv = ysdv, msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x48c0db0>
<environment: namespace:norm2>
--- function search by body ---
Function mcmcNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: mcmcNorm -> mcmcNorm.norm -> mcmcNorm.default
Execution halted
Flavor: r-devel-linux-x86_64-debian-clang
Version: 2.0.2
Check: examples
Result: ERROR
Running examples in ‘norm2-Ex.R’ failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: impNorm
> ### Title: Imputation and prediction for incomplete multivariate normal
> ### data
> ### Aliases: impNorm impNorm.default impNorm.formula impNorm.norm
> ### Keywords: multivariate NA
>
> ### ** Examples
>
>
> ## run EM for marijuana data with ridge prior
> data(marijuana)
> emResult <- emNorm(marijuana, prior="ridge", prior.df=0.5)
>
> ## generate 25 multiple imputations by running 25 chains
> ## of 100 iterations each, starting each chain at the
> ## posterior mode
> set.seed(456)
> imp.list <- as.list(NULL)
> for(m in 1:25){
+ mcmcResult <- mcmcNorm(emResult, iter=100)
+ imp.list[[m]] <- impNorm(mcmcResult)}
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
where 2: mcmcNorm.norm(emResult, iter = 100)
where 3: mcmcNorm(emResult, iter = 100)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, starting.values, iter = 1000,
multicycle = NULL, seeds = NULL, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, save.all.series = TRUE, save.worst.series = FALSE,
worst.linear.coef = NULL, impute.every = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len <- as.integer(255)
msg <- format("", width = msg.len)
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
length(prior.df) <- 1L
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1L
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter)) {
length(iter) <- 1L
}
storage.mode(iter) <- "integer"
if (iter < 1L)
stop("Argument iter must be positive.")
if (!is.null(multicycle)) {
length(multicycle) <- 1L
storage.mode(multicycle) <- "integer"
}
else {
multicycle <- 1L
}
if (multicycle < 1)
stop("Argument multicycle must be positive.")
if (!is.null(seeds)) {
if (length(seeds) != 2L)
stop("Two integer seeds were expected.")
}
else {
seeds = ceiling(runif(2) * .Machine$integer.max)
}
storage.mode(seeds) <- "integer"
if (is.null(impute.every)) {
impute.every <- as.integer(0)
}
else {
impute.every <- as.integer(impute.every)
length(impute.every) <- 1L
if (impute.every < 0L)
stop("Argument impute.every cannot be negative.")
if (impute.every > iter)
stop("Argument impute.every cannot exceed iter.")
}
if (impute.every == 0) {
nimps <- as.integer(0)
}
else {
nimps <- as.integer(floor(iter/impute.every))
}
if (mode(save.all.series) != "logical")
stop("Argument save.all.series should be of mode \"logical\".")
length(save.all.series) <- 1L
if (mode(save.worst.series) != "logical")
stop("Argument save.worst.series should be of mode \"logical\".")
length(save.worst.series) <- 1L
if (save.worst.series) {
if (is.null(worst.linear.coef)) {
stop("If save.worst.series, worst.linear.coef must be provided.")
}
}
if (is.null(worst.linear.coef))
worst.linear.coef <- rep(0, ncol(x) * ncol(y) + (ncol(y) *
(ncol(y) + 1L))/2L)
storage.mode(worst.linear.coef) <- "double"
y.imp <- y
if (save.all.series) {
series.length <- as.integer(iter)
series.beta <- array(0, c(ncol(x), ncol(y), iter))
dimnames(series.beta) <- list(rownames(beta), colnames(beta),
NULL)
series.sigma <- array(0, c(ncol(y), ncol(y), iter))
dimnames(series.sigma) <- list(rownames(sigma), colnames(sigma),
NULL)
}
else {
series.length <- as.integer(0)
series.beta <- array(0, c(ncol(x), ncol(y), 0))
series.sigma <- array(0, c(ncol(y), ncol(y), 0))
}
storage.mode(series.beta) <- "double"
storage.mode(series.sigma) <- "double"
imp.list <- array(0, c(nrow(y), ncol(y), nimps))
storage.mode(imp.list) <- "double"
series.worst <- numeric(iter)
storage.mode(series.worst) <- "double"
tmp <- .Fortran("norm_mcmc", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, impute.every = impute.every,
nimps = nimps, save.all.series = save.all.series, save.worst.series = save.worst.series,
worst.linear.coef = worst.linear.coef, series.length = series.length,
beta = beta, sigma = sigma, y.imp = y.imp, iter.actual = integer(1),
series.beta = series.beta, series.sigma = series.sigma,
loglik = numeric(iter), logpost = numeric(iter), series.worst = series.worst,
npatt = integer(1), mis = matrix(logical(1), nrow(y),
ncol(y)), n.in.patt = integer(nrow(y)), n.obs = integer(ncol(y)),
which.patt = integer(nrow(y)), ybar = numeric(ncol(y)),
ysdv = numeric(ncol(y)), imp.list = imp.list, nimps.actual = integer(1),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
iter <- tmp$iter.actual
if (iter >= 1) {
loglik <- tmp$loglik[1:iter]
logpost <- tmp$logpost[1:iter]
}
else {
loglik <- NULL
logpost <- NULL
}
nimps <- tmp$nimps.actual
if ((impute.every == 0) | (nimps == 0)) {
imp.list <- NULL
impute.every <- NULL
}
else {
imp.list <- as.list(1:nimps)
for (i in 1:nimps) {
imp.list[[i]] <- tmp$imp.list[, , i]
dimnames(imp.list[[i]]) <- list(rownames(y),
colnames(y))
}
}
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (save.all.series & (iter > 0)) {
beta.array <- tmp$series.beta[, , 1:iter, drop = FALSE]
series.beta <- matrix(numeric(1), iter, ncol(x) *
ncol(y))
beta.names <- character(ncol(x) * ncol(y))
posn <- 0
for (j in 1:ncol(x)) {
for (k in 1:ncol(y)) {
posn <- posn + 1
series.beta[, posn] <- beta.array[j, k, ]
beta.names[posn] <- paste(colnames(x)[j], ",",
colnames(y)[k], sep = "")
}
}
colnames(series.beta) <- beta.names
series.beta <- as.ts(series.beta)
sigma.array <- tmp$series.sigma[, , 1:iter, drop = FALSE]
series.sigma <- matrix(numeric(1), iter, ncol(y) *
(ncol(y) + 1)/2)
sigma.names <- character(ncol(y) * (ncol(y) + 1)/2)
posn <- 0
for (k in 1:ncol(y)) {
for (j in k:ncol(y)) {
posn <- posn + 1
series.sigma[, posn] <- sigma.array[j, k, ]
sigma.names[posn] <- paste(colnames(y)[j],
",", colnames(y)[k], sep = "")
}
}
colnames(series.sigma) <- sigma.names
series.sigma <- as.ts(series.sigma)
}
else {
series.beta <- NULL
series.sigma <- NULL
}
if (save.worst.series & (iter > 0)) {
series.worst <- as.ts(tmp$series.worst[1:iter])
}
else {
series.worst <- NULL
worst.linear.coef <- NULL
}
starting.values <- list(beta = beta, sigma = sigma)
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
result <- list(y = y, x = x, method = "MCMC", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, starting.values = starting.values,
param = list(beta = tmp$beta, sigma = tmp$sigma),
loglik = loglik, logpost = logpost, series.worst = series.worst,
series.beta = series.beta, series.sigma = series.sigma,
y.imp = tmp$y.imp, impute.every = impute.every, imp.list = imp.list,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, worst.linear.coef = worst.linear.coef,
ybar = ybar, ysdv = ysdv, msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x562cea390600>
<environment: namespace:norm2>
--- function search by body ---
Function mcmcNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: mcmcNorm -> mcmcNorm.norm -> mcmcNorm.default
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 2.0.2
Check: tests
Result: ERROR
Running ‘emTests.R’ [1s/1s]
Running ‘impTests.R’ [1s/1s]
Running ‘logpostTests.R’ [1s/1s]
Running ‘mcmcTests.R’ [1s/1s]
Running ‘miInferenceTests.R’ [1s/2s]
Running the tests in ‘tests/emTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with no missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(6), Y2=rnorm(6), Y3=rnorm(6), X1=rnorm(6) )
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x55a0f87c8998>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/impTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x557b79589e18>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/logpostTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x556cd9c8fe18>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/mcmcTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x55d9c6908e18>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/miInferenceTests.R’ failed.
Complete output:
> library(norm2)
>
> # generate M=25 imputations for cholesterol data
> data(cholesterol)
> emResult <- emNorm(cholesterol)
> set.seed(999)
> mcmcResult <- mcmcNorm(emResult, iter=2500, impute.every=100)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
where 2: mcmcNorm.norm(emResult, iter = 2500, impute.every = 100)
where 3: mcmcNorm(emResult, iter = 2500, impute.every = 100)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, starting.values, iter = 1000,
multicycle = NULL, seeds = NULL, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, save.all.series = TRUE, save.worst.series = FALSE,
worst.linear.coef = NULL, impute.every = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len <- as.integer(255)
msg <- format("", width = msg.len)
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
length(prior.df) <- 1L
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1L
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter)) {
length(iter) <- 1L
}
storage.mode(iter) <- "integer"
if (iter < 1L)
stop("Argument iter must be positive.")
if (!is.null(multicycle)) {
length(multicycle) <- 1L
storage.mode(multicycle) <- "integer"
}
else {
multicycle <- 1L
}
if (multicycle < 1)
stop("Argument multicycle must be positive.")
if (!is.null(seeds)) {
if (length(seeds) != 2L)
stop("Two integer seeds were expected.")
}
else {
seeds = ceiling(runif(2) * .Machine$integer.max)
}
storage.mode(seeds) <- "integer"
if (is.null(impute.every)) {
impute.every <- as.integer(0)
}
else {
impute.every <- as.integer(impute.every)
length(impute.every) <- 1L
if (impute.every < 0L)
stop("Argument impute.every cannot be negative.")
if (impute.every > iter)
stop("Argument impute.every cannot exceed iter.")
}
if (impute.every == 0) {
nimps <- as.integer(0)
}
else {
nimps <- as.integer(floor(iter/impute.every))
}
if (mode(save.all.series) != "logical")
stop("Argument save.all.series should be of mode \"logical\".")
length(save.all.series) <- 1L
if (mode(save.worst.series) != "logical")
stop("Argument save.worst.series should be of mode \"logical\".")
length(save.worst.series) <- 1L
if (save.worst.series) {
if (is.null(worst.linear.coef)) {
stop("If save.worst.series, worst.linear.coef must be provided.")
}
}
if (is.null(worst.linear.coef))
worst.linear.coef <- rep(0, ncol(x) * ncol(y) + (ncol(y) *
(ncol(y) + 1L))/2L)
storage.mode(worst.linear.coef) <- "double"
y.imp <- y
if (save.all.series) {
series.length <- as.integer(iter)
series.beta <- array(0, c(ncol(x), ncol(y), iter))
dimnames(series.beta) <- list(rownames(beta), colnames(beta),
NULL)
series.sigma <- array(0, c(ncol(y), ncol(y), iter))
dimnames(series.sigma) <- list(rownames(sigma), colnames(sigma),
NULL)
}
else {
series.length <- as.integer(0)
series.beta <- array(0, c(ncol(x), ncol(y), 0))
series.sigma <- array(0, c(ncol(y), ncol(y), 0))
}
storage.mode(series.beta) <- "double"
storage.mode(series.sigma) <- "double"
imp.list <- array(0, c(nrow(y), ncol(y), nimps))
storage.mode(imp.list) <- "double"
series.worst <- numeric(iter)
storage.mode(series.worst) <- "double"
tmp <- .Fortran("norm_mcmc", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, impute.every = impute.every,
nimps = nimps, save.all.series = save.all.series, save.worst.series = save.worst.series,
worst.linear.coef = worst.linear.coef, series.length = series.length,
beta = beta, sigma = sigma, y.imp = y.imp, iter.actual = integer(1),
series.beta = series.beta, series.sigma = series.sigma,
loglik = numeric(iter), logpost = numeric(iter), series.worst = series.worst,
npatt = integer(1), mis = matrix(logical(1), nrow(y),
ncol(y)), n.in.patt = integer(nrow(y)), n.obs = integer(ncol(y)),
which.patt = integer(nrow(y)), ybar = numeric(ncol(y)),
ysdv = numeric(ncol(y)), imp.list = imp.list, nimps.actual = integer(1),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
iter <- tmp$iter.actual
if (iter >= 1) {
loglik <- tmp$loglik[1:iter]
logpost <- tmp$logpost[1:iter]
}
else {
loglik <- NULL
logpost <- NULL
}
nimps <- tmp$nimps.actual
if ((impute.every == 0) | (nimps == 0)) {
imp.list <- NULL
impute.every <- NULL
}
else {
imp.list <- as.list(1:nimps)
for (i in 1:nimps) {
imp.list[[i]] <- tmp$imp.list[, , i]
dimnames(imp.list[[i]]) <- list(rownames(y),
colnames(y))
}
}
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (save.all.series & (iter > 0)) {
beta.array <- tmp$series.beta[, , 1:iter, drop = FALSE]
series.beta <- matrix(numeric(1), iter, ncol(x) *
ncol(y))
beta.names <- character(ncol(x) * ncol(y))
posn <- 0
for (j in 1:ncol(x)) {
for (k in 1:ncol(y)) {
posn <- posn + 1
series.beta[, posn] <- beta.array[j, k, ]
beta.names[posn] <- paste(colnames(x)[j], ",",
colnames(y)[k], sep = "")
}
}
colnames(series.beta) <- beta.names
series.beta <- as.ts(series.beta)
sigma.array <- tmp$series.sigma[, , 1:iter, drop = FALSE]
series.sigma <- matrix(numeric(1), iter, ncol(y) *
(ncol(y) + 1)/2)
sigma.names <- character(ncol(y) * (ncol(y) + 1)/2)
posn <- 0
for (k in 1:ncol(y)) {
for (j in k:ncol(y)) {
posn <- posn + 1
series.sigma[, posn] <- sigma.array[j, k, ]
sigma.names[posn] <- paste(colnames(y)[j],
",", colnames(y)[k], sep = "")
}
}
colnames(series.sigma) <- sigma.names
series.sigma <- as.ts(series.sigma)
}
else {
series.beta <- NULL
series.sigma <- NULL
}
if (save.worst.series & (iter > 0)) {
series.worst <- as.ts(tmp$series.worst[1:iter])
}
else {
series.worst <- NULL
worst.linear.coef <- NULL
}
starting.values <- list(beta = beta, sigma = sigma)
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
result <- list(y = y, x = x, method = "MCMC", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, starting.values = starting.values,
param = list(beta = tmp$beta, sigma = tmp$sigma),
loglik = loglik, logpost = logpost, series.worst = series.worst,
series.beta = series.beta, series.sigma = series.sigma,
y.imp = tmp$y.imp, impute.every = impute.every, imp.list = imp.list,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, worst.linear.coef = worst.linear.coef,
ybar = ybar, ysdv = ysdv, msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x55dc5daaa198>
<environment: namespace:norm2>
--- function search by body ---
Function mcmcNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: mcmcNorm -> mcmcNorm.norm -> mcmcNorm.default
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 2.0.2
Check: compiled code
Result: NOTE
File ‘norm2/libs/norm2.so’:
Found no calls to: ‘R_registerRoutines’, ‘R_useDynamicSymbols’
It is good practice to register native routines and to disable symbol
search.
See ‘Writing portable packages’ in the ‘Writing R Extensions’ manual.
Flavors: r-devel-linux-x86_64-fedora-clang, r-devel-linux-x86_64-fedora-gcc
Version: 2.0.2
Check: examples
Result: ERROR
Running examples in ‘norm2-Ex.R’ failed
The error most likely occurred in:
> ### Name: impNorm
> ### Title: Imputation and prediction for incomplete multivariate normal
> ### data
> ### Aliases: impNorm impNorm.default impNorm.formula impNorm.norm
> ### Keywords: multivariate NA
>
> ### ** Examples
>
>
> ## run EM for marijuana data with ridge prior
> data(marijuana)
> emResult <- emNorm(marijuana, prior="ridge", prior.df=0.5)
>
> ## generate 25 multiple imputations by running 25 chains
> ## of 100 iterations each, starting each chain at the
> ## posterior mode
> set.seed(456)
> imp.list <- as.list(NULL)
> for(m in 1:25){
+ mcmcResult <- mcmcNorm(emResult, iter=100)
+ imp.list[[m]] <- impNorm(mcmcResult)}
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
where 2: mcmcNorm.norm(emResult, iter = 100)
where 3: mcmcNorm(emResult, iter = 100)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, starting.values, iter = 1000,
multicycle = NULL, seeds = NULL, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, save.all.series = TRUE, save.worst.series = FALSE,
worst.linear.coef = NULL, impute.every = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len <- as.integer(255)
msg <- format("", width = msg.len)
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
length(prior.df) <- 1L
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1L
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter)) {
length(iter) <- 1L
}
storage.mode(iter) <- "integer"
if (iter < 1L)
stop("Argument iter must be positive.")
if (!is.null(multicycle)) {
length(multicycle) <- 1L
storage.mode(multicycle) <- "integer"
}
else {
multicycle <- 1L
}
if (multicycle < 1)
stop("Argument multicycle must be positive.")
if (!is.null(seeds)) {
if (length(seeds) != 2L)
stop("Two integer seeds were expected.")
}
else {
seeds = ceiling(runif(2) * .Machine$integer.max)
}
storage.mode(seeds) <- "integer"
if (is.null(impute.every)) {
impute.every <- as.integer(0)
}
else {
impute.every <- as.integer(impute.every)
length(impute.every) <- 1L
if (impute.every < 0L)
stop("Argument impute.every cannot be negative.")
if (impute.every > iter)
stop("Argument impute.every cannot exceed iter.")
}
if (impute.every == 0) {
nimps <- as.integer(0)
}
else {
nimps <- as.integer(floor(iter/impute.every))
}
if (mode(save.all.series) != "logical")
stop("Argument save.all.series should be of mode \"logical\".")
length(save.all.series) <- 1L
if (mode(save.worst.series) != "logical")
stop("Argument save.worst.series should be of mode \"logical\".")
length(save.worst.series) <- 1L
if (save.worst.series) {
if (is.null(worst.linear.coef)) {
stop("If save.worst.series, worst.linear.coef must be provided.")
}
}
if (is.null(worst.linear.coef))
worst.linear.coef <- rep(0, ncol(x) * ncol(y) + (ncol(y) *
(ncol(y) + 1L))/2L)
storage.mode(worst.linear.coef) <- "double"
y.imp <- y
if (save.all.series) {
series.length <- as.integer(iter)
series.beta <- array(0, c(ncol(x), ncol(y), iter))
dimnames(series.beta) <- list(rownames(beta), colnames(beta),
NULL)
series.sigma <- array(0, c(ncol(y), ncol(y), iter))
dimnames(series.sigma) <- list(rownames(sigma), colnames(sigma),
NULL)
}
else {
series.length <- as.integer(0)
series.beta <- array(0, c(ncol(x), ncol(y), 0))
series.sigma <- array(0, c(ncol(y), ncol(y), 0))
}
storage.mode(series.beta) <- "double"
storage.mode(series.sigma) <- "double"
imp.list <- array(0, c(nrow(y), ncol(y), nimps))
storage.mode(imp.list) <- "double"
series.worst <- numeric(iter)
storage.mode(series.worst) <- "double"
tmp <- .Fortran("norm_mcmc", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, impute.every = impute.every,
nimps = nimps, save.all.series = save.all.series, save.worst.series = save.worst.series,
worst.linear.coef = worst.linear.coef, series.length = series.length,
beta = beta, sigma = sigma, y.imp = y.imp, iter.actual = integer(1),
series.beta = series.beta, series.sigma = series.sigma,
loglik = numeric(iter), logpost = numeric(iter), series.worst = series.worst,
npatt = integer(1), mis = matrix(logical(1), nrow(y),
ncol(y)), n.in.patt = integer(nrow(y)), n.obs = integer(ncol(y)),
which.patt = integer(nrow(y)), ybar = numeric(ncol(y)),
ysdv = numeric(ncol(y)), imp.list = imp.list, nimps.actual = integer(1),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
iter <- tmp$iter.actual
if (iter >= 1) {
loglik <- tmp$loglik[1:iter]
logpost <- tmp$logpost[1:iter]
}
else {
loglik <- NULL
logpost <- NULL
}
nimps <- tmp$nimps.actual
if ((impute.every == 0) | (nimps == 0)) {
imp.list <- NULL
impute.every <- NULL
}
else {
imp.list <- as.list(1:nimps)
for (i in 1:nimps) {
imp.list[[i]] <- tmp$imp.list[, , i]
dimnames(imp.list[[i]]) <- list(rownames(y),
colnames(y))
}
}
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (save.all.series & (iter > 0)) {
beta.array <- tmp$series.beta[, , 1:iter, drop = FALSE]
series.beta <- matrix(numeric(1), iter, ncol(x) *
ncol(y))
beta.names <- character(ncol(x) * ncol(y))
posn <- 0
for (j in 1:ncol(x)) {
for (k in 1:ncol(y)) {
posn <- posn + 1
series.beta[, posn] <- beta.array[j, k, ]
beta.names[posn] <- paste(colnames(x)[j], ",",
colnames(y)[k], sep = "")
}
}
colnames(series.beta) <- beta.names
series.beta <- as.ts(series.beta)
sigma.array <- tmp$series.sigma[, , 1:iter, drop = FALSE]
series.sigma <- matrix(numeric(1), iter, ncol(y) *
(ncol(y) + 1)/2)
sigma.names <- character(ncol(y) * (ncol(y) + 1)/2)
posn <- 0
for (k in 1:ncol(y)) {
for (j in k:ncol(y)) {
posn <- posn + 1
series.sigma[, posn] <- sigma.array[j, k, ]
sigma.names[posn] <- paste(colnames(y)[j],
",", colnames(y)[k], sep = "")
}
}
colnames(series.sigma) <- sigma.names
series.sigma <- as.ts(series.sigma)
}
else {
series.beta <- NULL
series.sigma <- NULL
}
if (save.worst.series & (iter > 0)) {
series.worst <- as.ts(tmp$series.worst[1:iter])
}
else {
series.worst <- NULL
worst.linear.coef <- NULL
}
starting.values <- list(beta = beta, sigma = sigma)
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
result <- list(y = y, x = x, method = "MCMC", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, starting.values = starting.values,
param = list(beta = tmp$beta, sigma = tmp$sigma),
loglik = loglik, logpost = logpost, series.worst = series.worst,
series.beta = series.beta, series.sigma = series.sigma,
y.imp = tmp$y.imp, impute.every = impute.every, imp.list = imp.list,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, worst.linear.coef = worst.linear.coef,
ybar = ybar, ysdv = ysdv, msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x38d8e48>
<environment: namespace:norm2>
--- function search by body ---
Function mcmcNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: mcmcNorm -> mcmcNorm.norm -> mcmcNorm.default
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 2.0.2
Check: tests
Result: ERROR
Running ‘emTests.R’
Running ‘impTests.R’
Running ‘logpostTests.R’
Running ‘mcmcTests.R’
Running ‘miInferenceTests.R’
Running the tests in ‘tests/emTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with no missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(6), Y2=rnorm(6), Y3=rnorm(6), X1=rnorm(6) )
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x3a40f20>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/impTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x4322120>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/logpostTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x357a120>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/mcmcTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x45f6120>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/miInferenceTests.R’ failed.
Complete output:
> library(norm2)
>
> # generate M=25 imputations for cholesterol data
> data(cholesterol)
> emResult <- emNorm(cholesterol)
> set.seed(999)
> mcmcResult <- mcmcNorm(emResult, iter=2500, impute.every=100)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
where 2: mcmcNorm.norm(emResult, iter = 2500, impute.every = 100)
where 3: mcmcNorm(emResult, iter = 2500, impute.every = 100)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, starting.values, iter = 1000,
multicycle = NULL, seeds = NULL, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, save.all.series = TRUE, save.worst.series = FALSE,
worst.linear.coef = NULL, impute.every = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len <- as.integer(255)
msg <- format("", width = msg.len)
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
length(prior.df) <- 1L
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1L
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter)) {
length(iter) <- 1L
}
storage.mode(iter) <- "integer"
if (iter < 1L)
stop("Argument iter must be positive.")
if (!is.null(multicycle)) {
length(multicycle) <- 1L
storage.mode(multicycle) <- "integer"
}
else {
multicycle <- 1L
}
if (multicycle < 1)
stop("Argument multicycle must be positive.")
if (!is.null(seeds)) {
if (length(seeds) != 2L)
stop("Two integer seeds were expected.")
}
else {
seeds = ceiling(runif(2) * .Machine$integer.max)
}
storage.mode(seeds) <- "integer"
if (is.null(impute.every)) {
impute.every <- as.integer(0)
}
else {
impute.every <- as.integer(impute.every)
length(impute.every) <- 1L
if (impute.every < 0L)
stop("Argument impute.every cannot be negative.")
if (impute.every > iter)
stop("Argument impute.every cannot exceed iter.")
}
if (impute.every == 0) {
nimps <- as.integer(0)
}
else {
nimps <- as.integer(floor(iter/impute.every))
}
if (mode(save.all.series) != "logical")
stop("Argument save.all.series should be of mode \"logical\".")
length(save.all.series) <- 1L
if (mode(save.worst.series) != "logical")
stop("Argument save.worst.series should be of mode \"logical\".")
length(save.worst.series) <- 1L
if (save.worst.series) {
if (is.null(worst.linear.coef)) {
stop("If save.worst.series, worst.linear.coef must be provided.")
}
}
if (is.null(worst.linear.coef))
worst.linear.coef <- rep(0, ncol(x) * ncol(y) + (ncol(y) *
(ncol(y) + 1L))/2L)
storage.mode(worst.linear.coef) <- "double"
y.imp <- y
if (save.all.series) {
series.length <- as.integer(iter)
series.beta <- array(0, c(ncol(x), ncol(y), iter))
dimnames(series.beta) <- list(rownames(beta), colnames(beta),
NULL)
series.sigma <- array(0, c(ncol(y), ncol(y), iter))
dimnames(series.sigma) <- list(rownames(sigma), colnames(sigma),
NULL)
}
else {
series.length <- as.integer(0)
series.beta <- array(0, c(ncol(x), ncol(y), 0))
series.sigma <- array(0, c(ncol(y), ncol(y), 0))
}
storage.mode(series.beta) <- "double"
storage.mode(series.sigma) <- "double"
imp.list <- array(0, c(nrow(y), ncol(y), nimps))
storage.mode(imp.list) <- "double"
series.worst <- numeric(iter)
storage.mode(series.worst) <- "double"
tmp <- .Fortran("norm_mcmc", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, impute.every = impute.every,
nimps = nimps, save.all.series = save.all.series, save.worst.series = save.worst.series,
worst.linear.coef = worst.linear.coef, series.length = series.length,
beta = beta, sigma = sigma, y.imp = y.imp, iter.actual = integer(1),
series.beta = series.beta, series.sigma = series.sigma,
loglik = numeric(iter), logpost = numeric(iter), series.worst = series.worst,
npatt = integer(1), mis = matrix(logical(1), nrow(y),
ncol(y)), n.in.patt = integer(nrow(y)), n.obs = integer(ncol(y)),
which.patt = integer(nrow(y)), ybar = numeric(ncol(y)),
ysdv = numeric(ncol(y)), imp.list = imp.list, nimps.actual = integer(1),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
iter <- tmp$iter.actual
if (iter >= 1) {
loglik <- tmp$loglik[1:iter]
logpost <- tmp$logpost[1:iter]
}
else {
loglik <- NULL
logpost <- NULL
}
nimps <- tmp$nimps.actual
if ((impute.every == 0) | (nimps == 0)) {
imp.list <- NULL
impute.every <- NULL
}
else {
imp.list <- as.list(1:nimps)
for (i in 1:nimps) {
imp.list[[i]] <- tmp$imp.list[, , i]
dimnames(imp.list[[i]]) <- list(rownames(y),
colnames(y))
}
}
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (save.all.series & (iter > 0)) {
beta.array <- tmp$series.beta[, , 1:iter, drop = FALSE]
series.beta <- matrix(numeric(1), iter, ncol(x) *
ncol(y))
beta.names <- character(ncol(x) * ncol(y))
posn <- 0
for (j in 1:ncol(x)) {
for (k in 1:ncol(y)) {
posn <- posn + 1
series.beta[, posn] <- beta.array[j, k, ]
beta.names[posn] <- paste(colnames(x)[j], ",",
colnames(y)[k], sep = "")
}
}
colnames(series.beta) <- beta.names
series.beta <- as.ts(series.beta)
sigma.array <- tmp$series.sigma[, , 1:iter, drop = FALSE]
series.sigma <- matrix(numeric(1), iter, ncol(y) *
(ncol(y) + 1)/2)
sigma.names <- character(ncol(y) * (ncol(y) + 1)/2)
posn <- 0
for (k in 1:ncol(y)) {
for (j in k:ncol(y)) {
posn <- posn + 1
series.sigma[, posn] <- sigma.array[j, k, ]
sigma.names[posn] <- paste(colnames(y)[j],
",", colnames(y)[k], sep = "")
}
}
colnames(series.sigma) <- sigma.names
series.sigma <- as.ts(series.sigma)
}
else {
series.beta <- NULL
series.sigma <- NULL
}
if (save.worst.series & (iter > 0)) {
series.worst <- as.ts(tmp$series.worst[1:iter])
}
else {
series.worst <- NULL
worst.linear.coef <- NULL
}
starting.values <- list(beta = beta, sigma = sigma)
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
result <- list(y = y, x = x, method = "MCMC", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, starting.values = starting.values,
param = list(beta = tmp$beta, sigma = tmp$sigma),
loglik = loglik, logpost = logpost, series.worst = series.worst,
series.beta = series.beta, series.sigma = series.sigma,
y.imp = tmp$y.imp, impute.every = impute.every, imp.list = imp.list,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, worst.linear.coef = worst.linear.coef,
ybar = ybar, ysdv = ysdv, msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x3e436b0>
<environment: namespace:norm2>
--- function search by body ---
Function mcmcNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: mcmcNorm -> mcmcNorm.norm -> mcmcNorm.default
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 2.0.2
Check: examples
Result: ERROR
Running examples in ‘norm2-Ex.R’ failed
The error most likely occurred in:
> ### Name: impNorm
> ### Title: Imputation and prediction for incomplete multivariate normal
> ### data
> ### Aliases: impNorm impNorm.default impNorm.formula impNorm.norm
> ### Keywords: multivariate NA
>
> ### ** Examples
>
>
> ## run EM for marijuana data with ridge prior
> data(marijuana)
> emResult <- emNorm(marijuana, prior="ridge", prior.df=0.5)
>
> ## generate 25 multiple imputations by running 25 chains
> ## of 100 iterations each, starting each chain at the
> ## posterior mode
> set.seed(456)
> imp.list <- as.list(NULL)
> for(m in 1:25){
+ mcmcResult <- mcmcNorm(emResult, iter=100)
+ imp.list[[m]] <- impNorm(mcmcResult)}
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
where 2: mcmcNorm.norm(emResult, iter = 100)
where 3: mcmcNorm(emResult, iter = 100)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, starting.values, iter = 1000,
multicycle = NULL, seeds = NULL, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, save.all.series = TRUE, save.worst.series = FALSE,
worst.linear.coef = NULL, impute.every = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len <- as.integer(255)
msg <- format("", width = msg.len)
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
length(prior.df) <- 1L
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1L
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter)) {
length(iter) <- 1L
}
storage.mode(iter) <- "integer"
if (iter < 1L)
stop("Argument iter must be positive.")
if (!is.null(multicycle)) {
length(multicycle) <- 1L
storage.mode(multicycle) <- "integer"
}
else {
multicycle <- 1L
}
if (multicycle < 1)
stop("Argument multicycle must be positive.")
if (!is.null(seeds)) {
if (length(seeds) != 2L)
stop("Two integer seeds were expected.")
}
else {
seeds = ceiling(runif(2) * .Machine$integer.max)
}
storage.mode(seeds) <- "integer"
if (is.null(impute.every)) {
impute.every <- as.integer(0)
}
else {
impute.every <- as.integer(impute.every)
length(impute.every) <- 1L
if (impute.every < 0L)
stop("Argument impute.every cannot be negative.")
if (impute.every > iter)
stop("Argument impute.every cannot exceed iter.")
}
if (impute.every == 0) {
nimps <- as.integer(0)
}
else {
nimps <- as.integer(floor(iter/impute.every))
}
if (mode(save.all.series) != "logical")
stop("Argument save.all.series should be of mode \"logical\".")
length(save.all.series) <- 1L
if (mode(save.worst.series) != "logical")
stop("Argument save.worst.series should be of mode \"logical\".")
length(save.worst.series) <- 1L
if (save.worst.series) {
if (is.null(worst.linear.coef)) {
stop("If save.worst.series, worst.linear.coef must be provided.")
}
}
if (is.null(worst.linear.coef))
worst.linear.coef <- rep(0, ncol(x) * ncol(y) + (ncol(y) *
(ncol(y) + 1L))/2L)
storage.mode(worst.linear.coef) <- "double"
y.imp <- y
if (save.all.series) {
series.length <- as.integer(iter)
series.beta <- array(0, c(ncol(x), ncol(y), iter))
dimnames(series.beta) <- list(rownames(beta), colnames(beta),
NULL)
series.sigma <- array(0, c(ncol(y), ncol(y), iter))
dimnames(series.sigma) <- list(rownames(sigma), colnames(sigma),
NULL)
}
else {
series.length <- as.integer(0)
series.beta <- array(0, c(ncol(x), ncol(y), 0))
series.sigma <- array(0, c(ncol(y), ncol(y), 0))
}
storage.mode(series.beta) <- "double"
storage.mode(series.sigma) <- "double"
imp.list <- array(0, c(nrow(y), ncol(y), nimps))
storage.mode(imp.list) <- "double"
series.worst <- numeric(iter)
storage.mode(series.worst) <- "double"
tmp <- .Fortran("norm_mcmc", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, impute.every = impute.every,
nimps = nimps, save.all.series = save.all.series, save.worst.series = save.worst.series,
worst.linear.coef = worst.linear.coef, series.length = series.length,
beta = beta, sigma = sigma, y.imp = y.imp, iter.actual = integer(1),
series.beta = series.beta, series.sigma = series.sigma,
loglik = numeric(iter), logpost = numeric(iter), series.worst = series.worst,
npatt = integer(1), mis = matrix(logical(1), nrow(y),
ncol(y)), n.in.patt = integer(nrow(y)), n.obs = integer(ncol(y)),
which.patt = integer(nrow(y)), ybar = numeric(ncol(y)),
ysdv = numeric(ncol(y)), imp.list = imp.list, nimps.actual = integer(1),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
iter <- tmp$iter.actual
if (iter >= 1) {
loglik <- tmp$loglik[1:iter]
logpost <- tmp$logpost[1:iter]
}
else {
loglik <- NULL
logpost <- NULL
}
nimps <- tmp$nimps.actual
if ((impute.every == 0) | (nimps == 0)) {
imp.list <- NULL
impute.every <- NULL
}
else {
imp.list <- as.list(1:nimps)
for (i in 1:nimps) {
imp.list[[i]] <- tmp$imp.list[, , i]
dimnames(imp.list[[i]]) <- list(rownames(y),
colnames(y))
}
}
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (save.all.series & (iter > 0)) {
beta.array <- tmp$series.beta[, , 1:iter, drop = FALSE]
series.beta <- matrix(numeric(1), iter, ncol(x) *
ncol(y))
beta.names <- character(ncol(x) * ncol(y))
posn <- 0
for (j in 1:ncol(x)) {
for (k in 1:ncol(y)) {
posn <- posn + 1
series.beta[, posn] <- beta.array[j, k, ]
beta.names[posn] <- paste(colnames(x)[j], ",",
colnames(y)[k], sep = "")
}
}
colnames(series.beta) <- beta.names
series.beta <- as.ts(series.beta)
sigma.array <- tmp$series.sigma[, , 1:iter, drop = FALSE]
series.sigma <- matrix(numeric(1), iter, ncol(y) *
(ncol(y) + 1)/2)
sigma.names <- character(ncol(y) * (ncol(y) + 1)/2)
posn <- 0
for (k in 1:ncol(y)) {
for (j in k:ncol(y)) {
posn <- posn + 1
series.sigma[, posn] <- sigma.array[j, k, ]
sigma.names[posn] <- paste(colnames(y)[j],
",", colnames(y)[k], sep = "")
}
}
colnames(series.sigma) <- sigma.names
series.sigma <- as.ts(series.sigma)
}
else {
series.beta <- NULL
series.sigma <- NULL
}
if (save.worst.series & (iter > 0)) {
series.worst <- as.ts(tmp$series.worst[1:iter])
}
else {
series.worst <- NULL
worst.linear.coef <- NULL
}
starting.values <- list(beta = beta, sigma = sigma)
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
result <- list(y = y, x = x, method = "MCMC", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, starting.values = starting.values,
param = list(beta = tmp$beta, sigma = tmp$sigma),
loglik = loglik, logpost = logpost, series.worst = series.worst,
series.beta = series.beta, series.sigma = series.sigma,
y.imp = tmp$y.imp, impute.every = impute.every, imp.list = imp.list,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, worst.linear.coef = worst.linear.coef,
ybar = ybar, ysdv = ysdv, msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x2f67978>
<environment: namespace:norm2>
--- function search by body ---
Function mcmcNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: mcmcNorm -> mcmcNorm.norm -> mcmcNorm.default
Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc
Version: 2.0.2
Check: tests
Result: ERROR
Running ‘emTests.R’
Running ‘impTests.R’
Running ‘logpostTests.R’
Running ‘mcmcTests.R’
Running ‘miInferenceTests.R’
Running the tests in ‘tests/emTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with no missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(6), Y2=rnorm(6), Y3=rnorm(6), X1=rnorm(6) )
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x2c89610>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/impTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x2fbca90>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/logpostTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x3effa90>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/mcmcTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(10), Y2=rnorm(10), Y3=rnorm(10), X1=rnorm(10) )
> simdata$Y3[7] <- simdata$Y2[3]<- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values,
prior = prior, prior.df = prior.df, prior.sscp = prior.sscp,
...)
where 2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
where 3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, iter.max = 1000, criterion = NULL,
estimate.worst = TRUE, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, starting.values = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
if (is.null(starting.values)) {
beta <- matrix(0, ncol(x), ncol(y))
sigma <- matrix(0, ncol(y), ncol(y))
startval.present <- FALSE
}
else {
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
startval.present <- TRUE
}
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter.max)) {
length(iter.max) <- 1L
}
storage.mode(iter.max) <- "integer"
if (iter.max < 0)
stop("Argument iter.max cannot be negative.")
if (!is.null(criterion)) {
length(criterion) <- 1L
storage.mode(criterion) <- "double"
if (criterion < 0)
stop("Argument criterion cannot be negative.")
}
else {
criterion <- as.double(1e-05)
}
if (is.null(estimate.worst))
estimate.worst <- FALSE
storage.mode(estimate.worst) <- "logical"
length(estimate.worst) <- 1L
nparam <- as.integer(ncol(x) * ncol(y) + (ncol(y) * (ncol(y) +
1L))/2L)
tmp <- .Fortran("norm_em", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter.max = iter.max,
criterion = criterion, estimate.worst = estimate.worst,
startval.present = startval.present, iter = integer(1),
converged = logical(1), rel.diff = numeric(1), loglik = numeric(iter.max),
logpost = numeric(iter.max), beta = beta, sigma = sigma,
y.imp = y, npatt = integer(1), mis = matrix(logical(1),
nrow(y), ncol(y)), n.in.patt = integer(nrow(y)),
n.obs = integer(ncol(y)), which.patt = integer(nrow(y)),
ybar = numeric(ncol(y)), ysdv = numeric(ncol(y)), rate.beta = beta,
rate.sigma = sigma, em.worst.ok = logical(1), worst.frac = numeric(1),
nparam = nparam, worst.linear.coef = numeric(nparam),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
if ((!tmp$converged) & (iter.max > 0)) {
warning(paste("Algorithm did not converge by iteration",
format(tmp$iter)))
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
rel.diff <- tmp$rel.diff
if (tmp$iter == 0)
rel.diff <- NA
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (startval.present) {
starting.values <- list(beta = beta, sigma = sigma)
}
else {
starting.values <- NULL
}
if (tmp$iter >= 1) {
loglik <- tmp$loglik[1:tmp$iter]
logpost <- tmp$logpost[1:tmp$iter]
}
else {
loglik <- NULL
logpost <- NULL
}
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
em.worst.ok <- tmp$em.worst.ok
if (em.worst.ok) {
worst.frac <- tmp$worst.frac
if (worst.frac == 0) {
worst.linear.coef <- NULL
}
else {
worst.linear.coef <- tmp$worst.linear.coef
}
}
else {
worst.frac <- NULL
worst.linear.coef <- NULL
}
result <- list(y = y, x = x, method = "EM", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, starting.values = starting.values,
iter = tmp$iter, converged = tmp$converged, criterion = criterion,
estimate.worst = estimate.worst, loglik = loglik,
logpost = logpost, param = list(beta = tmp$beta,
sigma = tmp$sigma), param.rate = list(beta = tmp$rate.beta,
sigma = tmp$rate.sigma), y.mean.imp = tmp$y.imp,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, rel.diff = rel.diff,
ybar = ybar, ysdv = ysdv, em.worst.ok = em.worst.ok,
worst.frac = worst.frac, worst.linear.coef = worst.linear.coef,
msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x45bea90>
<environment: namespace:norm2>
--- function search by body ---
Function emNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: emNorm -> emNorm.formula -> emNorm.default
Execution halted
Running the tests in ‘tests/miInferenceTests.R’ failed.
Complete output:
> library(norm2)
>
> # generate M=25 imputations for cholesterol data
> data(cholesterol)
> emResult <- emNorm(cholesterol)
> set.seed(999)
> mcmcResult <- mcmcNorm(emResult, iter=2500, impute.every=100)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
norm2
--- call from context ---
mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
--- call from argument ---
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
--- R stacktrace ---
where 1: mcmcNorm.default(obj$y, x = obj$x, intercept = FALSE, starting.values = starting.values,
iter = iter, multicycle = multicycle, seeds = seeds, prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, save.all.series = save.all.series,
save.worst.series = save.worst.series, worst.linear.coef = worst.linear.coef,
impute.every = impute.every, ...)
where 2: mcmcNorm.norm(emResult, iter = 2500, impute.every = 100)
where 3: mcmcNorm(emResult, iter = 2500, impute.every = 100)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (obj, x = NULL, intercept = TRUE, starting.values, iter = 1000,
multicycle = NULL, seeds = NULL, prior = "uniform", prior.df = NULL,
prior.sscp = NULL, save.all.series = TRUE, save.worst.series = FALSE,
worst.linear.coef = NULL, impute.every = NULL, ...)
{
y <- obj
if (class(y) == "data.frame") {
status <- logical(length(y))
for (j in 1:length(y)) status[j] <- is.factor(y[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"y\"", "converted to mode \"numeric\".")
warning(msg)
}
}
y <- data.matrix(y)
storage.mode(y) <- "double"
if (is.null(colnames(y))) {
colnames(y) <- paste("Y", as.character(1:ncol(y)), sep = "")
}
if (is.null(x)) {
x <- matrix(1, nrow(y), 1)
colnames(x) <- "CONST"
}
else {
if (class(x) == "data.frame") {
status <- logical(length(x))
for (j in 1:length(x)) status[j] <- is.factor(x[[j]])
if (any(status)) {
msg <- paste("Factors in argument \"x\"", "converted to mode \"numeric\".")
warning(msg)
}
}
x <- data.matrix(x)
if (nrow(x) != nrow(y)) {
stop("Arguments x and y do not have the same number of rows.")
}
if (any(is.na(x))) {
stop("Missing values in x are not allowed.")
}
if (is.null(colnames(x))) {
colnames(x) <- paste("X", as.character(1:ncol(x)),
sep = "")
}
if (intercept) {
x <- cbind(CONST = 1, x)
}
}
rownames(x) <- rownames(y)
storage.mode(x) <- "double"
mvcode <- .Machine$double.xmax
y[is.na(y)] <- mvcode
msg.len <- as.integer(255)
msg <- format("", width = msg.len)
msg.len.max <- 40L
msg.codes <- matrix(0L, msg.len.max, 8L)
beta <- data.matrix(starting.values$beta)
sigma <- data.matrix(starting.values$sigma)
if ((nrow(beta) != ncol(x)) | (ncol(beta) != ncol(y)))
stop("Incorrect dimensions for argument beta.")
if ((nrow(sigma) != ncol(y)) | (ncol(sigma) != ncol(y)))
stop("Incorrect dimensions for argument sigma.")
storage.mode(beta) <- "double"
storage.mode(sigma) <- "double"
rownames(beta) <- colnames(x)
colnames(beta) <- colnames(y)
rownames(sigma) <- colnames(y)
colnames(sigma) <- colnames(y)
if (prior == "uniform") {
prior.type.int <- 1L
prior.df <- -(ncol(y) + 1)
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "jeffreys") {
prior.type.int <- 2L
prior.df <- 0
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "ridge") {
prior.type.int <- 3L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"ridge\".")
length(prior.df) <- 1L
prior.sscp <- matrix(0, ncol(y), ncol(y))
}
else if (prior == "invwish") {
prior.type.int <- 4L
if (is.null(prior.df))
stop("Argument prior.df must be provided when prior = \"invwish\".")
length(prior.df) <- 1L
if (is.null(prior.sscp))
stop("Argument prior.sscp must be provided when prior = \"invwish\".")
prior.sscp <- data.matrix(prior.sscp)
if ((nrow(prior.sscp) != ncol(y)) | (ncol(prior.sscp) !=
ncol(y)))
stop("Argument prior.sscp has incorrect dimensions.")
}
else {
msg <- paste("Prior type \"", prior, "\" not recognized.",
sep = "")
stop(msg)
}
storage.mode(prior.type.int) <- "integer"
length(prior.type.int) <- 1L
storage.mode(prior.df) <- "double"
length(prior.df) <- 1L
storage.mode(prior.sscp) <- "double"
if (!missing(iter)) {
length(iter) <- 1L
}
storage.mode(iter) <- "integer"
if (iter < 1L)
stop("Argument iter must be positive.")
if (!is.null(multicycle)) {
length(multicycle) <- 1L
storage.mode(multicycle) <- "integer"
}
else {
multicycle <- 1L
}
if (multicycle < 1)
stop("Argument multicycle must be positive.")
if (!is.null(seeds)) {
if (length(seeds) != 2L)
stop("Two integer seeds were expected.")
}
else {
seeds = ceiling(runif(2) * .Machine$integer.max)
}
storage.mode(seeds) <- "integer"
if (is.null(impute.every)) {
impute.every <- as.integer(0)
}
else {
impute.every <- as.integer(impute.every)
length(impute.every) <- 1L
if (impute.every < 0L)
stop("Argument impute.every cannot be negative.")
if (impute.every > iter)
stop("Argument impute.every cannot exceed iter.")
}
if (impute.every == 0) {
nimps <- as.integer(0)
}
else {
nimps <- as.integer(floor(iter/impute.every))
}
if (mode(save.all.series) != "logical")
stop("Argument save.all.series should be of mode \"logical\".")
length(save.all.series) <- 1L
if (mode(save.worst.series) != "logical")
stop("Argument save.worst.series should be of mode \"logical\".")
length(save.worst.series) <- 1L
if (save.worst.series) {
if (is.null(worst.linear.coef)) {
stop("If save.worst.series, worst.linear.coef must be provided.")
}
}
if (is.null(worst.linear.coef))
worst.linear.coef <- rep(0, ncol(x) * ncol(y) + (ncol(y) *
(ncol(y) + 1L))/2L)
storage.mode(worst.linear.coef) <- "double"
y.imp <- y
if (save.all.series) {
series.length <- as.integer(iter)
series.beta <- array(0, c(ncol(x), ncol(y), iter))
dimnames(series.beta) <- list(rownames(beta), colnames(beta),
NULL)
series.sigma <- array(0, c(ncol(y), ncol(y), iter))
dimnames(series.sigma) <- list(rownames(sigma), colnames(sigma),
NULL)
}
else {
series.length <- as.integer(0)
series.beta <- array(0, c(ncol(x), ncol(y), 0))
series.sigma <- array(0, c(ncol(y), ncol(y), 0))
}
storage.mode(series.beta) <- "double"
storage.mode(series.sigma) <- "double"
imp.list <- array(0, c(nrow(y), ncol(y), nimps))
storage.mode(imp.list) <- "double"
series.worst <- numeric(iter)
storage.mode(series.worst) <- "double"
tmp <- .Fortran("norm_mcmc", n = nrow(y), r = ncol(y), p = ncol(x),
x = x, y = y, mvcode = mvcode, prior.type.int = prior.type.int,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, impute.every = impute.every,
nimps = nimps, save.all.series = save.all.series, save.worst.series = save.worst.series,
worst.linear.coef = worst.linear.coef, series.length = series.length,
beta = beta, sigma = sigma, y.imp = y.imp, iter.actual = integer(1),
series.beta = series.beta, series.sigma = series.sigma,
loglik = numeric(iter), logpost = numeric(iter), series.worst = series.worst,
npatt = integer(1), mis = matrix(logical(1), nrow(y),
ncol(y)), n.in.patt = integer(nrow(y)), n.obs = integer(ncol(y)),
which.patt = integer(nrow(y)), ybar = numeric(ncol(y)),
ysdv = numeric(ncol(y)), imp.list = imp.list, nimps.actual = integer(1),
status = integer(1), msg.len.max = msg.len.max, msg.codes = msg.codes,
msg.len.actual = integer(1), PACKAGE = "norm2")
msg.lines <- msgNorm(tmp$msg.codes, tmp$msg.len.actual)
if (is.null(msg.lines)) {
msg <- "OK"
}
else {
msg <- paste0(msg.lines, collapse = "\n")
}
msg <- paste(msg, "\n", sep = "")
if (msg != "OK\n")
cat(paste("Note: ", msg, sep = ""))
if (tmp$status != 0) {
result <- invisible(NULL)
}
else {
iter <- tmp$iter.actual
if (iter >= 1) {
loglik <- tmp$loglik[1:iter]
logpost <- tmp$logpost[1:iter]
}
else {
loglik <- NULL
logpost <- NULL
}
nimps <- tmp$nimps.actual
if ((impute.every == 0) | (nimps == 0)) {
imp.list <- NULL
impute.every <- NULL
}
else {
imp.list <- as.list(1:nimps)
for (i in 1:nimps) {
imp.list[[i]] <- tmp$imp.list[, , i]
dimnames(imp.list[[i]]) <- list(rownames(y),
colnames(y))
}
}
y <- tmp$y
y[y == tmp$mvcode] <- NA
if (save.all.series & (iter > 0)) {
beta.array <- tmp$series.beta[, , 1:iter, drop = FALSE]
series.beta <- matrix(numeric(1), iter, ncol(x) *
ncol(y))
beta.names <- character(ncol(x) * ncol(y))
posn <- 0
for (j in 1:ncol(x)) {
for (k in 1:ncol(y)) {
posn <- posn + 1
series.beta[, posn] <- beta.array[j, k, ]
beta.names[posn] <- paste(colnames(x)[j], ",",
colnames(y)[k], sep = "")
}
}
colnames(series.beta) <- beta.names
series.beta <- as.ts(series.beta)
sigma.array <- tmp$series.sigma[, , 1:iter, drop = FALSE]
series.sigma <- matrix(numeric(1), iter, ncol(y) *
(ncol(y) + 1)/2)
sigma.names <- character(ncol(y) * (ncol(y) + 1)/2)
posn <- 0
for (k in 1:ncol(y)) {
for (j in k:ncol(y)) {
posn <- posn + 1
series.sigma[, posn] <- sigma.array[j, k, ]
sigma.names[posn] <- paste(colnames(y)[j],
",", colnames(y)[k], sep = "")
}
}
colnames(series.sigma) <- sigma.names
series.sigma <- as.ts(series.sigma)
}
else {
series.beta <- NULL
series.sigma <- NULL
}
if (save.worst.series & (iter > 0)) {
series.worst <- as.ts(tmp$series.worst[1:iter])
}
else {
series.worst <- NULL
worst.linear.coef <- NULL
}
starting.values <- list(beta = beta, sigma = sigma)
if (prior == "ridge") {
prior.sscp <- NULL
}
else if (prior == "uniform") {
prior.df <- NULL
prior.sscp <- NULL
}
else if (prior == "jeffreys") {
prior.df <- NULL
prior.sscp <- NULL
}
miss.patt <- tmp$mis[1:tmp$npatt, , drop = FALSE]
colnames(miss.patt) <- colnames(y)
miss.patt.freq <- tmp$n.in.patt[1:tmp$npatt]
n.obs <- tmp$n.obs
names(n.obs) <- colnames(y)
which.patt <- tmp$which.patt
names(which.patt) <- rownames(y)
ybar <- tmp$ybar
names(ybar) <- colnames(y)
ysdv <- tmp$ysdv
names(ysdv) <- colnames(y)
result <- list(y = y, x = x, method = "MCMC", prior = prior,
prior.df = prior.df, prior.sscp = prior.sscp, iter = iter,
multicycle = multicycle, seeds = seeds, starting.values = starting.values,
param = list(beta = tmp$beta, sigma = tmp$sigma),
loglik = loglik, logpost = logpost, series.worst = series.worst,
series.beta = series.beta, series.sigma = series.sigma,
y.imp = tmp$y.imp, impute.every = impute.every, imp.list = imp.list,
miss.patt = miss.patt, miss.patt.freq = miss.patt.freq,
n.obs = n.obs, which.patt = which.patt, worst.linear.coef = worst.linear.coef,
ybar = ybar, ysdv = ysdv, msg = msg)
class(result) <- "norm"
}
return(result)
}
<bytecode: 0x3a9d910>
<environment: namespace:norm2>
--- function search by body ---
Function mcmcNorm.default in namespace norm2 has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(y) == "data.frame") { : the condition has length > 1
Calls: mcmcNorm -> mcmcNorm.norm -> mcmcNorm.default
Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc
Version: 2.0.2
Check: examples
Result: ERROR
Running examples in ‘norm2-Ex.R’ failed
The error most likely occurred in:
> ### Name: emNorm
> ### Title: EM algorithm for incomplete multivariate normal data
> ### Aliases: emNorm emNorm.default emNorm.formula emNorm.norm
> ### Keywords: multivariate NA
>
> ### ** Examples
>
> ## run EM for marijuana data with strict convergence criterion
> data(marijuana)
> result <- emNorm(marijuana, criterion=1e-06)
*** caught segfault ***
address 0, cause 'memory not mapped'
Traceback:
1: emNorm.default(marijuana, criterion = 1e-06)
2: emNorm(marijuana, criterion = 1e-06)
An irrecoverable exception occurred. R is aborting now ...
Flavor: r-patched-solaris-x86
Version: 2.0.2
Check: tests
Result: ERROR
Running ‘emTests.R’
Running ‘impTests.R’
Comparing ‘impTests.Rout’ to ‘impTests.Rout.save’ ... OK
Running ‘logpostTests.R’
Comparing ‘logpostTests.Rout’ to ‘logpostTests.Rout.save’ ... OK
Running ‘mcmcTests.R’
Comparing ‘mcmcTests.Rout’ to ‘mcmcTests.Rout.save’ ... OK
Running ‘miInferenceTests.R’
Comparing ‘miInferenceTests.Rout’ to ‘miInferenceTests.Rout.save’ ... OK
Running the tests in ‘tests/emTests.R’ failed.
Complete output:
> library(norm2)
>
> ## run EM on fake data with no missing values
> set.seed(1234)
> simdata <- data.frame(
+ Y1=rnorm(6), Y2=rnorm(6), Y3=rnorm(6), X1=rnorm(6) )
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
> print( summary( emResult ) )
Predictor (X) variables:
Mean SD Observed Missing Pct.Missing
(Intercept) 1.0000000 0.000000 6 0 0
X1 0.2068512 1.178512 6 0 0
Response (Y) variables:
Mean SD Observed Missing Pct.Missing
Y1 -0.2092854 1.2953550 6 0 0
Y2 -0.6752401 0.2138685 6 0 0
Y3 -0.2141319 0.6864124 6 0 0
Missingness patterns for response (Y) variables
(. denotes observed value, m denotes missing value)
(variable names are displayed vertically)
(rightmost column is the frequency):
YYY
123
... 6
Method: EM
Prior: "uniform"
Convergence criterion: 1e-05
Iterations: 2
Converged: TRUE
Max. rel. difference: 0
-2 Loglikelihood: -8.3665
-2 Log-posterior density: -8.3665
Worst fraction missing information: 0
Estimated coefficients (beta):
Y1 Y2 Y3
(Intercept) -0.3069981 -0.6785479 -0.2457194
X1 0.4723816 0.0159912 0.1527063
Estimated covariance matrix (sigma):
Y1 Y2 Y3
Y1 1.14001811 0.06789367 0.13397509
Y2 0.06789367 0.03782049 0.03942553
Y3 0.13397509 0.03942553 0.36564511
>
> ## impose missing values and run again
> simdata$Y1[3] <- simdata$Y2[4] <- simdata$Y3[4] <- NA
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
> print( summary( emResult ) )
Predictor (X) variables:
Mean SD Observed Missing Pct.Missing
(Intercept) 1.0000000 0.000000 6 0 0
X1 0.2068512 1.178512 6 0 0
Response (Y) variables:
Mean SD Observed Missing Pct.Missing
Y1 -0.4680307 1.2630567 5 1 16.66667
Y2 -0.6322806 0.2081665 5 1 16.66667
Y3 -0.2349012 0.7653216 5 1 16.66667
Missingness patterns for response (Y) variables
(. denotes observed value, m denotes missing value)
(variable names are displayed vertically)
(rightmost column is the frequency):
YYY
123
... 4
.mm 1
m.. 1
Method: EM
Prior: "uniform"
Convergence criterion: 1e-05
Iterations: 114
Converged: TRUE
Max. rel. difference: 9.8836e-06
-2 Loglikelihood: -10.70465
-2 Log-posterior density: -10.70465
Worst fraction missing information: 0.9313
Estimated coefficients (beta):
Y1 Y2 Y3
(Intercept) -1.0195738 -0.61863204 -0.1874030
X1 0.5166042 -0.01611584 0.1214564
Estimated covariance matrix (sigma):
Y1 Y2 Y3
Y1 1.80844106 -0.06800034 -0.75502100
Y2 -0.06800034 0.03443290 0.05759078
Y3 -0.75502100 0.05759078 0.41668488
>
> ## redundant Y-variable
> simdata$Y3 <- simdata$Y1 + simdata$Y2
> emResult <- emNorm( cbind(Y1,Y2,Y3) ~ X1, data=simdata )
*** caught segfault ***
address 0, cause 'memory not mapped'
Traceback:
1: emNorm.default(y, x = x, intercept = FALSE, iter.max = iter.max, criterion = criterion, estimate.worst = estimate.worst, starting.values = starting.values, prior = prior, prior.df = prior.df, prior.sscp = prior.sscp, ...)
2: emNorm.formula(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
3: emNorm(cbind(Y1, Y2, Y3) ~ X1, data = simdata)
An irrecoverable exception occurred. R is aborting now ...
Flavor: r-patched-solaris-x86