Running ‘BLUP.tests.R’ [1s/2s]
  Running ‘regress.tests.R’ [1s/2s]
  Running ‘regressPaper.R’ [1s/2s]
Running the tests in ‘tests/BLUP.tests.R’ failed.
Complete output:
  > set.seed(1001)
  > library(regress)
  > n <- 101
  > x1 <- runif(n)
  > x2 <- seq(0,1,l=n)
  > z1 <- gl(4,10,n)
  > z2 <- gl(6,1,n)
  > 
  > X <- model.matrix(~1 + x1 + x2)
  > Z1 <- model.matrix(~z1-1)
  > Z2 <- model.matrix(~z2-1)
  > 
  > ## Create the individual random and fixed effects
  > beta <- c(1,2,3)
  > eta1 <- rnorm(ncol(Z1),0,10)
  > eta2 <- rnorm(ncol(Z2),0,10)
  > eps <- rnorm(n,0,3)
  > 
  > ## Combine them into a response
  > y <- X %*% beta + Z1 %*% eta1 + Z2 %*% eta2 + eps
  > 
  > ## Fit the same model again
  > model <- regress(y~1 + x1 + x2,~z1 + z2, identity=TRUE,verbose=2)
  1 sigma = 111.2472 111.2472 111.2472 resid llik = -193.5592 
  1 adjusted sigma = 14.95327 14.95327 14.95327
  2 sigma = 48.05741 34.5205 12.87801 resid llik = -185.4081 
  2 adjusted sigma = 43.27607 31.08598 11.59674 delta.llik = 8.151054 
  3 sigma = 70.09679 47.21641 11.03377 resid llik = -184.5278 
  3 adjusted sigma = 69.55275 46.84995 10.94813 delta.llik = 0.8802895 
  4 sigma = 80.45497 53.4053 10.80752 resid llik = -184.4768 
  4 adjusted sigma = 80.4402 53.3955 10.80554 delta.llik = 0.05097438 
  5 sigma = 81.55823 54.06055 10.79318 resid llik = -184.4764 
  5 adjusted sigma = 81.5581 54.06046 10.79316 delta.llik = 0.0004381634 
   ----------- FAILURE REPORT -------------- 
   --- failure: the condition has length > 1 ---
   --- srcref --- 
  : 
   --- package (from environment) --- 
  regress
   --- call from context --- 
  regress(y ~ 1 + x1 + x2, ~z1 + z2, identity = TRUE, verbose = 2)
   --- 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 ~ 1 + x1 + x2, ~z1 + z2, identity = TRUE, verbose = 2)
  
   --- 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: 0x55bc1c842448>
  <environment: namespace:regress>
   --- function search by body ---
  Function regress in namespace regress has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Fatal error: the condition has length > 1
Running the tests in ‘tests/regress.tests.R’ failed.
Complete output:
  > 
  > ## Test 1, Random Effects Model
  > 
  > library(nlme)
  > library(regress)
  > data(Oats)
  > 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)
  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: 0x557fc0d4fda0>
  <environment: namespace:regress>
   --- function search by body ---
  Function regress in namespace regress has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Fatal error: the condition has length > 1
Running the tests in ‘tests/regressPaper.R’ failed.
Complete output:
  > rm(list=ls())
  > if(FALSE) {
  +     files <- list.files("../R/",full.names=TRUE)
  +     files <- files[-grep("\\~",files)]
  +     for(ff in files) source(ff)
  + } else {
  +     library(regress)
  + }
  > 
  > library(MASS) ## needed for mvrnorm
  > n <- 100
  > mu <- c(1,2)
  > Sigma <- matrix(c(10,5,5,10),2,2)
  > Y <- mvrnorm(n,mu,Sigma)
  > 
  > ## this simulates multivariate normal rvs
  > y <- as.vector(t(Y))
  > X <- kronecker(rep(1,n),diag(1,2))
  > V1 <- matrix(c(1,0,0,0),2,2)
  > V2 <- matrix(c(0,0,0,1),2,2)
  > V3 <- matrix(c(0,1,1,0),2,2)
  > sig1 <- kronecker(diag(1,n),V1)
  > sig2 <- kronecker(diag(1,n),V2)
  > gam <- kronecker(diag(1,n),V3)
  > 
  > reg.obj <- regress(y~X-1,~sig1+sig2+gam,identity=FALSE,start=c(1,1,0.5),verbose=2)
  1 sigma = 1 1 0.5 resid llik = -323.7827 
  1 adjusted sigma = 11.18276 11.18276 5.591382
  2 sigma = 12.86956 11.0899 7.185317 resid llik = -322.0121 
  2 adjusted sigma = 12.82722 11.05341 7.161678 delta.llik = 1.77066 
  3 sigma = 13.03401 11.07697 7.342347 resid llik = -321.9935 
  3 adjusted sigma = 13.03356 11.07659 7.342096 delta.llik = 0.01856986 
  4 sigma = 13.05464 11.07928 7.360389 resid llik = -321.9933 
  4 adjusted sigma = 13.05464 11.07928 7.360386 delta.llik = 0.0001857913 
   ----------- FAILURE REPORT -------------- 
   --- failure: the condition has length > 1 ---
   --- srcref --- 
  : 
   --- package (from environment) --- 
  regress
   --- call from context --- 
  regress(y ~ X - 1, ~sig1 + sig2 + gam, identity = FALSE, start = c(1, 
      1, 0.5), verbose = 2)
   --- 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 ~ X - 1, ~sig1 + sig2 + gam, identity = FALSE, start = c(1, 
      1, 0.5), verbose = 2)
  
   --- 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: 0x562d56038200>
  <environment: namespace:regress>
   --- function search by body ---
  Function regress in namespace regress has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Fatal error: the condition has length > 1