Running examples in ‘sanon-Ex.R’ failed
The error most likely occurred in:

> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: coef.sanon
> ### Title: Extract Model Coefficients
> ### Aliases: coef.sanon
> 
> ### ** Examples
> 
> ##### Example 3.1 Randomized Clinical Trial of Chronic Pain #####
> data(cpain)
> out1 = sanon(response ~ grp(treat, ref="placebo") + strt(center) + strt(diagnosis), data=cpain)
 ----------- FAILURE REPORT -------------- 
 --- failure: the condition has length > 1 ---
 --- srcref --- 
: 
 --- package (from environment) --- 
sanon
 --- call from context --- 
sanon.default(outcome = outcome, group = group, strt = strt, 
    covar = covar, catecovar = catecovar, ref = ref, covref = covref, 
    ...)
 --- call from argument --- 
if (class(Y) != "matrix") Y = as.matrix(Y)
 --- R stacktrace ---
where 1: sanon.default(outcome = outcome, group = group, strt = strt, 
    covar = covar, catecovar = catecovar, ref = ref, covref = covref, 
    ...)
where 2: sanon.formula(response ~ grp(treat, ref = "placebo") + strt(center) + 
    strt(diagnosis), data = cpain)
where 3: sanon(response ~ grp(treat, ref = "placebo") + strt(center) + 
    strt(diagnosis), data = cpain)

 --- value of length: 2 type: logical ---
[1] FALSE  TRUE
 --- function from context --- 
function (outcome, group, strt = NULL, covar = NULL, catecovar = NULL, 
    ref = NULL, covref = NULL, P = NULL, res.na.action = "default", 
    ...) 
{
    if (missing(outcome)) 
        stop("outcome should be specified")
    if (missing(group)) 
        stop("group should be specified")
    if (class(P) != "matrix" & !is.null(P)) 
        stop("P should be the matrix class")
    if (!(res.na.action %in% c("default", "LOCF1", "LOCF2", "replace", 
        "remove"))) 
        stop(paste(res.na.action, "not a option for res.na.action"))
    Y = as.data.frame(outcome)
    outnames = colnames(Y)
    if (res.na.action == "remove") 
        Y = na.omit(Y)
    naY = na.action(Y)
    if (!is.null(ncol(Y))) {
        reslevels = lapply(1:ncol(Y), function(i) levels(as.factor(Y[, 
            i])))
        for (i in 1:ncol(Y)) {
            if (!(class(Y[, i]) %in% c("numeric", "integer"))) 
                Y[, i] = as.numeric(Y[, i])
        }
    }
    else {
        Y = as.matrix(Y, 1)
        reslevels = levels(as.factor(Y))
        if (class(Y) != "numeric") 
            Y = as.numeric(Y)
    }
    names(reslevels) = outnames
    if (res.na.action == "LOCF2") {
        if (ncol(Y) == 1) 
            stop("LOCF2 can not be applied to a single response")
        Y = t(apply(Y, 1, function(x) {
            for (i in 2:length(x)) {
                x[i] = ifelse(is.na(x[i]), x[i - 1], x[i])
            }
            x
        }))
    }
    N = ifelse(!is.null(ncol(Y)), nrow(Y), length(Y))
    Nna = length(naY)
    r = ncol(Y)
    if (is.null(strt)) {
        strtnames = NULL
        strt = rep(1, N + Nna)
    }
    else {
        strtnames = paste(colnames(strt), collapse = "*")
    }
    S = as.data.frame(strt)
    if (!is.null(na.action(na.omit(S)))) 
        stop("strata should not have missing values (NA)")
    if (!is.null(ncol(S))) 
        S = apply(S, 1, function(x) paste(x, collapse = "*"))
    if (!is.null(naY)) 
        S = S[-naY]
    strtout = as.factor(S)
    strtlevels = levels(strtout)
    S = as.numeric(as.factor(S))
    if (class(group) == "data.frame") {
        t = group[, 1]
    }
    else {
        t = group
    }
    if (!is.null(na.action(na.omit(t)))) 
        stop("group should not have missing values (NA)")
    if (!is.null(ncol(t))) {
        if (ncol(t) > 1) 
            stop("duplicated group variable")
    }
    if (length(unique(t)) != 2) 
        stop("group should have two categories")
    if (is.null(ref)) {
        grouplevels = levels(factor(t))
    }
    else {
        grouplevels = levels(relevel(factor(t), ref = ref))
    }
    t = ifelse(t == grouplevels[1], -1, 1)
    if (class(t) != "numeric") 
        t = as.numeric(t)
    if (!is.null(naY)) 
        t = t[-naY]
    X = NULL
    if (!is.null(covar)) {
        X = as.data.frame(covar)
    }
    if (!is.null(catecovar)) {
        catecovar = as.data.frame(catecovar)
        catecovar = do.call(cbind, lapply(1:ncol(catecovar), 
            function(x) {
                tmp = catecovar[, x]
                tmpcate = unique(tmp)
                tmpref = tmpcate[length(tmpcate)]
                if (!is.null(covref[x])) {
                  if (!any(tmpcate %in% covref[x])) 
                    stop(paste(covref[x], "not in categories"))
                  tmpref = covref[x]
                }
                tmpcate2 = tmpcate[!(tmpcate %in% tmpref)]
                out = sapply(tmpcate2, function(y) as.numeric(tmp == 
                  y))
                colnames(out) = paste(colnames(catecovar)[x], 
                  "[", tmpcate2, "/", tmpref, "]", sep = "")
                out
            }))
        X = cbind(X, catecovar)
    }
    if (!is.null(X)) {
        if (!is.null(na.action(na.omit(X)))) 
            stop("covariables should not have missing values (NA)")
        if (!is.null(ncol(X))) {
            for (i in 1:ncol(X)) {
                if (!(class(X[, i]) %in% c("numeric", "integer"))) 
                  X[, i] = as.numeric(X[, i])
            }
        }
        else {
            X = as.matrix(X, 1)
            if (class(X) != "numeric") 
                X = as.numeric(X)
        }
        covarnames = colnames(X)
        if (!is.null(naY)) 
            X = X[-naY, ]
        if (class(X) != "matrix") 
            X = as.matrix(X)
    }
    else {
        covarnames = NULL
    }
    M = ifelse(!is.null(X), ncol(X), 0)
    if (res.na.action == "LOCF1") {
        if (ncol(Y) == 1) 
            stop("LOCF1 can not be applied to a single response")
        Z = matrix(1, nrow(Y), ncol(Y))
        Z[, 1] = as.numeric(!is.na(Y[, 1]))
        Y[, 1] = ifelse(is.na(Y[, 1]), 0, Y[, 1])
    }
    else {
        Z = apply(Y, 2, function(x) as.numeric(!is.na(x)))
        Y = apply(Y, 2, function(x) ifelse(is.na(x), 0, x))
    }
    if (class(Y) != "matrix") 
        Y = as.matrix(Y)
    n0 = lapply(1:ncol(Z), function(x) table(factor(t)[Z[, x] == 
        1], factor(S)[Z[, x] == 1]))
    n = sapply(1:r, function(i) sapply(1:N, function(x) n0[[i]][as.character(t[x]), 
        as.character(S[x])]))
    n20 = sapply(1:N, function(x) table(t, S)[as.character(t[x]), 
        as.character(S[x])])
    n2 = n20 %o% rep(1, M)
    n2r = n20 %o% rep(1, r)
    names(n0) = outnames
    for (i in 1:length(n0)) dimnames(n0[[i]]) = list(grouplevels, 
        strtlevels)
    n1 = lapply(n0, function(x) rowSums(x))
    njudge1 = do.call(cbind, n1) < 50
    njudge2 = lapply(n0, function(x) x < 4)
    if (any(njudge1)) {
        njudge1wh = which(njudge1, arr.ind = TRUE)
        njudge1names = apply(njudge1wh, 1, function(x) c(rownames(njudge1)[x[1]], 
            colnames(njudge1)[x[2]]))
    }
    if (any(unlist(njudge2))) {
        njudge2wh = lapply(njudge2, function(x) which(x, arr.ind = TRUE))
        njudge2names = lapply(njudge2wh, function(y) apply(y, 
            1, function(x) c(rownames(njudge2[[1]])[x[1]], colnames(njudge2[[1]])[x[2]])))
        njudge2names = njudge2names[unlist(lapply(njudge2names, 
            class)) == "matrix"]
        njudge2names2 = lapply(njudge2names, function(y) paste(apply(y, 
            2, function(x) paste(x[1], "group in strata", x[2])), 
            collapse = ","))
    }
    U = function(j1) {
        j2 = (1:N)[-j1]
        tmpS = ifelse(S[j1] - S[j2] == 0, 1, 0)
        tmpY = ifelse((Y[rep(j1, N - 1), ] - Y[j2, ]) == 0, 1, 
            0)
        tmptY = ifelse((t[j1] - t[j2]) * (Y[rep(j1, N - 1), ] - 
            Y[j2, ]) * Z[rep(j1, N - 1), ] * Z[j2, ] > 0, 1, 
            0)
        tmpt = ifelse((t[j1] - t[j2])^2 * Z[rep(j1, N - 1), ] * 
            Z[j2, ] > 0, 1, 0)
        if (res.na.action %in% c("replace", "LOCF1", "LOCF2")) 
            tmpt2 = ifelse((t[j1] - t[j2])^2 * (1 - Z[rep(j1, 
                N - 1), ] * Z[j2, ]) > 0, 1, 0)
        if (res.na.action %in% c("replace", "LOCF1", "LOCF2")) {
            tmpU1 = (tmpS * (tmptY + 0.5 * tmpt * tmpY + 0.5 * 
                tmpt2))/(n2r[rep(j1, N - 1), ] + n2r[j2, ] + 
                1)
            tmpU2 = (tmpS * (tmpt + tmpt2))/(n2r[rep(j1, N - 
                1), ] + n2r[j2, ] + 1)
        }
        else {
            tmpU1 = (tmpS * (tmptY + 0.5 * tmpt * tmpY))/(n[rep(j1, 
                N - 1), ] + n[j2, ] + 1)
            tmpU2 = (tmpS * tmpt)/(n[rep(j1, N - 1), ] + n[j2, 
                ] + 1)
        }
        dim(tmpU1) = dim(tmpU2) = c(N - 1, r)
        if (res.na.action == "LOCF1") {
            misidx1 = apply(tmpU1, 1, function(x) any(is.na(x)))
            misidx2 = apply(tmpU2, 1, function(x) any(is.na(x)))
            tmpU1[misidx1, ] = t(apply(tmpU1[misidx1, ], 1, function(x) {
                for (i in 2:length(x)) {
                  x[i] = ifelse(is.na(x[i]), x[i - 1], x[i])
                }
                x
            }))
            tmpU2[misidx2, ] = t(apply(tmpU2[misidx2, ], 1, function(x) {
                for (i in 2:length(x)) {
                  x[i] = ifelse(is.na(x[i]), x[i - 1], x[i])
                }
                x
            }))
        }
        c(apply(tmpU1, 2, mean), apply(tmpU2, 2, mean))
    }
    tmpU = do.call(rbind, lapply(1:N, function(x) U(x)))
    U1j = tmpU[, 1:r]
    U2j = tmpU[, r + (1:r)]
    thetas = apply(tmpU, 2, mean)
    theta1 = thetas[1:r]
    theta2 = thetas[r + (1:r)]
    F = cbind(U1j, U2j)
    meanF = apply(F, 2, mean)
    tmpF = sweep(F, 2, meanF)
    VF = 4 * t(tmpF) %*% (tmpF)/(N * (N - 1))
    if (r > 1) {
        Dtheta1 = diag(c(theta1))
        Dtheta2 = diag(c(theta2))
        xi = solve(Dtheta2) %*% theta1
    }
    else {
        xi = theta1/theta2
    }
    if (r > 1) {
        Dxi = diag(c(xi))
        Dthetas = cbind(diag(r), -solve(Dtheta2)^2)
        Vxi = Dthetas %*% VF %*% t(Dthetas)
    }
    else {
        Dthetas = cbind(1/theta1, -1/theta2)
        Vxi = xi * Dthetas %*% VF %*% t(Dthetas) * xi
    }
    if (!is.null(X)) {
        tU = function(j1) {
            j2 = (1:N)[-j1]
            tmpS = ifelse(S[j1] - S[j2] == 0, 1, 0)
            tmpt1 = 0.5 * (t[j1] - t[j2])
            tmpt2 = ifelse(t[j1] - t[j2] != 0, 1, 0)
            tmpX = X[rep(j1, N - 1), ] - X[j2, ]
            tmpU1 = (tmpS * tmpt1 * tmpX)/(n2[rep(j1, N - 1), 
                ] + n2[j2, ])
            tmpU2 = (tmpS * tmpt2)/(n2[rep(j1, N - 1), ] + n2[j2, 
                ])
            dim(tmpU1) = dim(tmpU2) = c(N - 1, M)
            c(apply(tmpU1, 2, mean), apply(tmpU2, 2, mean))
        }
        tmptU = do.call(rbind, lapply(1:N, function(x) tU(x)))
        tU1j = tmptU[, 1:M]
        tU2j = tmptU[, M + (1:M)]
        phis = apply(tmptU, 2, mean)
        phi1 = phis[1:M]
        phi2 = phis[M + (1:M)]
        g = matrix(phi1/phi2, ncol = 1)
    }
    else {
        g = NULL
    }
    if (!is.null(X)) {
        G = cbind(U1j, tU1j, U2j, tU2j)
    }
    else {
        G = cbind(U1j, U2j)
    }
    meanG = apply(G, 2, mean)
    tmpG = sweep(G, 2, meanG)
    VG = 4 * t(tmpG) %*% (tmpG)/(N * (N - 1))
    if (!is.null(X)) {
        f = rbind(xi, g)
        tG = c(theta1, phi1, theta2, phi2)
    }
    else {
        f = xi
        tG = c(theta1, theta2)
    }
    if (length(f) > 1) {
        if (r > 1) {
            tmptheta = diag(theta2^(-2)) %*% diag(theta1)
            invDtheta2 = solve(Dtheta2)
        }
        else {
            tmptheta = theta1/theta2^2
            invDtheta2 = 1/theta2
        }
        if (M > 1) {
            tmpphi = diag(phi2^(-2)) %*% diag(phi1)
            invphi2 = diag(phi2^(-1))
        }
        else if (M > 0) {
            tmpphi = phi1/phi2^2
            invphi2 = 1/phi2
        }
        else {
            tmpphi = invphi2 = matrix(0, M, M)
        }
        H = rbind(cbind(invDtheta2, matrix(0, r, M), -tmptheta, 
            matrix(0, r, M)), cbind(matrix(0, M, r), invphi2, 
            matrix(0, M, r), -tmpphi))
    }
    else {
        H0 = cbind(1, -1)
        H = f %*% H0 %*% diag(1/tG)
    }
    Vf = H %*% VG %*% t(H)
    if (!is.null(X)) {
        if (is.null(P)) 
            P = rbind(diag(r), matrix(0, M, r))
        f = rbind(xi - 0.5, g)
    }
    else {
        if (is.null(P)) 
            P = diag(r)
        f = xi - 0.5
    }
    allnames = c(outnames, covarnames)
    if (nrow(P) != length(allnames)) {
        stop("The number of row of P should be r+M")
    }
    else {
        rownames(P) = allnames
    }
    advarnames = allnames[apply(P, 1, function(x) all(x == 0))]
    bnames = apply(P, 2, function(x) paste(rownames(P)[which(x == 
        1)], collapse = " + "))
    invVf = try(solve(Vf), silent = TRUE)
    if (class(invVf) == "try-error") {
        warning("Vf is computationally singular.")
        e = 1e-07
        while (class(invVf) == "try-error") {
            invVf = try(solve(Vf + e * diag(ncol(Vf))), silent = TRUE)
            e = 2 * e
        }
    }
    b = solve(t(P) %*% invVf %*% P) %*% (t(P) %*% invVf %*% f)
    Vb = solve(t(P) %*% invVf %*% P)
    if (r > 1) {
        se = sqrt(diag(Vb))
    }
    else {
        se = sqrt(Vb)
    }
    Q = (b/se)^2
    pQ = 1 - pchisq(Q, 1)
    out = list(N = N, Nna = Nna, nhik = n0, nik = n1, xi = xi, 
        g = g, f = f, Vf = Vf, b = b, Vb = Vb, se = se, Q = Q, 
        p = pQ, call = match.call(), outnames = outnames, covarnames = covarnames, 
        advarnames = advarnames, bnames = bnames, reslevels = reslevels, 
        grouplevels = grouplevels, strtout = strtout, strtlevels = strtlevels, 
        strtnames = strtnames, matP = P)
    class(out) = "sanon"
    out
}
<bytecode: 0x564aa477aec0>
<environment: namespace:sanon>
 --- function search by body ---
Function sanon.default in namespace sanon has this body.
 ----------- END OF FAILURE REPORT -------------- 
Fatal error: the condition has length > 1