Running 'bdstest.R' [3s/3s]
  Comparing 'bdstest.Rout' to 'bdstest.Rout.save' ... OK
  Running 'brfit.R' [7s/8s]
  Running 'lmekin1.R' [5s/6s]
  Running 'lmekin2.R' [3s/4s]
  Comparing 'lmekin2.Rout' to 'lmekin2.Rout.save' ...50,52d49
< Warning message:
< In if (ismat(kmat)) { :
<   the condition has length > 1 and only the first element will be used
  Running 'nest1.R' [3s/3s]
  Comparing 'nest1.Rout' to 'nest1.Rout.save' ... OK
  Running 'nest2.R' [3s/3s]
  Comparing 'nest2.Rout' to 'nest2.Rout.save' ... OK
  Running 'refine.R' [4s/4s]
  Comparing 'refine.Rout' to 'refine.Rout.save' ... OK
  Running 'refine1.R' [4s/5s]
  Comparing 'refine1.Rout' to 'refine1.Rout.save' ... OK
  Running 'refine2.R' [4s/5s]
  Comparing 'refine2.Rout' to 'refine2.Rout.save' ... OK
  Running 'refine3.R' [4s/4s]
  Comparing 'refine3.Rout' to 'refine3.Rout.save' ... OK
  Running 'shrink.R' [4s/4s]
  Running 'slope1.R' [17s/19s]
  Comparing 'slope1.Rout' to 'slope1.Rout.save' ... OK
  Running 'slope2.R' [4s/4s]
  Running 'slope3.R' [5s/6s]
  Comparing 'slope3.Rout' to 'slope3.Rout.save' ...85,87d84
< Warning message:
< In if (ismat(kmat)) { :
<   the condition has length > 1 and only the first element will be used
137,139d133
< Warning message:
< In if (ismat(kmat)) { :
<   the condition has length > 1 and only the first element will be used
  Running 'slope4.R' [12s/13s]
  Comparing 'slope4.Rout' to 'slope4.Rout.save' ...181,183d180
< Warning message:
< In if (ismat(kmat)) { :
<   the condition has length > 1 and only the first element will be used
198,200d194
< Warning message:
< In if (ismat(kmat)) { :
<   the condition has length > 1 and only the first element will be used
  Running 'strata1.R' [3s/3s]
  Comparing 'strata1.Rout' to 'strata1.Rout.save' ... OK
  Running 'test0.R' [4s/5s]
  Running 'test1.R' [3s/4s]
  Comparing 'test1.Rout' to 'test1.Rout.save' ... OK
  Running 'test2.R' [3s/4s]
  Comparing 'test2.Rout' to 'test2.Rout.save' ... OK
  Running 'test3.R' [3s/3s]
  Comparing 'test3.Rout' to 'test3.Rout.save' ... OK
  Running 'test4.R' [3s/4s]
  Comparing 'test4.Rout' to 'test4.Rout.save' ... OK
  Running 'twoterm.R' [3s/4s]
  Comparing 'twoterm.Rout' to 'twoterm.Rout.save' ... OK
Running the tests in 'tests/lmekin1.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > require(nlme)
  Loading required package: nlme
  > aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > tdata <- eortc
  > tdata$center2 <- factor(tdata$center)
  > tdata$trt2 <- factor(tdata$trt)
  > 
  > fit1 <- lme(y ~ trt, random= ~ 1|center2, data=tdata,
  +             method="ML")
  > fit2 <- lmekin(y ~ trt + (1|center), tdata)
  > 
  > #I get the same loglik but slightly different coefficients.
  > # (The loglik is flat on top, and we are using different starting
  > #  estimates.)  So we need to use a toloerance for the test.
  > vcomp <- function(m1, m2, ...) {
  +     temp1 <- as.numeric(VarCorr(m1)[,2])  #lme gives a matrix of text
  + #    temp2 <- m2$sigma^2 * c(unlist(VarCorr(m2)), 1) #lmekin the variance struct
  +     temp2 <- c(unlist(VarCorr(m2)), m2$sigma^2) #lmekin the variance struct
  +     aeq(temp1, sqrt(temp2), ...)
  + }
  > 
  > aeq(fit1$logLik, fit2$loglik)
  [1] TRUE
  > aeq(fit1$sigma, fit2$sigma, tol=1e-4)
  [1] TRUE
  > aeq(fixef(fit1), fixef(fit2), tol=1e-4)
  [1] TRUE
  > vcomp(fit1, fit2, tol=1e-3)
  [1] TRUE
  > 
  > # A second fit, using the matrix form
  > ucen <- sort(unique(tdata$center))
  > kmat <- diag(length(ucen))
  > dimnames(kmat) <- list(ucen, ucen)
  > fit3 <- lmekin(y ~ trt + (1|center), tdata, varlist=kmat)
   ----------- FAILURE REPORT -------------- 
   --- failure: the condition has length > 1 ---
   --- srcref --- 
  : 
   --- package (from environment) --- 
  coxme
   --- call from context --- 
  lmekin(y ~ trt + (1 | center), tdata, varlist = kmat)
   --- call from argument --- 
  if (class(varlist) == "coxmevar") varlist <- list(varlist) else if (ismat(varlist)) varlist <- list(coxmeMlist(list(varlist))) else {
      if (!is.list(varlist)) 
          stop("Invalid varlist argument")
      if (all(sapply(varlist, ismat))) {
          if (nrandom > 1) 
              stop(paste("An unlabeled list of matrices is", "ambiguous when there are multiple random terms"))
          else varlist <- list(coxmeMlist(varlist))
      }
      else {
          for (i in 1:length(varlist)) {
              if (is.function(varlist[[i]])) 
                  varlist[[i]] <- varlist[[i]]()
              if (ismat(varlist[[i]])) 
                  varlist[[i]] <- coxmeMlist(list(varlist[[i]]))
              if (class(varlist[[i]]) != "coxmevar") {
                  if (is.list(varlist[[i]])) {
                    if (all(sapply(varlist[[i]], ismat))) 
                      varlist[[i]] <- coxmeMlist(varlist[[i]])
                    else stop("Invalid varlist element")
                  }
                  else stop("Invalid varlist element")
              }
          }
      }
  }
   --- R stacktrace ---
  where 1: lmekin(y ~ trt + (1 | center), tdata, varlist = kmat)
  
   --- value of length: 2 type: logical ---
  [1] FALSE FALSE
   --- function from context --- 
  function (formula, data, weights, subset, na.action, control, 
      varlist, vfixed, vinit, method = c("ML", "REML"), x = FALSE, 
      y = FALSE, model = FALSE, random, fixed, variance, ...) 
  {
      Call <- match.call()
      sparse <- c(1, 0)
      if (!missing(fixed)) {
          if (missing(formula)) {
              formula <- fixed
              warning("The 'fixed' argument of lmekin is depreciated")
          }
          else stop("Both a fixed and a formula argument are present")
      }
      if (!missing(random)) {
          warning("The random argument of lmekin is depreciated")
          if (class(random) != "formula" || length(random) != 2) 
              stop("Invalid random formula")
          j <- length(formula)
          formula[[j]] <- call("+", formula[[j]], call("(", random[[2]]))
      }
      if (!missing(variance)) {
          warning("The variance argument of lmekin is depreciated")
          vfixed <- variance
      }
      method <- match.arg(method)
      temp <- call("model.frame", formula = subbar(formula))
      for (i in c("data", "subset", "weights", "na.action")) if (!is.null(Call[[i]])) 
          temp[[i]] <- Call[[i]]
      m <- eval.parent(temp)
      Y <- model.extract(m, "response")
      n <- length(Y)
      if (n == 0) 
          stop("data has no observations")
      weights <- model.weights(m)
      if (length(weights) == 0) 
          weights <- rep(1, n)
      else if (any(weights <= 0)) 
          stop("Negative or zero weights are not allowed")
      offset <- model.offset(m)
      if (length(offset) == 0) 
          offset <- rep(0, n)
      pterms <- sapply(m, inherits, "coxph.penalty")
      if (any(pterms)) {
          stop("You cannot have penalized terms in lmekin")
      }
      if (missing(control)) 
          control <- lmekin.control(...)
      flist <- formula1(formula)
      if (hasAbar(flist$fixed)) 
          stop("Invalid formula: a '|' outside of a valid random effects term")
      special <- c("strata", "cluster")
      Terms <- terms(flist$fixed, special)
      if (length(attr(Terms, "specials")$strata)) 
          stop("A strata term is invalid in lmekin")
      if (length(attr(Terms, "specials")$cluster)) 
          stop("A cluster term is invalid in lmekin")
      X <- model.matrix(Terms, m)
      nrandom <- length(flist$random)
      if (nrandom == 0) 
          stop("No random effects terms found")
      vparm <- vector("list", nrandom)
      is.variance <- rep(TRUE, nrandom)
      ismat <- function(x) {
          inherits(x, "matrix") || inherits(x, "bdsmatrix") | inherits(x, 
              "Matrix")
      }
      if (missing(varlist) || is.null(varlist)) {
          varlist <- vector("list", nrandom)
          for (i in 1:nrandom) varlist[[i]] <- coxmeFull()
      }
      else {
          if (is.function(varlist)) 
              varlist <- varlist()
          if (class(varlist) == "coxmevar") 
              varlist <- list(varlist)
          else if (ismat(varlist)) 
              varlist <- list(coxmeMlist(list(varlist)))
          else {
              if (!is.list(varlist)) 
                  stop("Invalid varlist argument")
              if (all(sapply(varlist, ismat))) {
                  if (nrandom > 1) 
                    stop(paste("An unlabeled list of matrices is", 
                      "ambiguous when there are multiple random terms"))
                  else varlist <- list(coxmeMlist(varlist))
              }
              else {
                  for (i in 1:length(varlist)) {
                    if (is.function(varlist[[i]])) 
                      varlist[[i]] <- varlist[[i]]()
                    if (ismat(varlist[[i]])) 
                      varlist[[i]] <- coxmeMlist(list(varlist[[i]]))
                    if (class(varlist[[i]]) != "coxmevar") {
                      if (is.list(varlist[[i]])) {
                        if (all(sapply(varlist[[i]], ismat))) 
                          varlist[[i]] <- coxmeMlist(varlist[[i]])
                        else stop("Invalid varlist element")
                      }
                      else stop("Invalid varlist element")
                    }
                  }
              }
          }
          while (length(varlist) < nrandom) varlist <- c(varlist, 
              list(coxmeFull()))
      }
      if (!is.null(names(varlist))) {
          vname <- names(varlist)
          stop("Cannot (yet) have a names varlist")
          indx <- pmatch(vname, names(random), nomatch = 0)
          if (any(indx == 0 & vname != "")) 
              stop(paste("Varlist element not matched:", vname[indx == 
                  0 & vname != ""]))
          if (any(indx > 0)) {
              temp <- vector("list", nrandom)
              temp[indx] <- varlist[indx > 0]
              temp[-indx] <- varlist[indx == 0]
              varlist <- temp
          }
      }
      check <- sapply(varlist, function(x) {
          fname <- c("initialize", "generate", "wrapup")
          indx <- match(fname, names(x))
          if (any(is.na(x))) 
              stop(paste("Member not found in variance function:", 
                  fname(is.na(indx))))
          if (length(x) != 3 || any(!sapply(x, is.function))) 
              stop("Varlist objects must consist of exaclty three functions")
      })
      getcmat <- function(x, mf) {
          if (is.null(x) || x == 1) 
              return(NULL)
          Terms <- terms(eval(call("~", x)))
          attr(Terms, "intercept") <- 0
          varnames <- attr(Terms, "term.labels")
          ftemp <- sapply(mf[varnames], is.factor)
          if (any(ftemp)) {
              clist <- lapply(mf[varnames[ftemp]], function(x) diag(length(levels(x))))
              model.matrix(Terms, mf, contrasts.arg = clist)
          }
          else model.matrix(Terms, mf)
      }
      getGroupNames <- function(x) {
          if (is.call(x) && x[[1]] == as.name("/")) 
              c(getGroupNames(x[[2]]), getGroupNames(x[[3]]))
          else deparse(x)
      }
      getgroups <- function(x, mf) {
          if (is.numeric(x)) 
              return(NULL)
          varname <- getGroupNames(x)
          indx <- match(varname, names(mf), nomatch = 0)
          if (any(indx == 0)) 
              stop(paste("Invalid grouping factor", varname[indx == 
                  0]))
          else data.frame(lapply(mf[indx], as.factor))
      }
      if (missing(vinit) || is.null(vinit)) 
          vinit <- vector("list", nrandom)
      else {
          if (nrandom == 1) {
              if (is.numeric(vinit)) 
                  vinit <- list(vinit)
              else if (is.list(vinit)) 
                  vinit <- list(unlist(vinit))
          }
          if (!is.list(vinit)) 
              stop("Invalid value for `vinit` parameter")
          if (length(vinit) > nrandom) 
              stop(paste("Vinit must be a list of length", nrandom))
          if (!all(sapply(vinit, function(x) (is.null(x) || is.numeric(x))))) 
              stop("Vinit must contain numeric values")
          if (length(vinit) < nrandom) 
              vinit <- c(vinit, vector("list", nrandom - length(vinit)))
          tname <- names(vinit)
          if (!is.null(tname)) {
              stop("Named initial values not yet supported")
          }
      }
      if (missing(vfixed) || is.null(vfixed)) 
          vfixed <- vector("list", nrandom)
      else {
          if (nrandom == 1) {
              if (is.numeric(vfixed)) 
                  vfixed <- list(vfixed)
              else if (is.list(vfixed)) 
                  vfixed <- list(unlist(vfixed))
          }
          if (!is.list(vfixed)) 
              stop("Invalid value for `vfixed` parameter")
          if (length(vfixed) > nrandom) 
              stop(paste("Vfixed must be a list of length", nrandom))
          if (!all(sapply(vfixed, function(x) (is.null(x) || is.numeric(x))))) 
              stop("Vfixed must contain numeric values")
          if (length(vfixed) < nrandom) 
              vfixed <- c(vfixed, vector("list", nrandom - length(vfixed)))
          tname <- names(vfixed)
          if (!is.null(tname)) {
              temp <- pmatch(tname, names(flist$random), nomatch = 0)
              temp <- c(temp, (1:nrandom)[-temp])
              vfixed <- vfixed[temp]
          }
      }
      newzmat <- function(xmat, xmap) {
          n <- nrow(xmap)
          newz <- matrix(0, nrow = n, ncol = max(xmap))
          for (i in 1:ncol(xmap)) newz[cbind(1:n, xmap[, i])] <- xmat[, 
              i]
          newz
      }
      fmat <- zmat <- matrix(0, nrow = n, ncol = 0)
      ntheta <- integer(nrandom)
      ncoef <- matrix(0L, nrandom, 2, dimnames = list(NULL, c("intercept", 
          "slope")))
      itheta <- NULL
      for (i in 1:nrandom) {
          f2 <- formula2(flist$random[[i]])
          if (f2$intercept & f2$group == 1) 
              stop(paste("Error in random term ", i, ": Random intercepts require a grouping variable", 
                  sep = ""))
          vfun <- varlist[[i]]
          if (!is.null(f2$interaction)) 
              stop("Interactions not yet written")
          cmat <- getcmat(f2$fixed, m)
          groups <- getgroups(f2$group, m)
          ifun <- vfun$initialize(vinit[[i]], vfixed[[i]], intercept = f2$intercept, 
              groups, cmat, control)
          if (!is.null(ifun$error)) 
              stop(paste("In random term ", i, ": ", ifun$error, 
                  sep = ""))
          vparm[[i]] <- ifun$parms
          if (!is.null(ifun$is.variance)) 
              is.variance[i] <- ifun$is.variance
          itheta <- c(itheta, ifun$theta)
          ntheta[i] <- length(ifun$theta)
          if (f2$intercept) {
              if (!is.matrix(ifun$imap) || nrow(ifun$imap) != n) 
                  stop(paste("In random term ", i, ": Invalid intercept matrix F", 
                    sep = ""))
              temp <- sort(unique(c(ifun$imap)))
              if (any(temp != 1:length(temp))) 
                  stop(paste("In random term ", i, ": intercept matrix has an invalid element", 
                    sep = ""))
              if (ncol(fmat) > 0) 
                  fmat <- cbind(fmat, ifun$imap + max(fmat))
              else fmat <- ifun$imap
              ncoef[i, 1] <- 1 + max(ifun$imap) - min(ifun$imap)
          }
          if (length(cmat) > 0) {
              if (is.null(ifun$xmap) || is.null(ifun$X) || !is.matrix(ifun$xmap) || 
                  !is.matrix(ifun$X) || nrow(ifun$xmap) != n || 
                  nrow(ifun$X) != n || ncol(ifun$xmap) != ncol(ifun$X)) 
                  stop(paste("In random term ", i, "invalid X/xmap pair"))
              if (f2$intercept) 
                  xmap <- ifun$xmap - max(ifun$imap)
              else xmap <- ifun$xmap
              if (any(sort(unique(c(xmap))) != 1:max(xmap))) 
                  stop(paste("In random term ", i, ": xmap matrix has an invalid element", 
                    sep = ""))
              temp <- newzmat(ifun$X, xmap)
              ncoef[i, 2] <- ncol(temp)
              zmat <- cbind(zmat, temp)
          }
      }
      if (any(is.variance) & !all(is.variance)) 
          stop("All variance terms must have the same is.variance setting")
      kfun <- function(theta, varlist, vparm, ntheta, ncoef) {
          nrandom <- length(varlist)
          sindex <- rep(1:nrandom, ntheta)
          tmat <- varlist[[1]]$generate(theta[sindex == 1], vparm[[1]])
          dd <- dim(tmat)
          if (length(dd) != 2 || any(dd != rep(ncoef[1, 1] + ncoef[1, 
              2], 2))) 
              stop("Incorrect dimensions for generated penalty matrix, term 1")
          if (!inherits(tmat, "bdsmatrix")) 
              tmat <- bdsmatrix(blocksize = integer(0), blocks = numeric(0), 
                  rmat = tmat)
          if (nrandom == 1) 
              return(tmat)
          nsparse <- sum(tmat@blocksize)
          nrow.R <- sum(ncoef)
          ncol.R <- nrow.R - nsparse
          R <- matrix(0, nrow.R, ncol.R)
          indx1 <- 0
          indx2 <- sum(ncoef[, 1]) - nsparse
          if (ncol(tmat) > nsparse) {
              if (ncoef[1, 1] > nsparse) {
                  irow <- 1:ncoef[1, 1]
                  j <- ncoef[1, 1] - nsparse
                  R[irow, 1:j] <- tmat@rmat[irow, 1:j]
                  indx1 <- j
                  if (ncoef[1, 2] > 0) {
                    k <- 1:ncoef[1, 2]
                    R[irow, k + indx2 - nsparse] <- tmat@rmat[irow, 
                      k + j]
                  }
              }
              else j <- 0
              if (ncoef[1, 2] > 0) {
                  k <- 1:ncoef[1, 2]
                  R[k + indx2 + nsparse, k + indx2] <- tmat@rmat[k + 
                    indx1, j + k]
                  indx2 <- indx2 + ncoef[1, 2]
              }
          }
          for (i in 2:nrandom) {
              temp <- as.matrix(varlist[[i]]$generate(theta[sindex == 
                  i], vparm[[i]]))
              if (any(dim(temp) != rep(ncoef[i, 1] + ncoef[i, 2], 
                  2))) 
                  stop(paste("Invalid dimension for generated penalty matrix, term", 
                    i))
              if (ncoef[i, 1] > 0) {
                  j <- ncoef[i, 1]
                  R[indx1 + 1:j + nsparse, indx1 + 1:j] <- temp[1:j, 
                    1:j]
                  if (ncoef[i, 2] > 0) {
                    k <- 1:ncoef[i, 2]
                    R[indx1 + 1:j + nsparse, indx2 + k] <- temp[1:j, 
                      k + j]
                    R[indx2 + k + nsparse, indx2 + k] <- temp[k + 
                      j, k + j]
                  }
              }
              else if (ncoef[i, 2] > 0) {
                  k <- 1:ncoef[i, 2]
                  R[indx2 + k + nsparse, indx2 + k] <- temp
              }
              indx1 <- indx1 + ncoef[i, 1]
              indx2 <- indx2 + ncoef[i, 2]
          }
          bdsmatrix(blocksize = tmat@blocksize, blocks = tmat@blocks, 
              rmat = R)
      }
      itemp <- split(row(fmat), fmat)
      zstar1 <- new("dgCMatrix", i = as.integer(unlist(itemp) - 
          1), p = as.integer(c(0, cumsum(unlist(lapply(itemp, length))))), 
          Dim = as.integer(c(n, max(fmat))), Dimnames = list(NULL, 
              NULL), x = rep(1, length(fmat)), factors = list())
      if (length(zmat) > 0) {
          zstar1 <- cbind(zstar1, as(Matrix(zmat), "dgCMatrix"))
      }
      nfrail <- ncol(zstar1)
      nvar <- ncol(X)
      if (nvar == 0) 
          xstar <- NULL
      else xstar <- rbind(Matrix(X, sparse = FALSE), matrix(0, 
          nrow = nfrail, ncol = ncol(X)))
      ystar <- c(Y, rep(0, nfrail))
      mydiag <- function(x) {
          if (class(x) == "sparseQR") 
              diag(x@R)
          else diag(qr.R(x))
      }
      logfun <- function(theta, best = 0) {
          vmat <- kfun(theta, varlist, vparm, ntheta, ncoef)
          Delta <- t(solve(chol(as(vmat, "dsCMatrix"), pivot = FALSE)))
          zstar <- rbind(zstar1, Delta)
          qr1 <- qr(zstar)
          dd <- mydiag(qr1)
          cvec <- as.vector(qr.qty(qr1, ystar))[-(1:nfrail)]
          if (nvar > 0) {
              qr2 <- qr(qr.qty(qr1, xstar)[-(1:nfrail), ])
              cvec <- qr.qty(qr2, cvec)[-(1:nvar)]
              if (method != "ML") 
                  dd <- c(dd, mydiag(qr2))
          }
          loglik <- sum(log(abs(diag(Delta)))) - sum(log(abs(dd)))
          if (method == "ML") 
              loglik <- loglik - 0.5 * n * log(sum(cvec^2))
          else loglik <- loglik - 0.5 * length(cvec) * log(sum(cvec^2))
          best - (loglik + 1)
      }
      nstart <- sapply(itheta, length)
      if (length(nstart) == 0) 
          theta <- NULL
      else {
          testvals <- do.call(expand.grid, itheta)
          bestlog <- NULL
          for (i in 1:nrow(testvals)) {
              ll <- logfun(as.numeric(testvals[i, ]))
              if (is.finite(ll)) {
                  if (is.null(bestlog) || ll < bestlog) {
                    bestlog <- ll
                    theta <- as.numeric(testvals[i, ])
                  }
              }
          }
          if (is.null(bestlog)) 
              stop("No starting estimate was successful")
          optpar <- control$optpar
          optpar$hessian <- TRUE
          mfit <- do.call("optim", c(list(par = theta, fn = logfun, 
              gr = NULL, best = bestlog), optpar))
          theta <- mfit$par
      }
      vmat <- kfun(theta, varlist, vparm, ntheta, ncoef)
      Delta <- t(solve(chol(as(vmat, "dsCMatrix"), pivot = FALSE)))
      zstar <- rbind(zstar1, Delta)
      qr1 <- qr(zstar)
      dd <- mydiag(qr1)
      ctemp <- as.vector(qr.qty(qr1, ystar))
      cvec <- ctemp[-(1:nfrail)]
      if (is.null(xstar)) {
          rcoef <- qr.coef(qr1, ystar)
          yhat <- qr.fitted(qr1, ystar)
      }
      else {
          qtx <- qr.qty(qr1, xstar)
          qr2 <- qr(qtx[-(1:nfrail), , drop = F])
          if (method != "ML") 
              dd <- c(dd, mydiag(qr2))
          fcoef <- qr.coef(qr2, cvec)
          yresid <- ystar - xstar %*% fcoef
          rcoef <- qr.coef(qr1, yresid)
          cvec <- qr.qty(qr2, cvec)[-(1:nvar)]
          if (class(qr2) == "sparseQR") 
              varmat <- chol2inv(qr2@R)
          else varmat <- chol2inv(qr.R(qr2))
          yhat <- as.vector(zstar1 %*% rcoef + X %*% fcoef)
      }
      if (method == "ML") {
          sigma2 <- sum(cvec^2)/n
          loglik <- sum(log(abs(diag(Delta)))) - (sum(log(abs(dd))) + 
              0.5 * n * (log(2 * pi) + 1 + log(sigma2)))
      }
      else {
          np <- length(cvec)
          sigma2 <- mean(cvec^2)
          loglik <- sum(log(abs(diag(Delta)))) - (sum(log(abs(dd))) + 
              0.5 * np * (log(2 * pi) + 1 + log(sigma2)))
      }
      if (FALSE) {
          zx <- cbind(zstar, as(xstar, class(zstar)))
          qr3 <- qr(zx)
          cvec3 <- qr.qty(qr3, ystar)[-(1:(nvar + nfrail))]
          if (method == "ML") 
              dd3 <- (diag(myqrr(qr3)))[1:nfrail]
          else dd3 <- (diag(myqrr(qr3)))[1:(nfrail + nvar)]
          acoef <- qr.coef(qr3, ystar)
          browser()
      }
      newtheta <- random.coef <- list()
      nrandom <- length(varlist)
      sindex <- rep(1:nrandom, ntheta)
      bindex <- rep(1:nrandom, rowSums(ncoef))
      for (i in 1:nrandom) {
          temp <- varlist[[i]]$wrapup(theta[sindex == i], rcoef[bindex == 
              i], vparm[[i]])
          newtheta <- c(newtheta, lapply(temp$theta, function(x) x * 
              sigma2))
          if (!is.list(temp$b)) {
              temp$b <- list(temp$b)
              names(temp$b) <- paste("Random", i, sep = "")
          }
          random.coef <- c(random.coef, temp$b)
      }
      if (length(fcoef) > 0) {
          nvar <- length(fcoef)
          fit <- list(coefficients = list(fixed = fcoef, random = random.coef), 
              var = varmat * sigma2, vcoef = newtheta, residuals = Y - 
                  yhat, method = method, loglik = loglik, sigma = sqrt(sigma2), 
              n = n, call = Call)
      }
      else fit <- list(coefficients = list(fixed = NULL, random = random.coef), 
          vcoef = newtheta, residuals = Y - yhat, method = method, 
          loglik = loglik, sigma = sqrt(sigma2), n = n, call = Call)
      if (!is.null(theta)) {
          fit$rvar <- mfit$hessian
          fit$iter <- mfit$counts
      }
      if (x) 
          fit$x <- X
      if (y) 
          fit$y <- Y
      if (model) 
          fit$model <- m
      na.action <- attr(m, "na.action")
      if (length(na.action)) 
          fit$na.action <- na.action
      class(fit) <- "lmekin"
      fit
  }
  <environment: namespace:coxme>
   --- function search by body ---
  Function lmekin in namespace coxme has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Error in if (class(varlist) == "coxmevar") varlist <- list(varlist) else if (ismat(varlist)) varlist <- list(coxmeMlist(list(varlist))) else { : 
    the condition has length > 1
  Calls: lmekin
  Execution halted
Running the tests in 'tests/shrink.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > #
  > # Variable shrinkage
  > #
  > ecog0 <- 1*(lung$ph.ecog==0)
  > ecog1 <- 1*(lung$ph.ecog==1)
  > ecog2 <- 1*(lung$ph.ecog==2)
  > ecog3 <- 1*(lung$ph.ecog==3)
  > 
  > fit1 <- coxph(Surv(time, status) ~ age + ridge(ecog0, ecog1, ecog2, ecog3,
  +                                                scale=FALSE, theta=2), lung)
  > 
  > fit2 <- coxme(Surv(time, status) ~ age + (ecog0+ecog1+ecog2+ecog3 |1), lung,
  +               vfixed=.5)
  > 
  > aeq(fit1$coef, c(fixef(fit2), unlist(ranef(fit2))))
  [1] TRUE
  > indx <- c(5,1,2,3,4) #in coxme, shrinkage variables are first
  > all.equal(fit1$var, as.matrix(fit2$var)[indx, indx])
  [1] TRUE
  > 
  > fit3 <- coxme(Surv(time, status) ~ age + (1|ph.ecog), lung, vfixed=.5)
  > all.equal(fit2$var, fit3$var)
  [1] TRUE
  > all.equal(fit2$loglik, fit3$loglik)
  [1] TRUE
  > 
  > fit4 <- coxme(Surv(time, status) ~ age + (1|ph.ecog), lung)
  > dname <- paste("ecog", 0:3, sep='')
  > dummy <- matrix(diag(4), 4, dimnames=list(dname,dname))
  > fit5 <- coxme(Surv(time, status) ~ age + (ecog0+ecog1+ecog2+ecog3 |1), lung,
  +               varlist=dummy)
   ----------- FAILURE REPORT -------------- 
   --- failure: the condition has length > 1 ---
   --- srcref --- 
  : 
   --- package (from environment) --- 
  coxme
   --- call from context --- 
  coxme(Surv(time, status) ~ age + (ecog0 + ecog1 + ecog2 + ecog3 | 
      1), lung, varlist = dummy)
   --- call from argument --- 
  if (class(varlist) == "coxmevar") varlist <- list(varlist) else if (ismat(varlist)) varlist <- list(coxmeMlist(list(varlist))) else {
      if (!is.list(varlist)) 
          stop("Invalid varlist argument")
      if (all(sapply(varlist, ismat))) {
          if (nrandom > 1) 
              stop(paste("An unlabeled list of matrices is", "ambiguous when there are multiple random terms"))
          else varlist <- list(coxmeMlist(varlist))
      }
      else {
          for (i in 1:length(varlist)) {
              if (is.function(varlist[[i]])) 
                  varlist[[i]] <- varlist[[i]]()
              if (ismat(varlist[[i]])) 
                  varlist[[i]] <- coxmeMlist(list(varlist[[i]]))
              if (class(varlist[[i]]) != "coxmevar") {
                  if (is.list(varlist[[i]])) {
                    if (all(sapply(varlist[[i]], ismat))) 
                      varlist[[i]] <- coxmeMlist(varlist[[i]])
                    else stop("Invalid varlist element")
                  }
                  else stop("Invalid varlist element")
              }
          }
      }
  }
   --- R stacktrace ---
  where 1: coxme(Surv(time, status) ~ age + (ecog0 + ecog1 + ecog2 + ecog3 | 
      1), lung, varlist = dummy)
  
   --- value of length: 2 type: logical ---
  [1] FALSE FALSE
   --- function from context --- 
  function (formula, data, weights, subset, na.action, init, control, 
      ties = c("efron", "breslow"), varlist, vfixed, vinit, x = FALSE, 
      y = TRUE, refine.n = 0, random, fixed, variance, ...) 
  {
      ties <- match.arg(ties)
      Call <- match.call()
      if (!missing(fixed)) {
          if (missing(formula)) {
              formula <- fixed
              warning("The 'fixed' argument of coxme is depreciated")
          }
          else stop("Both a fixed and a formula argument are present")
      }
      if (!missing(random)) {
          warning("The random argument of coxme is depreciated")
          if (class(random) != "formula" || length(random) != 2) 
              stop("Invalid random formula")
          j <- length(formula)
          formula[[j]] <- call("+", formula[[j]], call("(", random[[2]]))
      }
      if (!missing(variance)) {
          warning("The variance argument of coxme is depreciated")
          vfixed <- variance
      }
      temp <- call("model.frame", formula = subbar(formula))
      for (i in c("data", "subset", "weights", "na.action")) if (!is.null(Call[[i]])) 
          temp[[i]] <- Call[[i]]
      if (is.R()) 
          m <- eval.parent(temp)
      else m <- eval(temp, sys.parent())
      Y <- model.extract(m, "response")
      n <- nrow(Y)
      if (n == 0) 
          stop("No observations remain in the data set")
      if (!inherits(Y, "Surv")) 
          stop("Response must be a survival object")
      type <- attr(Y, "type")
      if (type != "right" && type != "counting") 
          stop(paste("Cox model doesn't support '", type, "' survival data", 
              sep = ""))
      weights <- model.weights(m)
      if (length(weights) == 0) 
          weights <- rep(1, n)
      else if (any(weights <= 0)) 
          stop("Negative or zero weights are not allowed")
      offset <- model.offset(m)
      if (length(offset) == 0) 
          offset <- rep(0, n)
      pterms <- sapply(m, inherits, "coxph.penalty")
      if (any(pterms)) {
          stop("You cannot have penalized terms in coxme")
      }
      if (missing(control)) 
          control <- coxme.control(...)
      if (missing(init)) 
          init <- NULL
      flist <- formula1(formula)
      if (hasAbar(flist$fixed)) 
          stop("Invalid formula: a '|' outside of a valid random effects term")
      special <- c("strata", "cluster")
      Terms <- terms(flist$fixed, special)
      attr(Terms, "intercept") <- 1
      strats <- attr(Terms, "specials")$strata
      cluster <- attr(Terms, "specials")$cluster
      if (length(cluster)) {
          stop("A cluster() statement is invalid in coxme")
      }
      if (length(strats)) {
          temp <- untangle.specials(Terms, "strata", 1)
          if (length(temp$vars) == 1) 
              strata.keep <- m[[temp$vars]]
          else strata.keep <- strata(m[, temp$vars], shortlabel = T)
          strats <- as.numeric(strata.keep)
          X <- model.matrix(Terms[-temp$terms], m)[, -1, drop = F]
      }
      else X <- model.matrix(Terms, m)[, -1, drop = F]
      nrandom <- length(flist$random)
      if (nrandom == 0) 
          stop("No random effects terms found")
      vparm <- vector("list", nrandom)
      is.variance <- rep(TRUE, nrandom)
      ismat <- function(x) {
          inherits(x, "matrix") || inherits(x, "bdsmatrix") | inherits(x, 
              "Matrix")
      }
      if (missing(varlist) || is.null(varlist)) {
          varlist <- vector("list", nrandom)
          for (i in 1:nrandom) varlist[[i]] <- coxmeFull()
      }
      else {
          if (is.function(varlist)) 
              varlist <- varlist()
          if (class(varlist) == "coxmevar") 
              varlist <- list(varlist)
          else if (ismat(varlist)) 
              varlist <- list(coxmeMlist(list(varlist)))
          else {
              if (!is.list(varlist)) 
                  stop("Invalid varlist argument")
              if (all(sapply(varlist, ismat))) {
                  if (nrandom > 1) 
                    stop(paste("An unlabeled list of matrices is", 
                      "ambiguous when there are multiple random terms"))
                  else varlist <- list(coxmeMlist(varlist))
              }
              else {
                  for (i in 1:length(varlist)) {
                    if (is.function(varlist[[i]])) 
                      varlist[[i]] <- varlist[[i]]()
                    if (ismat(varlist[[i]])) 
                      varlist[[i]] <- coxmeMlist(list(varlist[[i]]))
                    if (class(varlist[[i]]) != "coxmevar") {
                      if (is.list(varlist[[i]])) {
                        if (all(sapply(varlist[[i]], ismat))) 
                          varlist[[i]] <- coxmeMlist(varlist[[i]])
                        else stop("Invalid varlist element")
                      }
                      else stop("Invalid varlist element")
                    }
                  }
              }
          }
          while (length(varlist) < nrandom) varlist <- c(varlist, 
              list(coxmeFull()))
      }
      if (!is.null(names(varlist))) {
          vname <- names(varlist)
          stop("Cannot (yet) have a names varlist")
          indx <- pmatch(vname, names(random), nomatch = 0)
          if (any(indx == 0 & vname != "")) 
              stop(paste("Varlist element not matched:", vname[indx == 
                  0 & vname != ""]))
          if (any(indx > 0)) {
              temp <- vector("list", nrandom)
              temp[indx] <- varlist[indx > 0]
              temp[-indx] <- varlist[indx == 0]
              varlist <- temp
          }
      }
      check <- sapply(varlist, function(x) {
          fname <- c("initialize", "generate", "wrapup")
          indx <- match(fname, names(x))
          if (any(is.na(x))) 
              stop(paste("Member not found in variance function:", 
                  fname(is.na(indx))))
          if (length(x) != 3 || any(!sapply(x, is.function))) 
              stop("Varlist objects must consist of exaclty three functions")
      })
      getcmat <- function(x, mf) {
          if (is.null(x) || x == 1) 
              return(NULL)
          Terms <- terms(eval(call("~", x)))
          attr(Terms, "intercept") <- 0
          varnames <- attr(Terms, "term.labels")
          ftemp <- sapply(mf[varnames], is.factor)
          if (any(ftemp)) {
              clist <- lapply(mf[varnames[ftemp]], function(x) diag(length(levels(x))))
              model.matrix(Terms, mf, contrasts.arg = clist)
          }
          else model.matrix(Terms, mf)
      }
      getGroupNames <- function(x) {
          if (is.call(x) && x[[1]] == as.name("/")) 
              c(getGroupNames(x[[2]]), getGroupNames(x[[3]]))
          else deparse(x)
      }
      getgroups <- function(x, mf) {
          if (is.numeric(x)) 
              return(NULL)
          varname <- getGroupNames(x)
          indx <- match(varname, names(mf), nomatch = 0)
          if (any(indx == 0)) 
              stop(paste("Invalid grouping factor", varname[indx == 
                  0]))
          else data.frame(lapply(mf[indx], as.factor))
      }
      if (missing(vinit) || is.null(vinit)) 
          vinit <- vector("list", nrandom)
      else {
          if (nrandom == 1) {
              if (is.numeric(vinit)) 
                  vinit <- list(vinit)
              else if (is.list(vinit)) 
                  vinit <- list(unlist(vinit))
          }
          if (!is.list(vinit)) 
              stop("Invalid value for `vinit` parameter")
          if (length(vinit) > nrandom) 
              stop(paste("Vinit must be a list of length", nrandom))
          if (!all(sapply(vinit, function(x) (is.null(x) || is.numeric(x))))) 
              stop("Vinit must contain numeric values")
          if (length(vinit) < nrandom) 
              vinit <- c(vinit, vector("list", nrandom - length(vinit)))
          tname <- names(vinit)
          if (!is.null(tname)) {
              stop("Named initial values not yet supported")
          }
      }
      if (missing(vfixed) || is.null(vfixed)) 
          vfixed <- vector("list", nrandom)
      else {
          if (nrandom == 1) {
              if (is.numeric(vfixed)) 
                  vfixed <- list(vfixed)
              else if (is.list(vfixed)) 
                  vfixed <- list(unlist(vfixed))
          }
          if (!is.list(vfixed)) 
              stop("Invalid value for `vfixed` parameter")
          if (length(vfixed) > nrandom) 
              stop(paste("Vfixed must be a list of length", nrandom))
          if (!all(sapply(vfixed, function(x) (is.null(x) || is.numeric(x))))) 
              stop("Vfixed must contain numeric values")
          if (length(vfixed) < nrandom) 
              vfixed <- c(vfixed, vector("list", nrandom - length(vfixed)))
          tname <- names(vfixed)
          if (!is.null(tname)) {
              temp <- pmatch(tname, names(flist$random), nomatch = 0)
              temp <- c(temp, (1:nrandom)[-temp])
              vfixed <- vfixed[temp]
          }
      }
      newzmat <- function(xmat, xmap) {
          n <- nrow(xmap)
          newz <- matrix(0, nrow = n, ncol = max(xmap))
          for (i in 1:ncol(xmap)) newz[cbind(1:n, xmap[, i])] <- xmat[, 
              i]
          newz
      }
      fmat <- zmat <- matrix(0, nrow = n, ncol = 0)
      ntheta <- integer(nrandom)
      ncoef <- matrix(0L, nrandom, 2, dimnames = list(NULL, c("intercept", 
          "slope")))
      itheta <- NULL
      for (i in 1:nrandom) {
          f2 <- formula2(flist$random[[i]])
          if (f2$intercept & f2$group == 1) 
              stop(paste("Error in random term ", i, ": Random intercepts require a grouping variable", 
                  sep = ""))
          vfun <- varlist[[i]]
          if (!is.null(f2$interaction)) 
              stop("Interactions not yet written")
          cmat <- getcmat(f2$fixed, m)
          groups <- getgroups(f2$group, m)
          ifun <- vfun$initialize(vinit[[i]], vfixed[[i]], intercept = f2$intercept, 
              groups, cmat, control)
          if (!is.null(ifun$error)) 
              stop(paste("In random term ", i, ": ", ifun$error, 
                  sep = ""))
          vparm[[i]] <- ifun$parms
          if (!is.null(ifun$is.variance)) 
              is.variance[i] <- ifun$is.variance
          itheta <- c(itheta, ifun$theta)
          ntheta[i] <- length(ifun$theta)
          if (f2$intercept) {
              if (!is.matrix(ifun$imap) || nrow(ifun$imap) != n) 
                  stop(paste("In random term ", i, ": Invalid intercept matrix F", 
                    sep = ""))
              temp <- sort(unique(c(ifun$imap)))
              if (any(temp != 1:length(temp))) 
                  stop(paste("In random term ", i, ": intercept matrix has an invalid element", 
                    sep = ""))
              if (ncol(fmat) > 0) 
                  fmat <- cbind(fmat, ifun$imap + max(fmat))
              else fmat <- ifun$imap
              ncoef[i, 1] <- 1 + max(ifun$imap) - min(ifun$imap)
          }
          if (length(cmat) > 0) {
              if (is.null(ifun$xmap) || is.null(ifun$X) || !is.matrix(ifun$xmap) || 
                  !is.matrix(ifun$X) || nrow(ifun$xmap) != n || 
                  nrow(ifun$X) != n || ncol(ifun$xmap) != ncol(ifun$X)) 
                  stop(paste("In random term ", i, "invalid X/xmap pair"))
              if (f2$intercept) 
                  xmap <- ifun$xmap - max(ifun$imap)
              else xmap <- ifun$xmap
              if (any(sort(unique(c(xmap))) != 1:max(xmap))) 
                  stop(paste("In random term ", i, ": xmap matrix has an invalid element", 
                    sep = ""))
              temp <- newzmat(ifun$X, xmap)
              ncoef[i, 2] <- ncol(temp)
              zmat <- cbind(zmat, temp)
          }
      }
      if (any(is.variance) & !all(is.variance)) 
          stop("All variance terms must have the same is.variance setting")
      fit <- coxme.fit(X, Y, strats, offset, init, control, weights = weights, 
          ties = ties, row.names(m), fmat, zmat, varlist, vparm, 
          itheta, ntheta, ncoef, refine.n, is.variance = any(is.variance))
      if (is.character(fit)) {
          fit <- list(fail = fit)
          oldClass(fit) <- "coxme"
          return(fit)
      }
      fcoef <- fit$coefficients
      nvar <- length(fcoef)
      if (length(fcoef) > 0 && any(is.na(fcoef))) {
          vars <- (1:length(fcoef))[is.na(fcoef)]
          msg <- paste("X matrix deemed to be singular; variable", 
              paste(vars, collapse = " "))
          warning(msg)
      }
      if (length(fcoef) > 0) {
          names(fit$coefficients) <- dimnames(X)[[2]]
      }
      rlinear <- rep(0, nrow(Y))
      if (length(fmat)) {
          for (i in 1:ncol(fmat)) {
              rlinear <- rlinear + fit$frail[fmat[, i]]
          }
      }
      if (length(zmat)) {
          indx <- if (length(fmat) > 0) 
              max(fmat)
          else 0
          for (i in 1:ncol(zmat)) rlinear <- rlinear + fit$frail[indx + 
              i] * zmat[, i]
      }
      if (nvar == 0) 
          fit$linear.predictor <- rlinear
      else fit$linear.predictor <- as.vector(rlinear + c(X %*% 
          fit$coef))
      newtheta <- random.coef <- list()
      nrandom <- length(varlist)
      sindex <- rep(1:nrandom, ntheta)
      bindex <- rep(row(ncoef), ncoef)
      for (i in 1:nrandom) {
          temp <- varlist[[i]]$wrapup(fit$theta[sindex == i], fit$frail[bindex == 
              i], vparm[[i]])
          newtheta <- c(newtheta, temp$theta)
          if (!is.list(temp$b)) {
              temp$b <- list(temp$b)
              names(temp$b) <- paste("Random", i, sep = "")
          }
          random.coef <- c(random.coef, temp$b)
      }
      fit$frail <- random.coef
      fit$vcoef <- newtheta
      fit$theta <- NULL
      fit$n <- c(sum(Y[, ncol(Y)]), nrow(Y))
      fit$terms <- Terms
      fit$assign <- attr(X, "assign")
      fit$formulaList <- flist
      na.action <- attr(m, "na.action")
      if (length(na.action)) 
          fit$na.action <- na.action
      if (x) {
          fit$x <- X
          if (length(strats)) 
              fit$strata <- strata.keep
      }
      if (y) 
          fit$y <- Y
      if (!is.null(weights) && any(weights != 1)) 
          fit$weights <- weights
      fit$formula <- as.vector(attr(Terms, "formula"))
      fit$call <- Call
      fit$ties <- ties
      names(fit$loglik) <- c("NULL", "Integrated", "Penalized")
      oldClass(fit) <- "coxme"
      fit
  }
  <bytecode: 0xac78428>
  <environment: namespace:coxme>
   --- function search by body ---
  Function coxme in namespace coxme has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Error in if (class(varlist) == "coxmevar") varlist <- list(varlist) else if (ismat(varlist)) varlist <- list(coxmeMlist(list(varlist))) else { : 
    the condition has length > 1
  Calls: coxme
  Execution halted
Running the tests in 'tests/slope2.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > options(na.action='na.exclude', contrasts=c('contr.treatment', 'contr.poly'))
  > aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > 
  > #
  > # Test of fitting random slopes
  > #
  > # Same data set as slope1, refit using variance matrices
  > #
  > set.seed(56)
  > n.subject <- seq(180, by=21, length=9) # number of subjects
  > slope <- sort(-.5 + rnorm(9, sd=.5))         # true treament effects
  > 
  > inst <- rep(1:9, n.subject)
  > n <- length(inst)
  > simdata <- data.frame(id=1:n, inst=inst,
  +                       trt= rep(0:1, length=n),
  +                       age= runif(n, 40, 70))
  > #risk goes up 30%/decade of age
  > simdata$hazard <- .8* exp(simdata$trt * rep(slope, n.subject) +
  +                           (simdata$age-55) * .03)
  > 
  > rtime <- function(hazard, censor=c(1,2)) {
  +     stime <- rexp(length(hazard), rate=hazard)
  +     ctime <- runif(length(hazard), censor[1], censor[2])
  +     list(time= pmin(stime, ctime), status=1*(stime <=ctime))
  +     }
  > temp <- rtime(simdata$hazard)
  > simdata$time <- temp$time
  > simdata$status <- temp$status
  > 
  > #Check fit1
  > fit1 <- coxme(Surv(time, status) ~ age + trt + (1|inst), simdata)
  > fit1b <- coxme(Surv(time, status) ~ age + trt + (1|inst), simdata,
  +                varlist=diag(9))
   ----------- FAILURE REPORT -------------- 
   --- failure: the condition has length > 1 ---
   --- srcref --- 
  : 
   --- package (from environment) --- 
  coxme
   --- call from context --- 
  coxme(Surv(time, status) ~ age + trt + (1 | inst), simdata, varlist = diag(9))
   --- call from argument --- 
  if (class(varlist) == "coxmevar") varlist <- list(varlist) else if (ismat(varlist)) varlist <- list(coxmeMlist(list(varlist))) else {
      if (!is.list(varlist)) 
          stop("Invalid varlist argument")
      if (all(sapply(varlist, ismat))) {
          if (nrandom > 1) 
              stop(paste("An unlabeled list of matrices is", "ambiguous when there are multiple random terms"))
          else varlist <- list(coxmeMlist(varlist))
      }
      else {
          for (i in 1:length(varlist)) {
              if (is.function(varlist[[i]])) 
                  varlist[[i]] <- varlist[[i]]()
              if (ismat(varlist[[i]])) 
                  varlist[[i]] <- coxmeMlist(list(varlist[[i]]))
              if (class(varlist[[i]]) != "coxmevar") {
                  if (is.list(varlist[[i]])) {
                    if (all(sapply(varlist[[i]], ismat))) 
                      varlist[[i]] <- coxmeMlist(varlist[[i]])
                    else stop("Invalid varlist element")
                  }
                  else stop("Invalid varlist element")
              }
          }
      }
  }
   --- R stacktrace ---
  where 1: coxme(Surv(time, status) ~ age + trt + (1 | inst), simdata, varlist = diag(9))
  
   --- value of length: 2 type: logical ---
  [1] FALSE FALSE
   --- function from context --- 
  function (formula, data, weights, subset, na.action, init, control, 
      ties = c("efron", "breslow"), varlist, vfixed, vinit, x = FALSE, 
      y = TRUE, refine.n = 0, random, fixed, variance, ...) 
  {
      ties <- match.arg(ties)
      Call <- match.call()
      if (!missing(fixed)) {
          if (missing(formula)) {
              formula <- fixed
              warning("The 'fixed' argument of coxme is depreciated")
          }
          else stop("Both a fixed and a formula argument are present")
      }
      if (!missing(random)) {
          warning("The random argument of coxme is depreciated")
          if (class(random) != "formula" || length(random) != 2) 
              stop("Invalid random formula")
          j <- length(formula)
          formula[[j]] <- call("+", formula[[j]], call("(", random[[2]]))
      }
      if (!missing(variance)) {
          warning("The variance argument of coxme is depreciated")
          vfixed <- variance
      }
      temp <- call("model.frame", formula = subbar(formula))
      for (i in c("data", "subset", "weights", "na.action")) if (!is.null(Call[[i]])) 
          temp[[i]] <- Call[[i]]
      if (is.R()) 
          m <- eval.parent(temp)
      else m <- eval(temp, sys.parent())
      Y <- model.extract(m, "response")
      n <- nrow(Y)
      if (n == 0) 
          stop("No observations remain in the data set")
      if (!inherits(Y, "Surv")) 
          stop("Response must be a survival object")
      type <- attr(Y, "type")
      if (type != "right" && type != "counting") 
          stop(paste("Cox model doesn't support '", type, "' survival data", 
              sep = ""))
      weights <- model.weights(m)
      if (length(weights) == 0) 
          weights <- rep(1, n)
      else if (any(weights <= 0)) 
          stop("Negative or zero weights are not allowed")
      offset <- model.offset(m)
      if (length(offset) == 0) 
          offset <- rep(0, n)
      pterms <- sapply(m, inherits, "coxph.penalty")
      if (any(pterms)) {
          stop("You cannot have penalized terms in coxme")
      }
      if (missing(control)) 
          control <- coxme.control(...)
      if (missing(init)) 
          init <- NULL
      flist <- formula1(formula)
      if (hasAbar(flist$fixed)) 
          stop("Invalid formula: a '|' outside of a valid random effects term")
      special <- c("strata", "cluster")
      Terms <- terms(flist$fixed, special)
      attr(Terms, "intercept") <- 1
      strats <- attr(Terms, "specials")$strata
      cluster <- attr(Terms, "specials")$cluster
      if (length(cluster)) {
          stop("A cluster() statement is invalid in coxme")
      }
      if (length(strats)) {
          temp <- untangle.specials(Terms, "strata", 1)
          if (length(temp$vars) == 1) 
              strata.keep <- m[[temp$vars]]
          else strata.keep <- strata(m[, temp$vars], shortlabel = T)
          strats <- as.numeric(strata.keep)
          X <- model.matrix(Terms[-temp$terms], m)[, -1, drop = F]
      }
      else X <- model.matrix(Terms, m)[, -1, drop = F]
      nrandom <- length(flist$random)
      if (nrandom == 0) 
          stop("No random effects terms found")
      vparm <- vector("list", nrandom)
      is.variance <- rep(TRUE, nrandom)
      ismat <- function(x) {
          inherits(x, "matrix") || inherits(x, "bdsmatrix") | inherits(x, 
              "Matrix")
      }
      if (missing(varlist) || is.null(varlist)) {
          varlist <- vector("list", nrandom)
          for (i in 1:nrandom) varlist[[i]] <- coxmeFull()
      }
      else {
          if (is.function(varlist)) 
              varlist <- varlist()
          if (class(varlist) == "coxmevar") 
              varlist <- list(varlist)
          else if (ismat(varlist)) 
              varlist <- list(coxmeMlist(list(varlist)))
          else {
              if (!is.list(varlist)) 
                  stop("Invalid varlist argument")
              if (all(sapply(varlist, ismat))) {
                  if (nrandom > 1) 
                    stop(paste("An unlabeled list of matrices is", 
                      "ambiguous when there are multiple random terms"))
                  else varlist <- list(coxmeMlist(varlist))
              }
              else {
                  for (i in 1:length(varlist)) {
                    if (is.function(varlist[[i]])) 
                      varlist[[i]] <- varlist[[i]]()
                    if (ismat(varlist[[i]])) 
                      varlist[[i]] <- coxmeMlist(list(varlist[[i]]))
                    if (class(varlist[[i]]) != "coxmevar") {
                      if (is.list(varlist[[i]])) {
                        if (all(sapply(varlist[[i]], ismat))) 
                          varlist[[i]] <- coxmeMlist(varlist[[i]])
                        else stop("Invalid varlist element")
                      }
                      else stop("Invalid varlist element")
                    }
                  }
              }
          }
          while (length(varlist) < nrandom) varlist <- c(varlist, 
              list(coxmeFull()))
      }
      if (!is.null(names(varlist))) {
          vname <- names(varlist)
          stop("Cannot (yet) have a names varlist")
          indx <- pmatch(vname, names(random), nomatch = 0)
          if (any(indx == 0 & vname != "")) 
              stop(paste("Varlist element not matched:", vname[indx == 
                  0 & vname != ""]))
          if (any(indx > 0)) {
              temp <- vector("list", nrandom)
              temp[indx] <- varlist[indx > 0]
              temp[-indx] <- varlist[indx == 0]
              varlist <- temp
          }
      }
      check <- sapply(varlist, function(x) {
          fname <- c("initialize", "generate", "wrapup")
          indx <- match(fname, names(x))
          if (any(is.na(x))) 
              stop(paste("Member not found in variance function:", 
                  fname(is.na(indx))))
          if (length(x) != 3 || any(!sapply(x, is.function))) 
              stop("Varlist objects must consist of exaclty three functions")
      })
      getcmat <- function(x, mf) {
          if (is.null(x) || x == 1) 
              return(NULL)
          Terms <- terms(eval(call("~", x)))
          attr(Terms, "intercept") <- 0
          varnames <- attr(Terms, "term.labels")
          ftemp <- sapply(mf[varnames], is.factor)
          if (any(ftemp)) {
              clist <- lapply(mf[varnames[ftemp]], function(x) diag(length(levels(x))))
              model.matrix(Terms, mf, contrasts.arg = clist)
          }
          else model.matrix(Terms, mf)
      }
      getGroupNames <- function(x) {
          if (is.call(x) && x[[1]] == as.name("/")) 
              c(getGroupNames(x[[2]]), getGroupNames(x[[3]]))
          else deparse(x)
      }
      getgroups <- function(x, mf) {
          if (is.numeric(x)) 
              return(NULL)
          varname <- getGroupNames(x)
          indx <- match(varname, names(mf), nomatch = 0)
          if (any(indx == 0)) 
              stop(paste("Invalid grouping factor", varname[indx == 
                  0]))
          else data.frame(lapply(mf[indx], as.factor))
      }
      if (missing(vinit) || is.null(vinit)) 
          vinit <- vector("list", nrandom)
      else {
          if (nrandom == 1) {
              if (is.numeric(vinit)) 
                  vinit <- list(vinit)
              else if (is.list(vinit)) 
                  vinit <- list(unlist(vinit))
          }
          if (!is.list(vinit)) 
              stop("Invalid value for `vinit` parameter")
          if (length(vinit) > nrandom) 
              stop(paste("Vinit must be a list of length", nrandom))
          if (!all(sapply(vinit, function(x) (is.null(x) || is.numeric(x))))) 
              stop("Vinit must contain numeric values")
          if (length(vinit) < nrandom) 
              vinit <- c(vinit, vector("list", nrandom - length(vinit)))
          tname <- names(vinit)
          if (!is.null(tname)) {
              stop("Named initial values not yet supported")
          }
      }
      if (missing(vfixed) || is.null(vfixed)) 
          vfixed <- vector("list", nrandom)
      else {
          if (nrandom == 1) {
              if (is.numeric(vfixed)) 
                  vfixed <- list(vfixed)
              else if (is.list(vfixed)) 
                  vfixed <- list(unlist(vfixed))
          }
          if (!is.list(vfixed)) 
              stop("Invalid value for `vfixed` parameter")
          if (length(vfixed) > nrandom) 
              stop(paste("Vfixed must be a list of length", nrandom))
          if (!all(sapply(vfixed, function(x) (is.null(x) || is.numeric(x))))) 
              stop("Vfixed must contain numeric values")
          if (length(vfixed) < nrandom) 
              vfixed <- c(vfixed, vector("list", nrandom - length(vfixed)))
          tname <- names(vfixed)
          if (!is.null(tname)) {
              temp <- pmatch(tname, names(flist$random), nomatch = 0)
              temp <- c(temp, (1:nrandom)[-temp])
              vfixed <- vfixed[temp]
          }
      }
      newzmat <- function(xmat, xmap) {
          n <- nrow(xmap)
          newz <- matrix(0, nrow = n, ncol = max(xmap))
          for (i in 1:ncol(xmap)) newz[cbind(1:n, xmap[, i])] <- xmat[, 
              i]
          newz
      }
      fmat <- zmat <- matrix(0, nrow = n, ncol = 0)
      ntheta <- integer(nrandom)
      ncoef <- matrix(0L, nrandom, 2, dimnames = list(NULL, c("intercept", 
          "slope")))
      itheta <- NULL
      for (i in 1:nrandom) {
          f2 <- formula2(flist$random[[i]])
          if (f2$intercept & f2$group == 1) 
              stop(paste("Error in random term ", i, ": Random intercepts require a grouping variable", 
                  sep = ""))
          vfun <- varlist[[i]]
          if (!is.null(f2$interaction)) 
              stop("Interactions not yet written")
          cmat <- getcmat(f2$fixed, m)
          groups <- getgroups(f2$group, m)
          ifun <- vfun$initialize(vinit[[i]], vfixed[[i]], intercept = f2$intercept, 
              groups, cmat, control)
          if (!is.null(ifun$error)) 
              stop(paste("In random term ", i, ": ", ifun$error, 
                  sep = ""))
          vparm[[i]] <- ifun$parms
          if (!is.null(ifun$is.variance)) 
              is.variance[i] <- ifun$is.variance
          itheta <- c(itheta, ifun$theta)
          ntheta[i] <- length(ifun$theta)
          if (f2$intercept) {
              if (!is.matrix(ifun$imap) || nrow(ifun$imap) != n) 
                  stop(paste("In random term ", i, ": Invalid intercept matrix F", 
                    sep = ""))
              temp <- sort(unique(c(ifun$imap)))
              if (any(temp != 1:length(temp))) 
                  stop(paste("In random term ", i, ": intercept matrix has an invalid element", 
                    sep = ""))
              if (ncol(fmat) > 0) 
                  fmat <- cbind(fmat, ifun$imap + max(fmat))
              else fmat <- ifun$imap
              ncoef[i, 1] <- 1 + max(ifun$imap) - min(ifun$imap)
          }
          if (length(cmat) > 0) {
              if (is.null(ifun$xmap) || is.null(ifun$X) || !is.matrix(ifun$xmap) || 
                  !is.matrix(ifun$X) || nrow(ifun$xmap) != n || 
                  nrow(ifun$X) != n || ncol(ifun$xmap) != ncol(ifun$X)) 
                  stop(paste("In random term ", i, "invalid X/xmap pair"))
              if (f2$intercept) 
                  xmap <- ifun$xmap - max(ifun$imap)
              else xmap <- ifun$xmap
              if (any(sort(unique(c(xmap))) != 1:max(xmap))) 
                  stop(paste("In random term ", i, ": xmap matrix has an invalid element", 
                    sep = ""))
              temp <- newzmat(ifun$X, xmap)
              ncoef[i, 2] <- ncol(temp)
              zmat <- cbind(zmat, temp)
          }
      }
      if (any(is.variance) & !all(is.variance)) 
          stop("All variance terms must have the same is.variance setting")
      fit <- coxme.fit(X, Y, strats, offset, init, control, weights = weights, 
          ties = ties, row.names(m), fmat, zmat, varlist, vparm, 
          itheta, ntheta, ncoef, refine.n, is.variance = any(is.variance))
      if (is.character(fit)) {
          fit <- list(fail = fit)
          oldClass(fit) <- "coxme"
          return(fit)
      }
      fcoef <- fit$coefficients
      nvar <- length(fcoef)
      if (length(fcoef) > 0 && any(is.na(fcoef))) {
          vars <- (1:length(fcoef))[is.na(fcoef)]
          msg <- paste("X matrix deemed to be singular; variable", 
              paste(vars, collapse = " "))
          warning(msg)
      }
      if (length(fcoef) > 0) {
          names(fit$coefficients) <- dimnames(X)[[2]]
      }
      rlinear <- rep(0, nrow(Y))
      if (length(fmat)) {
          for (i in 1:ncol(fmat)) {
              rlinear <- rlinear + fit$frail[fmat[, i]]
          }
      }
      if (length(zmat)) {
          indx <- if (length(fmat) > 0) 
              max(fmat)
          else 0
          for (i in 1:ncol(zmat)) rlinear <- rlinear + fit$frail[indx + 
              i] * zmat[, i]
      }
      if (nvar == 0) 
          fit$linear.predictor <- rlinear
      else fit$linear.predictor <- as.vector(rlinear + c(X %*% 
          fit$coef))
      newtheta <- random.coef <- list()
      nrandom <- length(varlist)
      sindex <- rep(1:nrandom, ntheta)
      bindex <- rep(row(ncoef), ncoef)
      for (i in 1:nrandom) {
          temp <- varlist[[i]]$wrapup(fit$theta[sindex == i], fit$frail[bindex == 
              i], vparm[[i]])
          newtheta <- c(newtheta, temp$theta)
          if (!is.list(temp$b)) {
              temp$b <- list(temp$b)
              names(temp$b) <- paste("Random", i, sep = "")
          }
          random.coef <- c(random.coef, temp$b)
      }
      fit$frail <- random.coef
      fit$vcoef <- newtheta
      fit$theta <- NULL
      fit$n <- c(sum(Y[, ncol(Y)]), nrow(Y))
      fit$terms <- Terms
      fit$assign <- attr(X, "assign")
      fit$formulaList <- flist
      na.action <- attr(m, "na.action")
      if (length(na.action)) 
          fit$na.action <- na.action
      if (x) {
          fit$x <- X
          if (length(strats)) 
              fit$strata <- strata.keep
      }
      if (y) 
          fit$y <- Y
      if (!is.null(weights) && any(weights != 1)) 
          fit$weights <- weights
      fit$formula <- as.vector(attr(Terms, "formula"))
      fit$call <- Call
      fit$ties <- ties
      names(fit$loglik) <- c("NULL", "Integrated", "Penalized")
      oldClass(fit) <- "coxme"
      fit
  }
  <bytecode: 0x92adeb0>
  <environment: namespace:coxme>
   --- function search by body ---
  Function coxme in namespace coxme has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Error in if (class(varlist) == "coxmevar") varlist <- list(varlist) else if (ismat(varlist)) varlist <- list(coxmeMlist(list(varlist))) else { : 
    the condition has length > 1
  Calls: coxme
  Execution halted
Running the tests in 'tests/test0.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > options(na.action='na.exclude', contrasts=c('contr.treatment', 'contr.poly'))
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > 
  > # Really simple dataset -- covariate x1 is our old friend from the 
  > #  validation section at the back of the book.  (Well, as yet only
  > #  my online copy has the examples with weights).
  > #
  > tdata0 <- data.frame(time  =c(5,4,1,1,2,2,2,2,3), 
  +                      status=c(0,1,1,0,1,1,1,0,0),
  +                      x1    =c(0,1,2,0,1,1,0,1,0),
  +                      wt    =c(1,2,1,2,3,4,3,2,1),
  +                      x2    =c(1,3,5,2,3,6,4,3,1),
  +                      grp   =c(1,1,2,2,1,1,2,2,1))
  > 
  > # these 3 functions give the loglik, and the u/imat results for variable
  > #   x1
  > lfun <- function(beta, efron=T) {
  +     r <- exp(beta)
  +     a <- 7*r +3
  +     b <- 4*r +2
  +     temp1 <- 11*beta - (log(r^2 + 11*r +7) + 2*log(2*r +1))
  +     if (efron) temp2 <- (10/3)*(log(a+b) + log(a*2/3 +b) + log(a/3 +b))
  +     else       temp2 <- 10 * log(a+b)
  +     temp1 - temp2
  +     }
  > 
  > ufun <- function(beta, efron=T) {
  +     r <- exp(beta)
  +     a <- 7*r +3
  +     b <- 4*r +2
  +     xbar1 <- (2*r^2 + 11*r)/(r^2 + 11*r + 7)
  +     xbar2 <- 11*r/(11*r +5)
  +     xbar3 <- 2*r/(2*r +1)
  +     xbar2b <- (7*r*2/3 + 4*r)/(a*2/3 +b)
  +     xbar2c <- (7*r/3 + 4*r)/(a/3 + b)
  + 
  +     temp1 <- 11 - (xbar1 + 2*xbar3)
  +     if (efron) temp2 <- (10/3)*(xbar2 + xbar2b + xbar2c)
  +     else       temp2 <- 10*xbar2
  +     temp1 - temp2
  +     }
  > 
  > ifun <- function(beta, efron=T) {
  +     r <- exp(beta)
  +     a <- 7*r +3
  +     b <- 4*r +2
  +     xbar1 <- (2*r^2 + 11*r)/(r^2 + 11*r + 7)
  +     xbar2 <- 11*r/(11*r +5)
  +     xbar3 <- 2*r/(2*r +1)
  +     xbar2b <- (7*r*2/3 + 4*r)/(a*2/3 +b)
  +     xbar2c <- (7*r/3 + 4*r)/(a/3 + b)
  + 
  +     temp1 <- (4*r^2 + 11*r)/(r^2 + 11*r +7) - xbar1^2
  +     if (efron) temp2 <- (10/3)*((xbar2 - xbar2^2) + (xbar2b - xbar2b^2) +
  +                                 (xbar2c- xbar2c^2))
  +     else       temp2 <- 10*(xbar2 - xbar2^2)
  +     temp1 + temp2 + 2*(xbar3 - xbar3^2)
  +     }
  > 
  > tfit <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
  +               vfixed=.5, weight=wt, iter=0,
  +               ties='breslow', sparse.calc=0)
  > 
  > aeq(tfit$loglik[1], lfun(0,F))
  [1] TRUE
  > aeq(tfit$u[3], ufun(0,F))
  [1] TRUE
  > aeq((solve(tfit$var))[3,3], ifun(0,F))
  [1] TRUE
  > 
  > tfit1 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
  +               vfixed=.5, weight=wt, iter=0,
  +               ties='breslow', sparse.calc=1)
  > aeq(tfit$u, tfit1$u)
  [1] TRUE
  > all.equal(tfit$var, tfit1$var)
  [1] TRUE
  > aeq(tfit$loglik, tfit1$loglik)
  [1] TRUE
  > 
  > # Do the matrix form, using coxmeMlist
  > dmat <- diag(2)
  > dimnames(dmat) <- list(1:2, 1:2)
  > tfit2 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
  +               vfixed=.5, weight=wt, iter=0,
  +               ties='breslow', varlist=dmat)
   ----------- FAILURE REPORT -------------- 
   --- failure: the condition has length > 1 ---
   --- srcref --- 
  : 
   --- package (from environment) --- 
  coxme
   --- call from context --- 
  coxme(Surv(time, status) ~ x1 + x2 + (1 | grp), data = tdata0, 
      vfixed = 0.5, weight = wt, iter = 0, ties = "breslow", varlist = dmat)
   --- call from argument --- 
  if (class(varlist) == "coxmevar") varlist <- list(varlist) else if (ismat(varlist)) varlist <- list(coxmeMlist(list(varlist))) else {
      if (!is.list(varlist)) 
          stop("Invalid varlist argument")
      if (all(sapply(varlist, ismat))) {
          if (nrandom > 1) 
              stop(paste("An unlabeled list of matrices is", "ambiguous when there are multiple random terms"))
          else varlist <- list(coxmeMlist(varlist))
      }
      else {
          for (i in 1:length(varlist)) {
              if (is.function(varlist[[i]])) 
                  varlist[[i]] <- varlist[[i]]()
              if (ismat(varlist[[i]])) 
                  varlist[[i]] <- coxmeMlist(list(varlist[[i]]))
              if (class(varlist[[i]]) != "coxmevar") {
                  if (is.list(varlist[[i]])) {
                    if (all(sapply(varlist[[i]], ismat))) 
                      varlist[[i]] <- coxmeMlist(varlist[[i]])
                    else stop("Invalid varlist element")
                  }
                  else stop("Invalid varlist element")
              }
          }
      }
  }
   --- R stacktrace ---
  where 1: coxme(Surv(time, status) ~ x1 + x2 + (1 | grp), data = tdata0, 
      vfixed = 0.5, weight = wt, iter = 0, ties = "breslow", varlist = dmat)
  
   --- value of length: 2 type: logical ---
  [1] FALSE FALSE
   --- function from context --- 
  function (formula, data, weights, subset, na.action, init, control, 
      ties = c("efron", "breslow"), varlist, vfixed, vinit, x = FALSE, 
      y = TRUE, refine.n = 0, random, fixed, variance, ...) 
  {
      ties <- match.arg(ties)
      Call <- match.call()
      if (!missing(fixed)) {
          if (missing(formula)) {
              formula <- fixed
              warning("The 'fixed' argument of coxme is depreciated")
          }
          else stop("Both a fixed and a formula argument are present")
      }
      if (!missing(random)) {
          warning("The random argument of coxme is depreciated")
          if (class(random) != "formula" || length(random) != 2) 
              stop("Invalid random formula")
          j <- length(formula)
          formula[[j]] <- call("+", formula[[j]], call("(", random[[2]]))
      }
      if (!missing(variance)) {
          warning("The variance argument of coxme is depreciated")
          vfixed <- variance
      }
      temp <- call("model.frame", formula = subbar(formula))
      for (i in c("data", "subset", "weights", "na.action")) if (!is.null(Call[[i]])) 
          temp[[i]] <- Call[[i]]
      if (is.R()) 
          m <- eval.parent(temp)
      else m <- eval(temp, sys.parent())
      Y <- model.extract(m, "response")
      n <- nrow(Y)
      if (n == 0) 
          stop("No observations remain in the data set")
      if (!inherits(Y, "Surv")) 
          stop("Response must be a survival object")
      type <- attr(Y, "type")
      if (type != "right" && type != "counting") 
          stop(paste("Cox model doesn't support '", type, "' survival data", 
              sep = ""))
      weights <- model.weights(m)
      if (length(weights) == 0) 
          weights <- rep(1, n)
      else if (any(weights <= 0)) 
          stop("Negative or zero weights are not allowed")
      offset <- model.offset(m)
      if (length(offset) == 0) 
          offset <- rep(0, n)
      pterms <- sapply(m, inherits, "coxph.penalty")
      if (any(pterms)) {
          stop("You cannot have penalized terms in coxme")
      }
      if (missing(control)) 
          control <- coxme.control(...)
      if (missing(init)) 
          init <- NULL
      flist <- formula1(formula)
      if (hasAbar(flist$fixed)) 
          stop("Invalid formula: a '|' outside of a valid random effects term")
      special <- c("strata", "cluster")
      Terms <- terms(flist$fixed, special)
      attr(Terms, "intercept") <- 1
      strats <- attr(Terms, "specials")$strata
      cluster <- attr(Terms, "specials")$cluster
      if (length(cluster)) {
          stop("A cluster() statement is invalid in coxme")
      }
      if (length(strats)) {
          temp <- untangle.specials(Terms, "strata", 1)
          if (length(temp$vars) == 1) 
              strata.keep <- m[[temp$vars]]
          else strata.keep <- strata(m[, temp$vars], shortlabel = T)
          strats <- as.numeric(strata.keep)
          X <- model.matrix(Terms[-temp$terms], m)[, -1, drop = F]
      }
      else X <- model.matrix(Terms, m)[, -1, drop = F]
      nrandom <- length(flist$random)
      if (nrandom == 0) 
          stop("No random effects terms found")
      vparm <- vector("list", nrandom)
      is.variance <- rep(TRUE, nrandom)
      ismat <- function(x) {
          inherits(x, "matrix") || inherits(x, "bdsmatrix") | inherits(x, 
              "Matrix")
      }
      if (missing(varlist) || is.null(varlist)) {
          varlist <- vector("list", nrandom)
          for (i in 1:nrandom) varlist[[i]] <- coxmeFull()
      }
      else {
          if (is.function(varlist)) 
              varlist <- varlist()
          if (class(varlist) == "coxmevar") 
              varlist <- list(varlist)
          else if (ismat(varlist)) 
              varlist <- list(coxmeMlist(list(varlist)))
          else {
              if (!is.list(varlist)) 
                  stop("Invalid varlist argument")
              if (all(sapply(varlist, ismat))) {
                  if (nrandom > 1) 
                    stop(paste("An unlabeled list of matrices is", 
                      "ambiguous when there are multiple random terms"))
                  else varlist <- list(coxmeMlist(varlist))
              }
              else {
                  for (i in 1:length(varlist)) {
                    if (is.function(varlist[[i]])) 
                      varlist[[i]] <- varlist[[i]]()
                    if (ismat(varlist[[i]])) 
                      varlist[[i]] <- coxmeMlist(list(varlist[[i]]))
                    if (class(varlist[[i]]) != "coxmevar") {
                      if (is.list(varlist[[i]])) {
                        if (all(sapply(varlist[[i]], ismat))) 
                          varlist[[i]] <- coxmeMlist(varlist[[i]])
                        else stop("Invalid varlist element")
                      }
                      else stop("Invalid varlist element")
                    }
                  }
              }
          }
          while (length(varlist) < nrandom) varlist <- c(varlist, 
              list(coxmeFull()))
      }
      if (!is.null(names(varlist))) {
          vname <- names(varlist)
          stop("Cannot (yet) have a names varlist")
          indx <- pmatch(vname, names(random), nomatch = 0)
          if (any(indx == 0 & vname != "")) 
              stop(paste("Varlist element not matched:", vname[indx == 
                  0 & vname != ""]))
          if (any(indx > 0)) {
              temp <- vector("list", nrandom)
              temp[indx] <- varlist[indx > 0]
              temp[-indx] <- varlist[indx == 0]
              varlist <- temp
          }
      }
      check <- sapply(varlist, function(x) {
          fname <- c("initialize", "generate", "wrapup")
          indx <- match(fname, names(x))
          if (any(is.na(x))) 
              stop(paste("Member not found in variance function:", 
                  fname(is.na(indx))))
          if (length(x) != 3 || any(!sapply(x, is.function))) 
              stop("Varlist objects must consist of exaclty three functions")
      })
      getcmat <- function(x, mf) {
          if (is.null(x) || x == 1) 
              return(NULL)
          Terms <- terms(eval(call("~", x)))
          attr(Terms, "intercept") <- 0
          varnames <- attr(Terms, "term.labels")
          ftemp <- sapply(mf[varnames], is.factor)
          if (any(ftemp)) {
              clist <- lapply(mf[varnames[ftemp]], function(x) diag(length(levels(x))))
              model.matrix(Terms, mf, contrasts.arg = clist)
          }
          else model.matrix(Terms, mf)
      }
      getGroupNames <- function(x) {
          if (is.call(x) && x[[1]] == as.name("/")) 
              c(getGroupNames(x[[2]]), getGroupNames(x[[3]]))
          else deparse(x)
      }
      getgroups <- function(x, mf) {
          if (is.numeric(x)) 
              return(NULL)
          varname <- getGroupNames(x)
          indx <- match(varname, names(mf), nomatch = 0)
          if (any(indx == 0)) 
              stop(paste("Invalid grouping factor", varname[indx == 
                  0]))
          else data.frame(lapply(mf[indx], as.factor))
      }
      if (missing(vinit) || is.null(vinit)) 
          vinit <- vector("list", nrandom)
      else {
          if (nrandom == 1) {
              if (is.numeric(vinit)) 
                  vinit <- list(vinit)
              else if (is.list(vinit)) 
                  vinit <- list(unlist(vinit))
          }
          if (!is.list(vinit)) 
              stop("Invalid value for `vinit` parameter")
          if (length(vinit) > nrandom) 
              stop(paste("Vinit must be a list of length", nrandom))
          if (!all(sapply(vinit, function(x) (is.null(x) || is.numeric(x))))) 
              stop("Vinit must contain numeric values")
          if (length(vinit) < nrandom) 
              vinit <- c(vinit, vector("list", nrandom - length(vinit)))
          tname <- names(vinit)
          if (!is.null(tname)) {
              stop("Named initial values not yet supported")
          }
      }
      if (missing(vfixed) || is.null(vfixed)) 
          vfixed <- vector("list", nrandom)
      else {
          if (nrandom == 1) {
              if (is.numeric(vfixed)) 
                  vfixed <- list(vfixed)
              else if (is.list(vfixed)) 
                  vfixed <- list(unlist(vfixed))
          }
          if (!is.list(vfixed)) 
              stop("Invalid value for `vfixed` parameter")
          if (length(vfixed) > nrandom) 
              stop(paste("Vfixed must be a list of length", nrandom))
          if (!all(sapply(vfixed, function(x) (is.null(x) || is.numeric(x))))) 
              stop("Vfixed must contain numeric values")
          if (length(vfixed) < nrandom) 
              vfixed <- c(vfixed, vector("list", nrandom - length(vfixed)))
          tname <- names(vfixed)
          if (!is.null(tname)) {
              temp <- pmatch(tname, names(flist$random), nomatch = 0)
              temp <- c(temp, (1:nrandom)[-temp])
              vfixed <- vfixed[temp]
          }
      }
      newzmat <- function(xmat, xmap) {
          n <- nrow(xmap)
          newz <- matrix(0, nrow = n, ncol = max(xmap))
          for (i in 1:ncol(xmap)) newz[cbind(1:n, xmap[, i])] <- xmat[, 
              i]
          newz
      }
      fmat <- zmat <- matrix(0, nrow = n, ncol = 0)
      ntheta <- integer(nrandom)
      ncoef <- matrix(0L, nrandom, 2, dimnames = list(NULL, c("intercept", 
          "slope")))
      itheta <- NULL
      for (i in 1:nrandom) {
          f2 <- formula2(flist$random[[i]])
          if (f2$intercept & f2$group == 1) 
              stop(paste("Error in random term ", i, ": Random intercepts require a grouping variable", 
                  sep = ""))
          vfun <- varlist[[i]]
          if (!is.null(f2$interaction)) 
              stop("Interactions not yet written")
          cmat <- getcmat(f2$fixed, m)
          groups <- getgroups(f2$group, m)
          ifun <- vfun$initialize(vinit[[i]], vfixed[[i]], intercept = f2$intercept, 
              groups, cmat, control)
          if (!is.null(ifun$error)) 
              stop(paste("In random term ", i, ": ", ifun$error, 
                  sep = ""))
          vparm[[i]] <- ifun$parms
          if (!is.null(ifun$is.variance)) 
              is.variance[i] <- ifun$is.variance
          itheta <- c(itheta, ifun$theta)
          ntheta[i] <- length(ifun$theta)
          if (f2$intercept) {
              if (!is.matrix(ifun$imap) || nrow(ifun$imap) != n) 
                  stop(paste("In random term ", i, ": Invalid intercept matrix F", 
                    sep = ""))
              temp <- sort(unique(c(ifun$imap)))
              if (any(temp != 1:length(temp))) 
                  stop(paste("In random term ", i, ": intercept matrix has an invalid element", 
                    sep = ""))
              if (ncol(fmat) > 0) 
                  fmat <- cbind(fmat, ifun$imap + max(fmat))
              else fmat <- ifun$imap
              ncoef[i, 1] <- 1 + max(ifun$imap) - min(ifun$imap)
          }
          if (length(cmat) > 0) {
              if (is.null(ifun$xmap) || is.null(ifun$X) || !is.matrix(ifun$xmap) || 
                  !is.matrix(ifun$X) || nrow(ifun$xmap) != n || 
                  nrow(ifun$X) != n || ncol(ifun$xmap) != ncol(ifun$X)) 
                  stop(paste("In random term ", i, "invalid X/xmap pair"))
              if (f2$intercept) 
                  xmap <- ifun$xmap - max(ifun$imap)
              else xmap <- ifun$xmap
              if (any(sort(unique(c(xmap))) != 1:max(xmap))) 
                  stop(paste("In random term ", i, ": xmap matrix has an invalid element", 
                    sep = ""))
              temp <- newzmat(ifun$X, xmap)
              ncoef[i, 2] <- ncol(temp)
              zmat <- cbind(zmat, temp)
          }
      }
      if (any(is.variance) & !all(is.variance)) 
          stop("All variance terms must have the same is.variance setting")
      fit <- coxme.fit(X, Y, strats, offset, init, control, weights = weights, 
          ties = ties, row.names(m), fmat, zmat, varlist, vparm, 
          itheta, ntheta, ncoef, refine.n, is.variance = any(is.variance))
      if (is.character(fit)) {
          fit <- list(fail = fit)
          oldClass(fit) <- "coxme"
          return(fit)
      }
      fcoef <- fit$coefficients
      nvar <- length(fcoef)
      if (length(fcoef) > 0 && any(is.na(fcoef))) {
          vars <- (1:length(fcoef))[is.na(fcoef)]
          msg <- paste("X matrix deemed to be singular; variable", 
              paste(vars, collapse = " "))
          warning(msg)
      }
      if (length(fcoef) > 0) {
          names(fit$coefficients) <- dimnames(X)[[2]]
      }
      rlinear <- rep(0, nrow(Y))
      if (length(fmat)) {
          for (i in 1:ncol(fmat)) {
              rlinear <- rlinear + fit$frail[fmat[, i]]
          }
      }
      if (length(zmat)) {
          indx <- if (length(fmat) > 0) 
              max(fmat)
          else 0
          for (i in 1:ncol(zmat)) rlinear <- rlinear + fit$frail[indx + 
              i] * zmat[, i]
      }
      if (nvar == 0) 
          fit$linear.predictor <- rlinear
      else fit$linear.predictor <- as.vector(rlinear + c(X %*% 
          fit$coef))
      newtheta <- random.coef <- list()
      nrandom <- length(varlist)
      sindex <- rep(1:nrandom, ntheta)
      bindex <- rep(row(ncoef), ncoef)
      for (i in 1:nrandom) {
          temp <- varlist[[i]]$wrapup(fit$theta[sindex == i], fit$frail[bindex == 
              i], vparm[[i]])
          newtheta <- c(newtheta, temp$theta)
          if (!is.list(temp$b)) {
              temp$b <- list(temp$b)
              names(temp$b) <- paste("Random", i, sep = "")
          }
          random.coef <- c(random.coef, temp$b)
      }
      fit$frail <- random.coef
      fit$vcoef <- newtheta
      fit$theta <- NULL
      fit$n <- c(sum(Y[, ncol(Y)]), nrow(Y))
      fit$terms <- Terms
      fit$assign <- attr(X, "assign")
      fit$formulaList <- flist
      na.action <- attr(m, "na.action")
      if (length(na.action)) 
          fit$na.action <- na.action
      if (x) {
          fit$x <- X
          if (length(strats)) 
              fit$strata <- strata.keep
      }
      if (y) 
          fit$y <- Y
      if (!is.null(weights) && any(weights != 1)) 
          fit$weights <- weights
      fit$formula <- as.vector(attr(Terms, "formula"))
      fit$call <- Call
      fit$ties <- ties
      names(fit$loglik) <- c("NULL", "Integrated", "Penalized")
      oldClass(fit) <- "coxme"
      fit
  }
  <bytecode: 0xad012f8>
  <environment: namespace:coxme>
   --- function search by body ---
  Function coxme in namespace coxme has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Error in if (class(varlist) == "coxmevar") varlist <- list(varlist) else if (ismat(varlist)) varlist <- list(coxmeMlist(list(varlist))) else { : 
    the condition has length > 1
  Calls: coxme
  Execution halted