Running examples in 'regress-Ex.R' failed
The error most likely occurred in:

> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: regress
> ### Title: Fit a Gaussian Linear Model with Linear Covariance Structure
> ### Aliases: regress print.regress summary.regress BLUP
> ### Keywords: regression models multivariate
> 
> ### ** Examples
> 
> 
>   ######################
>   ## Comparison with lme
>   ######################
> 
>   ## Example of Random Effects model from Venables and Ripley, page 205
>   library(nlme)
>   library(regress)
> 
>   citation("regress")

To cite package 'regress' in publications use:

  David Clifford, Peter McCullagh (2014). The regress package R package
  version 1.3-15.

  David Clifford, Peter McCullagh (2006), The regress function, R News
  6:2, 6-10

To see these entries in BibTeX format, use 'print(<citation>,
bibtex=TRUE)', 'toBibtex(.)', or set
'options(citation.bibtex.max=999)'.

> 
>   names(Oats) <- c("B","V","N","Y")
>   Oats$N <- as.factor(Oats$N)
> 
>   ## Using regress
>   oats.reg <- regress(Y~N+V,~B+I(B:V),identity=TRUE,verbose=1,data=Oats)
Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
  Consider formula(paste(x, collapse = " ")) instead.
1 sigma = 732.1964 732.1964 732.1964 resid llik = -215.1927 
1 adjusted sigma = 157.9849 157.9849 157.9849
2 sigma = 186.231 133.8389 160.2718 resid llik = -215.0362 
2 adjusted sigma = 185.6373 133.4123 159.761 delta.llik = 0.1565299 
3 sigma = 205.8252 116.8087 161.7195 resid llik = -214.981 
3 adjusted sigma = 205.6644 116.7175 161.5931 delta.llik = 0.05518008 
4 sigma = 213.5958 110.3954 162.4623 resid llik = -214.975 
4 adjusted sigma = 213.5887 110.3917 162.4569 delta.llik = 0.005981576 
5 sigma = 214.3882 109.7628 162.5486 resid llik = -214.9749 
5 adjusted sigma = 214.3882 109.7628 162.5486 delta.llik = 6.213704e-05 
 ----------- FAILURE REPORT -------------- 
 --- failure: the condition has length > 1 ---
 --- srcref --- 
: 
 --- package (from environment) --- 
regress
 --- call from context --- 
regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1, 
    data = Oats)
 --- call from argument --- 
if (error1) {
    cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
    sigma.cov <- matrix(NA, k, k)
}
 --- R stacktrace ---
where 1: regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1, 
    data = Oats)

 --- value of length: 2 type: logical ---
[1] FALSE FALSE
 --- function from context --- 
function (formula, Vformula, identity = TRUE, kernel = NULL, 
    start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL, 
    maxcyc = 50, tol = 1e-04, data) 
{
    if (verbose > 9) 
        cat("Extracting objects from call\n")
    if (missing(data)) 
        data <- environment(formula)
    mf <- model.frame(formula, data = data, na.action = na.pass)
    mf <- eval(mf, parent.frame())
    y <- model.response(mf)
    model <- list()
    model <- c(model, mf)
    if (missing(Vformula)) 
        Vformula <- NULL
    isNA <- apply(is.na(mf), 1, any)
    if (!is.null(Vformula)) {
        V <- model.frame(Vformula, data = data, na.action = na.pass)
        V <- eval(V, parent.frame())
        mfr <- is.na(V)
        if (ncol(mfr) == 1) {
            isNA <- isNA | mfr
        }
        else {
            isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)], 
                1, any)
        }
        rm(mfr)
        Vcoef.names <- names(V)
        V <- as.list(V)
        k <- length(V)
    }
    else {
        V <- NULL
        k <- 0
        Vcoef.names = NULL
    }
    if (ncol(mf) == 1) 
        mf <- cbind(mf, 1)
    X <- model.matrix(formula, mf[!isNA, ])
    y <- y[!isNA]
    n <- length(y)
    Xcolnames <- dimnames(X)[[2]]
    if (is.null(Xcolnames)) {
        Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]), 
            sep = "")
    }
    X <- matrix(X, n, length(X)/n)
    qr <- qr(X)
    rankQ <- n - qr$rank
    if (qr$rank) {
        X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
        Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
    }
    else {
        cat("\nERROR: X has rank 0\n\n")
    }
    if (verbose > 9) 
        cat("Setting up kernel\n")
    if (missing(kernel)) {
        K <- X
        colnames(K) <- Xcolnames
        reml <- TRUE
        kernel <- NULL
    }
    else {
        if (length(kernel) == 1 && kernel > 0) {
            K <- matrix(rep(1, n), n, 1)
            colnames(K) <- c("1")
        }
        if (length(kernel) == 1 && kernel <= 0) {
            K <- Kcolnames <- NULL
            KX <- X
            rankQK <- n
        }
        if (length(kernel) > 1) {
            if (is.matrix(kernel)) {
                K <- kernel[!isNA, ]
            }
            else {
                K <- model.frame(kernel, data = data, na.action = na.pass)
                K <- eval(K, parent.frame())
                if (ncol(K) == 1) {
                  dimNamesK <- dimnames(K)
                  K <- K[!isNA, ]
                  dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
                  K <- data.frame(V1 = K)
                  dimnames(K) <- dimNamesK
                }
                else {
                  K <- K[!isNA, ]
                }
                K <- model.matrix(kernel, K)
            }
        }
        reml <- FALSE
    }
    if (!is.null(K)) {
        Kcolnames <- colnames(K)
        qr <- qr(K)
        rankQK <- n - qr$rank
        if (qr$rank == 0) 
            K <- NULL
        else {
            K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
            Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
            KX <- cbind(K, X)
            qr <- qr(KX)
            KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
        }
    }
    if (missing(maxcyc)) 
        maxcyc <- 50
    if (missing(tol)) 
        tol <- 1e-04
    delta <- 1
    if (verbose > 9) 
        cat("Removing parts of random effects corresponding to missing values\n")
    for (i in 1:k) {
        if (is.matrix(V[[i]])) {
            V[[i]] <- V[[i]][!isNA, !isNA]
        }
        if (is.factor(V[[i]])) {
            V[[i]] <- V[[i]][!isNA]
        }
    }
    In <- diag(rep(1, n), n, n)
    if (identity) {
        V[[k + 1]] <- as.factor(1:n)
        names(V)[k + 1] <- "In"
        k <- k + 1
        Vcoef.names <- c(Vcoef.names, "In")
        Vformula <- as.character(Vformula)
        Vformula[1] <- "~"
        Vformula[2] <- paste(Vformula[2], "+In")
        Vformula <- as.formula(Vformula)
    }
    model <- c(model, V)
    model$formula <- formula
    model$Vformula <- Vformula
    if (!missing(pos)) 
        pos <- as.logical(pos)
    if (missing(pos)) 
        pos <- rep(FALSE, k)
    if (length(pos) < k) 
        cat("Warning: argument pos is only partially specified; additional terms (n=", 
            k - length(pos), ") set to FALSE internally.\n", 
            sep = "")
    pos <- c(pos, rep(FALSE, k))
    pos <- pos[1:k]
    if (verbose > 9) 
        cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
    if (all(sapply(V, is.factor)) & k > 2) {
        SWsolveINDICATOR <- TRUE
    }
    else SWsolveINDICATOR <- FALSE
    Z <- list()
    for (i in 1:length(V)) {
        if (is.factor(V[[i]])) {
            Vi <- model.matrix(~V[[i]] - 1)
            colnames(Vi) <- levels(V[[i]])
            Z[[i]] <- Vi
            V[[i]] <- tcrossprod(Vi)
        }
        else {
            Z[[i]] <- V[[i]]
        }
    }
    names(Z) <- names(V)
    A <- matrix(rep(0, k^2), k, k)
    entries <- expand.grid(1:k, 1:k)
    x <- rep(0, k)
    sigma <- c(1, rep(0, k - 1))
    stats <- rep(0, 0)
    if (missing(taper)) {
        taper <- rep(0.9, maxcyc)
        if (missing(start) && k > 1) 
            taper[1:2] <- c(0.5, 0.7)
    }
    else {
        taper <- pmin(abs(taper), 1)
        if ((l <- length(taper)) < maxcyc) 
            taper <- c(taper, rep(taper[l], maxcyc - l))
    }
    if (!is.null(start)) {
        start <- c(start, rep(1, k))
        start <- start[1:k]
    }
    if (k > 2 && is.null(start)) 
        start <- rep(var(y, na.rm = TRUE), k)
    if (k == 1 && is.null(start)) 
        start <- var(y, na.rm = TRUE)
    if (is.null(start) && k == 2) {
        if (missing(gamVals)) {
            gamVals <- seq(0.01, 0.02, length = 3)^2
            gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3), 
                1 - gamVals))
            gamVals <- 0.5
        }
        if (length(gamVals) > 1) {
            if (verbose >= 1) 
                cat("Evaluating the llik at gamma = \n")
            if (verbose >= 1) 
                cat(gamVals)
            if (verbose >= 1) 
                cat("\n")
            reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
            llik <- reg.obj$llik
            llik <- as.double(llik)
            if (verbose >= 2) 
                cat(llik, "\n")
            gam <- gamVals[llik == max(llik)]
            gam <- gam[1]
            if (verbose >= 2) 
                cat("MLE is near", gam, "and llik =", max(llik), 
                  "there\n")
        }
        if (length(gamVals) == 1) {
            gam <- gamVals[1]
            reg.obj <- list(rms = var(y))
        }
        start <- c(1 - gam, gam) * reg.obj$rms
        if (gam == 0.9999) {
            taper[1] <- taper[1]/100
            maxcyc <- maxcyc * 10
        }
        if (verbose >= 1) 
            cat(c("start algorithm at", round(start, 4), "\n"))
    }
    if (is.null(start) & k > 2) {
        LLvals <- NULL
        V2 <- V[[2]]
        for (ii in 3:k) V2 <- V2 + V[[ii]]
        LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
        V2 <- V[[1]] + V2
        for (ii in 1:k) {
            V2 <- V2 - V[[ii]]
            LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
        }
        best <- which.max(LLvals)
        if (verbose) {
            cat("Checking starting points\n")
            cat("llik values of", LLvals, "\n")
        }
        if (best == 1) {
            start <- rep(var(y, na.rm = TRUE), k)
        }
        else {
            start <- rep(0.25, k)
            start[best] <- 0.75
        }
    }
    sigma <- coef <- start
    coef[pos] <- log(sigma[pos])
    coef[!pos] <- sigma[!pos]
    T <- vector("list", length = k)
    for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
    for (cycle in 1:maxcyc) {
        ind <- which(pos)
        if (length(ind)) {
            coef[ind] <- pmin(coef[ind], 20)
            coef[ind] <- pmax(coef[ind], -20)
            sigma[ind] <- exp(coef[ind])
        }
        if (verbose) {
            cat(cycle, "sigma =", sigma)
        }
        if (!SWsolveINDICATOR) {
            Sigma <- 0
            for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
            W <- solve(Sigma, In)
        }
        else {
            W <- SWsolve2(Z[1:(k - 1)], sigma)
        }
        if (is.null(K)) 
            WQK <- W
        else {
            WK <- W %*% K
            WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
        }
        if (reml) 
            WQX <- WQK
        else {
            WX <- W %*% KX
            WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
        }
        rss <- as.numeric(t(y) %*% WQX %*% y)
        sigma <- sigma * rss/rankQK
        coef[!pos] <- sigma[!pos]
        coef[pos] <- log(sigma[pos])
        WQK <- WQK * rankQK/rss
        WQX <- WQX * rankQK/rss
        rss <- rankQK
        eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values, 
            decreasing = TRUE)[1:rankQK]
        if (any(eig < 0)) {
            cat("error: Sigma is not positive definite on contrasts: range(eig)=", 
                range(eig), "\n")
            WQK <- WQK + (tol - min(eig)) * diag(nobs)
            eig <- eig + tol - min(eig)
        }
        ldet <- sum(log(eig))
        llik <- ldet/2 - rss/2
        if (cycle == 1) 
            llik0 <- llik
        delta.llik <- llik - llik0
        llik0 <- llik
        if (verbose) {
            if (reml) 
                cat(" resid llik =", llik, "\n")
            else cat(" llik =", llik, "\n")
            cat(cycle, "adjusted sigma =", sigma)
            if (cycle > 1) {
                if (reml) 
                  cat(" delta.llik =", delta.llik, "\n")
                else cat(" delta.llik =", delta.llik, "\n")
            }
            else cat("\n")
        }
        x <- NULL
        var.components <- rep(1, k)
        ind <- which(pos)
        if (length(ind)) 
            var.components[ind] <- sigma[ind]
        if (!SWsolveINDICATOR) {
            if (identity) {
                T[[k]] <- WQK
                if (k > 1) {
                  for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
                }
            }
            else {
                for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
            }
        }
        else {
            if (identity) {
                T[[k]] <- WQK
                if (k > 1) {
                  for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*% 
                    Z[[ii]], Z[[ii]])
                }
            }
            else {
                for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*% 
                  Z[[ii]], Z[[ii]])
            }
        }
        x <- sapply(T, function(x) as.numeric(t(y) %*% x %*% 
            WQX %*% y - sum(diag(x))))
        x <- x * var.components
        ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] * 
            var.components[x[2]]
        aa <- apply(entries, 1, ff)
        A[as.matrix(entries)] <- aa
        stats <- c(stats, llik, sigma[1:k], x[1:k])
        if (verbose >= 9) {
        }
        A.svd <- ginv(A)
        x <- A.svd %*% x
        if (qr(A)$rank < k) {
            if (cycle == 1) {
                if (verbose) {
                  cat("Warning: Non identifiable dispersion model\n")
                  cat(sigma)
                  cat("\n")
                }
            }
        }
        coef <- coef + taper[cycle] * x
        sigma[!pos] <- coef[!pos]
        sigma[pos] <- exp(coef[pos])
        if (cycle > 1 & abs(delta.llik) < tol * 10) 
            break
        if (max(abs(x)) < tol) 
            break
    }
    if (!SWsolveINDICATOR) {
        Sigma <- 0
        for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
        W <- solve(Sigma, In)
    }
    else {
        W <- SWsolve2(Z[1:(k - 1)], sigma)
    }
    if (cycle == maxcyc) {
        if (verbose) 
            cat("WARNING:  maximum number of cycles reached before convergence\n")
    }
    stats <- as.numeric(stats)
    stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
    colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""), 
        paste("der_", Vcoef.names, sep = ""))
    WX <- W %*% X
    XtWX <- crossprod(X, WX)
    cov <- XtWX
    cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
    beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]], 
        dim(X)[2], dim(X)[2])
    cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
    beta <- cov %*% y
    beta <- matrix(beta, length(beta), 1)
    row.names(beta) <- Xcolnames
    beta.se <- sqrt(abs(diag(beta.cov)))
    pos.cov <- (diag(beta.cov) < 0)
    beta.se[pos.cov] <- NA
    beta.se <- matrix(beta.se, length(beta.se), 1)
    row.names(beta.se) <- Xcolnames
    rms <- rss/rankQ
    fitted.values <- X %*% beta
    Q <- In - X %*% cov
    predicted <- NULL
    predictedVariance <- NULL
    predictedVariance2 <- NULL
    if (identity) {
        gam <- sigma[k]
        if (SWsolveINDICATOR) {
            Sigma <- 0
            for (i in 1:k) {
                Sigma <- Sigma + V[[i]] * sigma[i]
            }
        }
        predicted <- fitted.values + (Sigma - gam * In) %*% W %*% 
            (y - fitted.values)
        predictedVariance <- 2 * gam - gam^2 * diag(W)
        predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*% 
            t(WX))
    }
    sigma.cov <- (A.svd[1:k, 1:k] * 2)
    FI <- A/2
    FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
    FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
    names(sigma) <- Vcoef.names
    sigma.cov <- try(ginv(FI.c), silent = TRUE)
    error1 <- (class(sigma.cov) == "try-error")
    if (error1) {
        cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
        sigma.cov <- matrix(NA, k, k)
    }
    rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
    if (!error1) {
        if (any(sigma[pos]^2 < 1e-04)) 
            cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
    }
    result <- list(trace = stats, llik = llik, cycle = cycle, 
        rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se, 
        sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k], 
        W = W, Q = Q, fitted = fitted.values, predicted = predicted, 
        predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2, 
        pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula, 
        Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
    class(result) <- "regress"
    return(result)
}
<bytecode: 0x2fd0480>
<environment: namespace:regress>
 --- function search by body ---
Function regress in namespace regress has this body.
 ----------- END OF FAILURE REPORT -------------- 
Error in if (error1) { : the condition has length > 1
Calls: regress
Execution halted