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

> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: ExampleBSseq
> ### Title: BSseq example dataset
> ### Aliases: ExampleBSseq
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(ExampleBSseq)
> attach(ExampleBSseq)
> model_DNA=pqlseq(RawCountDataSet=mcount, Phenotypes=predictor, 
+   RelatednessMatrix=relatednessmatrix, LibSize=totalcount,
+   fit.model="BMM",numCore=1)
## number of total individuals:  100 
## number of total genes/sites:  5 
## number of adjusted covariates:  0 
# fitting binomial mixed model ... 
 ----------- FAILURE REPORT -------------- 
 --- failure: the condition has length > 1 ---
 --- srcref --- 
: 
 --- package (from environment) --- 
PQLseq
 --- call from context --- 
PQLseq.AI(model0, RelatednessMatrix, maxiter = maxiter, tol = tol, 
    verbose = verbose)
 --- call from argument --- 
if (class(P) == "try-error") {
    stop("Error in P matrix calculation!")
}
 --- R stacktrace ---
where 1: PQLseq.AI(model0, RelatednessMatrix, maxiter = maxiter, tol = tol, 
    verbose = verbose)
where 2: PQLseq.fit(model0, tmpRelatednessMatrix)
where 3: doTryCatch(return(expr), name, parentenv, handler)
where 4: tryCatchOne(expr, names, parentenv, handlers[[1L]])
where 5: tryCatchList(expr, classes, parentenv, handlers)
where 6: tryCatch(expr, error = function(e) {
    call <- conditionCall(e)
    if (!is.null(call)) {
        if (identical(call[[1L]], quote(doTryCatch))) 
            call <- sys.call(-4L)
        dcall <- deparse(call)[1L]
        prefix <- paste("Error in", dcall, ": ")
        LONG <- 75L
        sm <- strsplit(conditionMessage(e), "\n")[[1L]]
        w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w")
        if (is.na(w)) 
            w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L], 
                type = "b")
        if (w > LONG) 
            prefix <- paste0(prefix, "\n  ")
    }
    else prefix <- "Error : "
    msg <- paste0(prefix, conditionMessage(e), "\n")
    .Internal(seterrmessage(msg[1L]))
    if (!silent && isTRUE(getOption("show.error.messages"))) {
        cat(msg, file = outFile)
        .Internal(printDeferredWarnings())
    }
    invisible(structure(msg, class = "try-error", condition = e))
})
where 7: try(PQLseq.fit(model0, tmpRelatednessMatrix))
where 8: eval(c.expr, envir = args, enclos = envir)
where 9: eval(c.expr, envir = args, enclos = envir)
where 10: doTryCatch(return(expr), name, parentenv, handler)
where 11: tryCatchOne(expr, names, parentenv, handlers[[1L]])
where 12: tryCatchList(expr, classes, parentenv, handlers)
where 13: tryCatch(eval(c.expr, envir = args, enclos = envir), error = function(e) e)
where 14: FUN(X[[i]], ...)
where 15: lapply(X = X, FUN = FUN, ...)
where 16: mclapply(argsList, FUN, mc.preschedule = preschedule, mc.set.seed = set.seed, 
    mc.silent = silent, mc.cores = cores)
where 17: e$fun(obj, substitute(ex), parent.frame(), e$data)
where 18: foreach(iVar = 1:numVar, .combine = rbind) %dopar% {
    numAnalysis <- beta <- tau1 <- tau2 <- se_beta <- pvalue <- converged <- h2 <- sigma2 <- overdisp <- NA
    if (verbose) {
        cat(paste("NO. Gene/Site = ", iVar, "\n"))
    }
    if (sum(dim(LibSize) == dim(CountData)) != 2) {
        stop("PQLseq::BMM::ERROR: the dimensions of read counts and total read counts do not match!")
    }
    LibSize <- as.matrix(LibSize)
    if (numCov == 0) {
        model0 <- glm(formula = CountData[iVar, ]/LibSize[iVar, 
            ] ~ Phenotypes, family = binomial(link = "logit"), 
            weights = LibSize[iVar, ])
        idx <- match(rownames(model.frame(formula = CountData[iVar, 
            ]/LibSize[iVar, ] ~ Phenotypes, na.action = na.omit)), 
            rownames(model.frame(formula = CountData[iVar, ]/LibSize[iVar, 
                ] ~ Phenotypes, na.action = na.pass)))
    }
    else {
        model0 <- glm(formula = CountData[iVar, ]/LibSize[iVar, 
            ] ~ Covariates + Phenotypes, family = binomial(link = "logit"), 
            weights = LibSize[iVar, ])
        idx <- match(rownames(model.frame(formula = CountData[iVar, 
            ]/LibSize[iVar, ] ~ Covariates + Phenotypes, na.action = na.omit)), 
            rownames(model.frame(formula = CountData[iVar, ]/LibSize[iVar, 
                ] ~ Covariates + Phenotypes, na.action = na.pass)))
    }
    model0$numTotal <- LibSize[iVar, idx]
    model0$numSucc <- CountData[iVar, idx]
    redflag <- FALSE
    for (ierr in c(2:dim(model.matrix(model0))[2])) {
        if (length(unique(model.matrix(model0)[, ierr])) == 1) {
            warning(paste("PQLseq::BMM::the ", ierr - 1, "-th column of covariates are the same for gene/site ", 
                rownames(CountData)[iVar], "!", sep = ""))
            redflag <- TRUE
        }
    }
    if (!redflag) {
        tmpRelatednessMatrix <- RelatednessMatrix
        if (class(tmpRelatednessMatrix) == "matrix") {
            tmpRelatednessMatrix <- tmpRelatednessMatrix[idx, 
                idx]
        }
        else {
            for (ik in seq_len(length(tmpRelatednessMatrix))) {
                tmpRelatednessMatrix[[ik]] <- tmpRelatednessMatrix[[ik]][idx, 
                  idx]
            }
        }
        names(tmpRelatednessMatrix) <- paste("kins", 1:length(tmpRelatednessMatrix), 
            sep = "")
        model1 <- try(PQLseq.fit(model0, tmpRelatednessMatrix))
        if (class(model1) != "try-error" & !is.null(model1)) {
            if (verbose) {
                cat(paste("PQLseq::BMM::tau = ", model1$theta, 
                  "\n"))
            }
            numAnalysis <- length(idx)
            beta <- model1$coefficients[length(model1$coefficients)]
            se_beta <- sqrt(diag(model1$cov)[length(model1$coefficients)])
            pvalue <- pchisq((beta/se_beta)^2, 1, lower.tail = F)
            sigma2 <- model1$theta[2] + model1$theta[3]
            h2 <- model1$theta[2]/(sigma2)
            tau1 <- model1$theta[2]
            tau2 <- model1$theta[3]
            converged <- model1$converged
        }
        else {
            converged <- FALSE
        }
        res <- data.frame(numIDV = numAnalysis, beta = beta, 
            se_beta = se_beta, pvalue = pvalue, h2 = h2, sigma2 = sigma2, 
            converged = converged)
    }
}
where 19: pqlseq(RawCountDataSet = mcount, Phenotypes = predictor, RelatednessMatrix = relatednessmatrix, 
    LibSize = totalcount, fit.model = "BMM", numCore = 1)

 --- value of length: 2 type: logical ---
[1] FALSE FALSE
 --- function from context --- 
function (model0, RelatednessMatrix, tau = rep(0, length(RelatednessMatrix) + 
    1), fixtau = rep(0, length(RelatednessMatrix) + 1), maxiter = 500, 
    tol = 1e-05, verbose = FALSE) 
{
    if (model0$family$family %in% c("binomial")) {
        y <- model0$numSucc
    }
    else {
        y <- model0$y
    }
    numIDV <- length(y)
    offset <- model0$offset
    if (is.null(offset)) {
        offset <- rep(0, numIDV)
    }
    family <- model0$family
    eta <- model0$linear.predictors
    mu <- model0$fitted.values
    mu.eta <- family$mu.eta(eta)
    D <- mu.eta/sqrt(model0$family$variance(mu))
    if (family$family %in% c("binomial")) {
        mu.eta <- model0$numTotal * mu.eta
        D <- mu.eta/sqrt(model0$numTotal * model0$family$variance(mu))
        mu <- model0$numTotal * mu
    }
    Y <- eta - offset + (y - mu)/mu.eta
    X <- model.matrix(model0)
    alpha <- model0$coef
    if (family$family %in% c("poisson", "binomial")) {
        tau[1] <- 1
        fixtau[1] <- 1
    }
    numK <- length(RelatednessMatrix)
    idxtau <- which(fixtau == 0)
    numK2 <- sum(fixtau == 0)
    if (numK2 > 0) {
        tau[fixtau == 0] <- rep(min(0.9, var(Y)/(numK + 1)), 
            numK2)
        H <- tau[1] * diag(1/D^2)
        for (ik in 1:numK) {
            H <- H + tau[ik + 1] * RelatednessMatrix[[ik]]
        }
        Hinv <- chol2inv(chol(H))
        HinvX <- crossprod(Hinv, X)
        XHinvX <- crossprod(X, HinvX)
        P <- try(Hinv - tcrossprod(tcrossprod(HinvX, chol2inv(chol(XHinvX))), 
            HinvX))
        if (class(P) == "try-error") {
            stop("Error in P matrix calculation!")
        }
        PY <- crossprod(P, Y)
        tau0 <- tau
        for (ik in 1:numK2) {
            if (ik == 1 && fixtau[1] == 0) 
                tau[1] <- max(0, tau0[1] + tau0[1]^2 * (sum((PY/D)^2) - 
                  sum(diag(P)/D^2))/numIDV)
            else {
                PAPY <- crossprod(P, crossprod(RelatednessMatrix[[idxtau[ik] - 
                  1]], PY))
                tau[idxtau[ik]] <- max(0, tau0[idxtau[ik]] + 
                  tau0[idxtau[ik]]^2 * (crossprod(Y, PAPY) - 
                    sum(P * RelatednessMatrix[[idxtau[ik] - 1]]))/numIDV)
            }
        }
    }
    for (iter in seq_len(maxiter)) {
        alpha0 <- alpha
        tau0 <- tau
        model1 <- AI(Y, X, length(RelatednessMatrix), RelatednessMatrix, 
            D^2, tau, fixtau, tol)
        tau <- as.numeric(model1$tau)
        cov <- as.matrix(model1$cov)
        alpha <- as.numeric(model1$alpha)
        eta <- as.numeric(model1$eta) + offset
        mu <- family$linkinv(eta)
        mu.eta <- family$mu.eta(eta)
        D <- mu.eta/sqrt(family$variance(mu))
        if (family$family %in% c("binomial")) {
            mu.eta <- model0$numTotal * mu.eta
            D <- mu.eta/sqrt(model0$numTotal * family$variance(mu))
            mu <- model0$numTotal * mu
        }
        Y <- eta - offset + (y - mu)/mu.eta
        if (2 * max(abs(alpha - alpha0)/(abs(alpha) + abs(alpha0) + 
            tol), abs(tau - tau0)/(abs(tau) + abs(tau0) + tol)) < 
            tol) {
            break
        }
        if (max(tau) > tol^(-2) | any(is.infinite(D)) | any(is.infinite(mu)) | 
            any(is.infinite(eta))) {
            iter <- maxiter
            break
        }
    }
    converged <- ifelse(iter < maxiter, TRUE, FALSE)
    res <- y - mu
    P <- model1$P
    return(list(theta = tau, coefficients = alpha, linear.predictors = eta, 
        fitted.values = mu, Y = Y, P = P, residuals = res, cov = cov, 
        converged = converged))
}
<bytecode: 0xbb1f090>
<environment: namespace:PQLseq>
 --- function search by body ---
Function PQLseq.AI in namespace PQLseq has this body.
 ----------- END OF FAILURE REPORT -------------- 
Fatal error: the condition has length > 1