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

> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: bGLMM
> ### Title: Fit Generalized Mixed-Effects Models
> ### Aliases: bGLMM
> ### Keywords: models methods
> 
> ### ** Examples
>     
> 
> data("soccer")
> ## linear mixed models
> lm1 <- bGLMM(points ~ transfer.spendings + I(transfer.spendings^2)
+        + ave.unfair.score + transfer.receits + ball.possession
+        + tackles + ave.attend + sold.out, rnd = list(team=~1), data = soccer)
[1] "Iteration 1"
[1] "Iteration 2"
[1] "Iteration 3"
[1] "Iteration 4"
[1] "Iteration 5"
[1] "Iteration 6"
[1] "Iteration 7"
[1] "Iteration 8"
[1] "Iteration 9"
[1] "Iteration 10"
[1] "Iteration 11"
[1] "Iteration 12"
[1] "Iteration 13"
[1] "Iteration 14"
[1] "Iteration 15"
[1] "Iteration 16"
[1] "Iteration 17"
[1] "Iteration 18"
[1] "Iteration 19"
[1] "Iteration 20"
[1] "Iteration 21"
[1] "Iteration 22"
[1] "Iteration 23"
[1] "Iteration 24"
[1] "Iteration 25"
[1] "Iteration 26"
[1] "Iteration 27"
[1] "Iteration 28"
[1] "Iteration 29"
[1] "Iteration 30"
[1] "Iteration 31"
[1] "Iteration 32"
[1] "Iteration 33"
[1] "Iteration 34"
[1] "Iteration 35"
[1] "Iteration 36"
[1] "Iteration 37"
[1] "Iteration 38"
[1] "Iteration 39"
[1] "Iteration 40"
[1] "Iteration 41"
[1] "Iteration 42"
[1] "Iteration 43"
[1] "Iteration 44"
[1] "Iteration 45"
[1] "Iteration 46"
[1] "Iteration 47"
[1] "Iteration 48"
[1] "Iteration 49"
[1] "Iteration 50"
[1] "Iteration 51"
[1] "Iteration 52"
[1] "Iteration 53"
[1] "Iteration 54"
[1] "Iteration 55"
[1] "Iteration 56"
[1] "Iteration 57"
[1] "Iteration 58"
[1] "Iteration 59"
[1] "Iteration 60"
[1] "Iteration 61"
[1] "Iteration 62"
[1] "Iteration 63"
[1] "Iteration 64"
[1] "Iteration 65"
[1] "Iteration 66"
[1] "Iteration 67"
[1] "Iteration 68"
[1] "Iteration 69"
[1] "Iteration 70"
[1] "Iteration 71"
[1] "Iteration 72"
[1] "Iteration 73"
[1] "Iteration 74"
[1] "Iteration 75"
[1] "Iteration 76"
[1] "Iteration 77"
[1] "Iteration 78"
[1] "Iteration 79"
[1] "Iteration 80"
[1] "Iteration 81"
[1] "Iteration 82"
[1] "Iteration 83"
[1] "Iteration 84"
[1] "Iteration 85"
[1] "Iteration 86"
[1] "Iteration 87"
[1] "Iteration 88"
[1] "Iteration 89"
[1] "Iteration 90"
[1] "Iteration 91"
[1] "Iteration 92"
[1] "Iteration 93"
[1] "Iteration 94"
[1] "Iteration 95"
[1] "Iteration 96"
[1] "Iteration 97"
[1] "Iteration 98"
[1] "Iteration 99"
[1] "Iteration 100"
[1] "Iteration 101"
[1] "Iteration 102"
[1] "Iteration 103"
[1] "Iteration 104"
[1] "Iteration 105"
[1] "Iteration 106"
[1] "Iteration 107"
[1] "Iteration 108"
[1] "Iteration 109"
[1] "Iteration 110"
[1] "Iteration 111"
[1] "Iteration 112"
[1] "Iteration 113"
[1] "Iteration 114"
[1] "Iteration 115"
[1] "Iteration 116"
[1] "Iteration 117"
[1] "Iteration 118"
[1] "Iteration 119"
[1] "Iteration 120"
[1] "Iteration 121"
[1] "Iteration 122"
[1] "Iteration 123"
[1] "Iteration 124"
[1] "Iteration 125"
[1] "Iteration 126"
[1] "Iteration 127"
[1] "Iteration 128"
[1] "Iteration 129"
[1] "Iteration 130"
[1] "Iteration 131"
[1] "Iteration 132"
[1] "Iteration 133"
[1] "Iteration 134"
[1] "Iteration 135"
[1] "Iteration 136"
[1] "Iteration 137"
[1] "Iteration 138"
[1] "Iteration 139"
[1] "Iteration 140"
[1] "Iteration 141"
[1] "Iteration 142"
[1] "Iteration 143"
[1] "Iteration 144"
[1] "Iteration 145"
[1] "Iteration 146"
[1] "Iteration 147"
[1] "Iteration 148"
[1] "Iteration 149"
[1] "Iteration 150"
[1] "Iteration 151"
[1] "Iteration 152"
[1] "Iteration 153"
[1] "Iteration 154"
[1] "Iteration 155"
[1] "Iteration 156"
[1] "Iteration 157"
[1] "Iteration 158"
[1] "Iteration 159"
[1] "Iteration 160"
[1] "Iteration 161"
[1] "Iteration 162"
[1] "Iteration 163"
[1] "Iteration 164"
[1] "Iteration 165"
[1] "Iteration 166"
[1] "Iteration 167"
[1] "Iteration 168"
[1] "Iteration 169"
[1] "Iteration 170"
[1] "Iteration 171"
[1] "Iteration 172"
[1] "Iteration 173"
[1] "Iteration 174"
[1] "Iteration 175"
[1] "Iteration 176"
[1] "Iteration 177"
[1] "Iteration 178"
[1] "Iteration 179"
[1] "Iteration 180"
[1] "Iteration 181"
[1] "Iteration 182"
[1] "Iteration 183"
[1] "Iteration 184"
[1] "Iteration 185"
[1] "Iteration 186"
[1] "Iteration 187"
[1] "Iteration 188"
[1] "Iteration 189"
[1] "Iteration 190"
[1] "Iteration 191"
[1] "Iteration 192"
[1] "Iteration 193"
[1] "Iteration 194"
[1] "Iteration 195"
[1] "Iteration 196"
[1] "Iteration 197"
[1] "Iteration 198"
[1] "Iteration 199"
[1] "Iteration 200"
[1] "Iteration 201"
[1] "Iteration 202"
[1] "Iteration 203"
[1] "Iteration 204"
[1] "Iteration 205"
[1] "Iteration 206"
[1] "Iteration 207"
[1] "Iteration 208"
[1] "Iteration 209"
[1] "Iteration 210"
[1] "Iteration 211"
[1] "Iteration 212"
[1] "Iteration 213"
[1] "Iteration 214"
[1] "Iteration 215"
[1] "Iteration 216"
[1] "Iteration 217"
[1] "Iteration 218"
[1] "Iteration 219"
[1] "Iteration 220"
[1] "Iteration 221"
[1] "Iteration 222"
[1] "Iteration 223"
[1] "Iteration 224"
[1] "Iteration 225"
[1] "Iteration 226"
[1] "Iteration 227"
[1] "Iteration 228"
[1] "Iteration 229"
[1] "Iteration 230"
[1] "Iteration 231"
[1] "Iteration 232"
[1] "Iteration 233"
[1] "Iteration 234"
[1] "Iteration 235"
[1] "Iteration 236"
[1] "Iteration 237"
[1] "Iteration 238"
[1] "Iteration 239"
[1] "Iteration 240"
[1] "Iteration 241"
[1] "Iteration 242"
[1] "Iteration 243"
[1] "Iteration 244"
[1] "Iteration 245"
[1] "Iteration 246"
[1] "Iteration 247"
[1] "Iteration 248"
[1] "Iteration 249"
[1] "Iteration 250"
[1] "Iteration 251"
[1] "Iteration 252"
[1] "Iteration 253"
[1] "Iteration 254"
[1] "Iteration 255"
[1] "Iteration 256"
[1] "Iteration 257"
[1] "Iteration 258"
[1] "Iteration 259"
[1] "Iteration 260"
[1] "Iteration 261"
[1] "Iteration 262"
[1] "Iteration 263"
[1] "Iteration 264"
[1] "Iteration 265"
[1] "Iteration 266"
[1] "Iteration 267"
[1] "Iteration 268"
[1] "Iteration 269"
[1] "Iteration 270"
[1] "Iteration 271"
[1] "Iteration 272"
[1] "Iteration 273"
[1] "Iteration 274"
[1] "Iteration 275"
[1] "Iteration 276"
[1] "Iteration 277"
[1] "Iteration 278"
[1] "Iteration 279"
[1] "Iteration 280"
[1] "Iteration 281"
[1] "Iteration 282"
[1] "Iteration 283"
[1] "Iteration 284"
[1] "Iteration 285"
[1] "Iteration 286"
[1] "Iteration 287"
[1] "Iteration 288"
[1] "Iteration 289"
[1] "Iteration 290"
[1] "Iteration 291"
[1] "Iteration 292"
[1] "Iteration 293"
[1] "Iteration 294"
[1] "Iteration 295"
[1] "Iteration 296"
[1] "Iteration 297"
[1] "Iteration 298"
[1] "Iteration 299"
[1] "Iteration 300"
[1] "Iteration 301"
[1] "Iteration 302"
[1] "Iteration 303"
[1] "Iteration 304"
[1] "Iteration 305"
[1] "Iteration 306"
[1] "Iteration 307"
[1] "Iteration 308"
[1] "Iteration 309"
[1] "Iteration 310"
[1] "Iteration 311"
[1] "Iteration 312"
[1] "Iteration 313"
[1] "Iteration 314"
[1] "Iteration 315"
[1] "Iteration 316"
[1] "Iteration 317"
[1] "Iteration 318"
[1] "Iteration 319"
[1] "Iteration 320"
[1] "Iteration 321"
[1] "Iteration 322"
[1] "Iteration 323"
[1] "Iteration 324"
[1] "Iteration 325"
[1] "Iteration 326"
[1] "Iteration 327"
[1] "Iteration 328"
[1] "Iteration 329"
[1] "Iteration 330"
[1] "Iteration 331"
[1] "Iteration 332"
[1] "Iteration 333"
[1] "Iteration 334"
[1] "Iteration 335"
[1] "Iteration 336"
[1] "Iteration 337"
[1] "Iteration 338"
[1] "Iteration 339"
[1] "Iteration 340"
[1] "Iteration 341"
[1] "Iteration 342"
[1] "Iteration 343"
[1] "Iteration 344"
[1] "Iteration 345"
[1] "Iteration 346"
[1] "Iteration 347"
[1] "Iteration 348"
[1] "Iteration 349"
[1] "Iteration 350"
[1] "Iteration 351"
[1] "Iteration 352"
[1] "Iteration 353"
[1] "Iteration 354"
[1] "Iteration 355"
[1] "Iteration 356"
[1] "Iteration 357"
[1] "Iteration 358"
[1] "Iteration 359"
[1] "Iteration 360"
[1] "Iteration 361"
[1] "Iteration 362"
[1] "Iteration 363"
[1] "Iteration 364"
[1] "Iteration 365"
[1] "Iteration 366"
[1] "Iteration 367"
[1] "Iteration 368"
[1] "Iteration 369"
>       
> lm2 <- bGLMM(points~transfer.spendings + I(transfer.spendings^2)
+        + ave.unfair.score + transfer.receits + ball.possession
+        + tackles + ave.attend + sold.out, rnd = list(team=~1 + ave.attend), 
+        data = soccer, control = list(steps=10, lin=c("(Intercept)","ave.attend"), 
+        method="REML", nue=1, sel.method="bic"))
[1] "Iteration 1"
 ----------- FAILURE REPORT -------------- 
 --- failure: the condition has length > 1 ---
 --- srcref --- 
: 
 --- package (from environment) --- 
GMMBoost
 --- call from context --- 
est.bGLMM(fix, rnd, ...)
 --- call from argument --- 
if (all(eigen(Q1)$values > 0) & class(Q_solvetest) != "try-error") break
 --- R stacktrace ---
where 1: est.bGLMM(fix, rnd, ...)
where 2: bGLMM.formula(points ~ transfer.spendings + I(transfer.spendings^2) + 
    ave.unfair.score + transfer.receits + ball.possession + tackles + 
    ave.attend + sold.out, rnd = list(team = ~1 + ave.attend), 
    data = soccer, control = list(steps = 10, lin = c("(Intercept)", 
        "ave.attend"), method = "REML", nue = 1, sel.method = "bic"))
where 3: bGLMM(points ~ transfer.spendings + I(transfer.spendings^2) + 
    ave.unfair.score + transfer.receits + ball.possession + tackles + 
    ave.attend + sold.out, rnd = list(team = ~1 + ave.attend), 
    data = soccer, control = list(steps = 10, lin = c("(Intercept)", 
        "ave.attend"), method = "REML", nue = 1, sel.method = "bic"))

 --- value of length: 2 type: logical ---
[1] TRUE TRUE
 --- function from context --- 
function (fix, rnd, data, family = gaussian(link = "identity"), 
    control = list()) 
{
    very.old.names <- attr(terms(fix), "term.labels")
    y <- model.response(model.frame(fix, data))
    X <- model.matrix(fix, data)
    if (dim(X)[2] == 1 || (dim(X)[2] == 2 && colnames(X)[1] == 
        "(Intercept)")) 
        stop("No terms to select! Use glmer, glmmPQL or glmmML!")
    icept.present <- is.element("(Intercept)", colnames(X))
    old.names <- attr(X, "dimnames")[[2]]
    control <- do.call(bGLMMControl, control)
    if (control$print.iter) 
        print(paste("Iteration ", 1, sep = ""))
    add.icept <- FALSE
    if (icept.present && !is.element("(Intercept)", control$lin)) {
        add.icept <- TRUE
        control$lin <- c("(Intercept)", control$lin)
    }
    if (all(X == 0)) {
        X <- rep(1, N)
        X <- as.matrix(X)
    }
    warn.ind <- FALSE
    if (!all(X[, 1] == 1) && is.element("(Intercept)", control$lin)) {
        fix <- update(fix, ~. + 1)
        warn.ind <- TRUE
        very.old.names <- attr(terms(fix), "term.labels")
        X <- model.matrix(fix, data)
        old.names <- attr(X, "dimnames")[[2]]
    }
    rndformula <- as.character(rnd)
    trmsrnd <- terms(rnd[[1]])
    newrndfrml <- "~ -1"
    newrndfrml <- paste(newrndfrml, if (attr(trmsrnd, "intercept")) 
        names(rnd)[1]
    else "", sep = " + ")
    if (length(attr(trmsrnd, "variables")) > 1) {
        newrndfrml <- paste(newrndfrml, paste(sapply(attr(trmsrnd, 
            "term.labels"), function(lbl) {
            paste(lbl, names(rnd)[1], sep = ":")
        }), collapse = " + "), sep = "+")
    }
    W_start <- model.matrix(formula(newrndfrml), data)
    rnlabels <- terms(formula(newrndfrml))
    random.labels <- attr(rnlabels, "term.labels")
    k <- table(data[, colnames(data) == (names(rnd)[1])])
    n <- length(k)
    s <- dim(W_start)[2]/n
    if (s > 1) {
        W <- W_start[, seq(from = 1, to = 1 + (s - 1) * n, by = n)]
        for (i in 2:n) W <- cbind(W, W_start[, seq(from = i, 
            to = i + (s - 1) * n, by = n)])
    }
    else {
        W <- W_start
    }
    if (sum(substr(control$lin, 1, 9) == "as.factor") > 0) {
        group0 <- substr(control$lin, 1, 9) == "as.factor"
        spl0 <- strsplit(control$lin[group0], "\\)")
        control$lin <- control$lin[!group0]
        for (ur in 1:length(spl0)) control$lin <- c(control$lin, 
            old.names[is.element(substr(old.names, 1, nchar(spl0[[ur]])), 
                spl0[[ur]])])
    }
    lin.out <- !is.element(very.old.names, control$lin[1:length(control$lin)])
    very.old.names <- very.old.names[lin.out]
    group <- substr(very.old.names, 1, 9) == "as.factor"
    if (sum(group) > 0) {
        spl <- strsplit(very.old.names[group], "\\(")
        categ.names <- character()
        for (uz in 1:length(spl)) categ.names <- c(categ.names, 
            spl[[uz]][2])
        spl2 <- strsplit(categ.names, "\\)")
        categ.names2 <- character()
        for (uz in 1:length(spl2)) categ.names2 <- c(categ.names2, 
            spl2[[uz]])
        block <- numeric()
        posi <- 1
        for (ip in 1:length(group)) {
            if (!group[ip]) {
                block <- c(block, 1)
            }
            else {
                block <- c(block, length(levels(as.factor(data[, 
                  categ.names2[posi]]))) - 1)
                posi <- posi + 1
            }
        }
    }
    else {
        block <- rep(1, length(group))
    }
    if (length(control$start) == 0) 
        control$start <- c(rep(0, (length(control$lin) + n * 
            s)))
    if (length(control$q_start) == 0) {
        control$q_start <- rep(0.1, s)
        if (s > 1) 
            control$q_start <- diag(control$q_start)
    }
    beta_null <- control$start[1:length(control$lin)]
    ranef_null <- control$start[(length(control$lin) + 1):(length(control$lin) + 
        n * s)]
    q_start <- control$q_start
    N <- length(y)
    Z_fastalles <- X
    no.sel <- is.element(attr(X, "dimnames")[[2]], control$lin)
    U <- as.matrix(X[, !no.sel])
    colnames(U) <- old.names[!no.sel]
    X <- as.matrix(X[, no.sel])
    colnames(X) <- old.names[no.sel]
    lin <- length(control$lin)
    m <- length(block)
    m_alt <- ncol(U)
    lin0 <- sum(beta_null != 0)
    if (s == 1) {
        Q_start <- diag(q_start^2, s)
        p_start <- c(rep(0, lin0), rep((q_start^2)^(-1), n * 
            s))
        P_start <- diag(p_start)
    }
    else {
        Q_start <- q_start
        P_start <- matrix(0, lin0 + n * s, lin0 + n * s)
        for (jf in 1:n) P_start[(lin0 + (jf - 1) * s + 1):(lin0 + 
            jf * s), (lin0 + (jf - 1) * s + 1):(lin0 + jf * s)] <- chol2inv(chol(q_start))
    }
    if (lin > 1) {
        Eta_start <- X %*% beta_null + W %*% ranef_null
    }
    else {
        Eta_start <- rep(beta_null, N) + W %*% ranef_null
    }
    D_start <- as.vector(family$mu.eta(Eta_start))
    Mu_start <- as.vector(family$linkinv(Eta_start))
    Sigma_start <- as.vector(family$variance(Mu_start))
    W0_inv <- D_start * 1/Sigma_start * D_start
    Z_start <- cbind(X[, beta_null != 0], W)
    M0 <- Z_start %*% chol2inv(chol((t(Z_start) %*% (Z_start * 
        W0_inv)) + P_start)) %*% t(Z_start * W0_inv)
    blocksum <- rep(0, m)
    for (jg in 1:m) blocksum[jg] <- sum(block[1:jg])
    Mu <- as.vector(family$linkinv(Eta_start))
    Z_alles <- cbind(X, U, W)
    Z_r <- list()
    Z_r[[1]] <- cbind(X, as.matrix(U[, 1:block[1]]))
    for (r in 2:m) Z_r[[r]] <- cbind(X, as.matrix(U[, (blocksum[r - 
        1] + 1):blocksum[r]]))
    Z_r_gross <- list()
    for (r in 1:m) Z_r_gross[[r]] <- cbind(Z_r[[r]], W)
    Sigma <- as.vector(family$variance(Mu))
    if (control$overdispersion == T) {
        phi_start <- (sum((y - Mu)^2/Mu))/(N - sum(diag(M0)))
        Sigma <- Sigma * phi_start
    }
    Delta <- matrix(0, control$steps, (lin + m_alt + s * n))
    Standard <- matrix(0, control$steps, (lin + m_alt + s * n))
    Standard_ma <- list()
    for (ijk in 1:m) Standard_ma[[ijk]] <- as.numeric()
    Delta[1, 1:lin] <- beta_null
    Delta[1, (lin + m_alt + 1):(lin + m_alt + s * n)] <- t(ranef_null)
    Delta_start <- Delta[1, ]
    IC <- matrix(0, nrow = control$steps, ncol = m)
    IC_stop <- rep(0, control$steps)
    komp <- rep(0, control$steps)
    E <- diag(N)
    M1 <- list()
    diag_H1 <- rep(0, m)
    FinalHat <- list()
    phi <- rep(0, control$steps)
    Q <- array(0, c(s, s, control$steps + 1))
    Q[, , 1] <- Q_start
    D <- as.vector(family$mu.eta(Eta_start))
    Delta_r <- list()
    for (r in 1:m) Delta_r[[r]] <- rep(0, lin + block[r] + s * 
        n)
    Eta_r <- matrix(0, N, m)
    Mu_r <- matrix(0, N, m)
    grp <- rep(1, m_alt)
    for (r in 2:m_alt) grp[((r - 1) + 1):(r)] <- r
    if (s == 1) {
        P_alles <- c(rep(0, lin + m_alt), rep((Q_start^(-1)), 
            n * s))
        P_alles <- diag(P_alles)
    }
    else {
        P_alles <- matrix(0, lin + m_alt + n * s, lin + m_alt + 
            n * s)
        for (jf in 1:n) P_alles[(lin + m_alt + (jf - 1) * s + 
            1):(lin + m_alt + jf * s), (lin + m_alt + (jf - 1) * 
            s + 1):(lin + m_alt + jf * s)] <- chol2inv(chol(Q_start))
    }
    score_alles <- t(Z_alles) %*% ((y - Mu) * D * 1/Sigma)
    F_alles <- t(Z_alles) %*% (Z_alles * D * 1/Sigma * D) + P_alles
    Inv_r <- list()
    l = 1
    for (r in 1:m) {
        if (r == 1) {
            Inv_r[[1]] <- chol2inv(chol(F_alles[c(1:(lin + block[1]), 
                (lin + m_alt + 1):(lin + m_alt + n * s)), c(1:(lin + 
                block[1]), (lin + m_alt + 1):(lin + m_alt + n * 
                s))]))
            Delta_r[[r]] <- Inv_r[[1]] %*% score_alles[c(1:(lin + 
                block[1]), (lin + m_alt + 1):(lin + m_alt + n * 
                s))]
        }
        else {
            Inv_r[[r]] <- chol2inv(chol(F_alles[c(1:lin, (lin + 
                sum(block[1:(r - 1)]) + 1):(lin + sum(block[1:r])), 
                (lin + m_alt + 1):(lin + m_alt + n * s)), c(1:lin, 
                (lin + sum(block[1:(r - 1)]) + 1):(lin + sum(block[1:r])), 
                (lin + m_alt + 1):(lin + m_alt + n * s))]))
            Delta_r[[r]] <- Inv_r[[r]] %*% score_alles[c(1:lin, 
                (lin + sum(block[1:(r - 1)]) + 1):(lin + sum(block[1:r])), 
                (lin + m_alt + 1):(lin + m_alt + n * s))]
        }
        Eta_r[, r] <- Eta_start + Z_r_gross[[r]] %*% Delta_r[[r]]
        Mu_r[, r] <- as.vector(family$linkinv(Eta_r[, r]))
        M1[[r]] <- (Z_r_gross[[r]] * sqrt(Sigma * D * 1/Sigma * 
            D)) %*% Inv_r[[r]] %*% t(Z_r_gross[[r]] * sqrt(D * 
            1/Sigma * D * 1/Sigma))
        diag_H1[r] <- sum(diag(M0 + M1[[r]] %*% (E - M0)))
        if (family$family == "poisson") 
            IC[1, r] <- sum(y * log(Mu_r[, r]) - Mu_r[, r])
        if (family$family == "binomial") 
            IC[1, r] <- sum(y * log(Mu_r[, r]) + (1 - y) * log(1 - 
                Mu_r[, r]))
        if (family$family == "gaussian") 
            IC[1, r] <- sum(y * Mu_r[, r] - 0.5 * (Mu_r[, r]^2))
        if (control$sel.method == "aic") 
            IC[1, r] <- (-IC[1, r]) + diag_H1[r]
        if (control$sel.method == "bic") 
            IC[1, r] <- (-IC[1, r]) + 0.5 * diag_H1[r] * log(n)
    }
    mi_ic <- match(min(IC[1, ]), IC[1, ])
    if (mi_ic == 1) {
        g <- c(rep(T, lin), is.element(grp, 1:blocksum[1]), rep(T, 
            s * n))
    }
    else {
        g <- c(rep(T, lin), is.element(grp, (blocksum[mi_ic - 
            1] + 1):blocksum[mi_ic]), rep(T, s * n))
    }
    Delta_r[[mi_ic]][(lin + 1):(lin + block[mi_ic])] <- control$nue * 
        Delta_r[[mi_ic]][(lin + 1):(lin + block[mi_ic])]
    Delta[1, g] <- Delta[1, g] + t(Delta_r[[mi_ic]])
    Nue_ma <- c(rep(1, lin), rep(control$nue, block[mi_ic]), 
        rep(1, s * n))
    Uu <- (E - (cbind(Z_r[[mi_ic]], W) * sqrt(Sigma * D * 1/Sigma * 
        D)) %*% (Inv_r[[mi_ic]] * Nue_ma) %*% t(cbind(Z_r[[mi_ic]], 
        W) * sqrt(D * 1/Sigma * D * 1/Sigma))) %*% (E - M0)
    FinalHat[[1]] <- E - Uu
    defre <- sum(diag(FinalHat[[1]]))
    komp[1] <- mi_ic
    if (mi_ic == 1) {
        X1 <- cbind(X, as.matrix(U[, 1:block[1]]))
    }
    else {
        X1 <- cbind(X, as.matrix(U[, (blocksum[mi_ic - 1] + 1):blocksum[mi_ic]]))
    }
    if (s == 1) {
        PP <- rep((Q_start^(-1)), n * s)
        PP <- diag(PP)
    }
    else {
        PP <- matrix(0, n * s, n * s)
        for (jf in 1:n) PP[((jf - 1) * s + 1):(jf * s), ((jf - 
            1) * s + 1):(jf * s)] <- chol2inv(chol(Q_start))
    }
    R1 <- chol2inv(chol((t(X1) %*% (X1 * D * 1/Sigma * D)) - 
        (t(X1) %*% (W * D * 1/Sigma * D) %*% chol2inv(chol(t(W) %*% 
            (W * D * 1/Sigma * D) + PP)) %*% t(t(X1) %*% (W * 
            D * 1/Sigma * D)))))
    R2 <- (t(X1 * D * 1/Sigma) - (t(X1) %*% (W * D * 1/Sigma * 
        D) %*% chol2inv(chol(t(W) %*% (W * D * 1/Sigma * D) + 
        PP))) %*% t(W * D * 1/Sigma))
    Q_mi <- control$nue * R1[(lin + block[mi_ic]):(lin + block[mi_ic]), 
        ] %*% R2 %*% (E - M0)
    Standard_ma[[mi_ic]] <- Q_mi
    if (mi_ic == 1) {
        Standard[1, (lin + 1):(lin + blocksum[1])] <- sqrt(diag(Standard_ma[[1]] %*% 
            t(Standard_ma[[1]] * Sigma)))
    }
    else {
        Standard[1, (lin + blocksum[mi_ic - 1] + 1):(lin + blocksum[mi_ic])] <- sqrt(diag(Standard_ma[[mi_ic]] %*% 
            t(Standard_ma[[mi_ic]] * Sigma)))
    }
    Eta <- Z_alles %*% Delta[1, ]
    Mu <- as.vector(family$linkinv(Eta))
    Sigma <- as.vector(family$variance(Mu))
    D <- as.vector(family$mu.eta(Eta))
    if (control$method == "EM") {
        aktuell_vec <- !is.element(Delta[1, ], 0)
        X_aktuell <- Z_alles[, aktuell_vec]
        m_aktuell <- sum(block[is.element(1:m, komp)])
        if (s == 1) {
            P_gross <- c(rep(0, lin + m_aktuell), rep((Q_start^(-1)), 
                n * s))
            P_gross <- diag(P_gross)
        }
        else {
            P_gross <- matrix(0, lin + m_aktuell + n * s, lin + 
                m_aktuell + n * s)
            for (jf in 1:n) P_gross[(lin + m_aktuell + (jf - 
                1) * s + 1):(lin + m_aktuell + jf * s), (lin + 
                m_aktuell + (jf - 1) * s + 1):(lin + m_aktuell + 
                jf * s)] <- chol2inv(chol(Q_start))
        }
        F_gross <- t(X_aktuell) %*% (X_aktuell * D * 1/Sigma * 
            D) + P_gross
        InvFisher <- chol2inv(chol(F_gross))
        Q1 <- InvFisher[(lin + m_aktuell + 1):(lin + m_aktuell + 
            s), (lin + m_aktuell + 1):(lin + m_aktuell + s)] + 
            Delta[1, (lin + m + 1):(lin + m + s)] %*% t(Delta[1, 
                (lin + m + 1):(lin + m + s)])
        for (i in 2:n) Q1 <- Q1 + InvFisher[(lin + m_aktuell + 
            (i - 1) * s + 1):(lin + m_aktuell + i * s), (lin + 
            m_aktuell + (i - 1) * s + 1):(lin + m_aktuell + i * 
            s)] + Delta[1, (lin + m + (i - 1) * s + 1):(lin + 
            m + i * s)] %*% t(Delta[1, (lin + m + (i - 1) * s + 
            1):(lin + m + i * s)])
        Q1 <- 1/n * Q1
    }
    else {
        Eta_tilde <- Eta + (y - Mu) * 1/D
        Betadach <- Delta[1, 1:(lin + m_alt)]
        aktuell_vec <- !is.element(Delta[1, 1:(lin + m_alt)], 
            0)
        X_aktuell <- Z_fastalles[, aktuell_vec]
        if (s == 1) {
            optim.obj <- nlminb(sqrt(Q_start), likelihood_bobyqa, 
                D = D, Sigma = Sigma, X = Z_fastalles, X_aktuell = X_aktuell, 
                Eta_tilde = Eta_tilde, n = n, Betadach = Betadach, 
                W = W, lower = 1e-14, upper = 20)
            Q1 <- as.matrix(optim.obj$par)^2
        }
        else {
            q_start_vec <- c(diag(q_start), q_start[lower.tri(q_start)])
            up1 <- min(20, 50 * max(q_start_vec))
            upp <- rep(up1, length(q_start_vec))
            low <- c(rep(0, s), rep(-up1, 0.5 * (s^2 - s)))
            kkk_vec <- c(rep(-1, s), rep(0.5, 0.5 * (s^2 - s)))
            optim.obj <- try(bobyqa(q_start_vec, likelihood, 
                D = D, Sigma = Sigma, X = Z_fastalles, X_aktuell = X_aktuell, 
                Eta_tilde = Eta_tilde, n = n, s = s, k = k, Betadach = Betadach, 
                W = W, lower = low, upper = upp))
            Q1 <- matrix(0, s, s)
            Q1[lower.tri(Q1)] <- optim.obj$par[(s + 1):(s * (s + 
                1) * 0.5)]
            Q1 <- Q1 + t(Q1)
            diag(Q1) <- (optim.obj$par[1:s])
            for (ttt in 0:100) {
                Q1[lower.tri(Q1)] <- ((0.5)^ttt) * Q1[lower.tri(Q1)]
                Q1[upper.tri(Q1)] <- ((0.5)^ttt) * Q1[upper.tri(Q1)]
                Q_solvetest <- try(solve(Q1))
                if (all(eigen(Q1)$values > 0) & class(Q_solvetest) != 
                  "try-error") 
                  break
            }
        }
    }
    if (family$family == "poisson") 
        IC_stop[1] <- sum(y * log(Mu) - Mu)
    if (family$family == "binomial") 
        IC_stop[1] <- sum(y * log(Mu) + (1 - y) * log(1 - Mu))
    if (family$family == "gaussian") 
        IC_stop[1] <- sum(y * Mu - 0.5 * (Mu^2))
    if (control$sel.method == "aic") 
        IC_stop[1] <- (-IC_stop[1]) + defre
    if (control$sel.method == "bic") 
        IC_stop[1] <- (-IC_stop[1]) + 0.5 * defre * log(n)
    if (control$overdispersion == T) {
        phi[1] <- (sum((y - Mu)^2/Mu))/(N - defre)
        Sigma <- Sigma * phi[1]
    }
    Q[, , 2] <- Q1
    if (control$steps != 1) {
        for (l in 2:control$steps) {
            if (control$print.iter) 
                print(paste("Iteration ", l, sep = ""))
            if (s == 1) {
                P_alles <- c(rep(0, lin + m_alt), rep((Q1^(-1)), 
                  n * s))
                P_alles <- diag(P_alles)
            }
            else {
                P_alles <- matrix(0, lin + m_alt + n * s, lin + 
                  m_alt + n * s)
                for (jf in 1:n) P_alles[(lin + m_alt + (jf - 
                  1) * s + 1):(lin + m_alt + jf * s), (lin + 
                  m_alt + (jf - 1) * s + 1):(lin + m_alt + jf * 
                  s)] <- chol2inv(chol(Q1))
            }
            score_alles <- t(Z_alles) %*% ((y - Mu) * D * 1/Sigma)
            F_alles <- t(Z_alles) %*% (Z_alles * D * 1/Sigma * 
                D) + P_alles
            for (r in 1:m) {
                if (r == 1) {
                  Inv_r[[1]] <- chol2inv(chol(F_alles[c(1:(lin + 
                    block[1]), (lin + m_alt + 1):(lin + m_alt + 
                    n * s)), c(1:(lin + block[1]), (lin + m_alt + 
                    1):(lin + m_alt + n * s))]))
                  Delta_r[[r]] <- Inv_r[[1]] %*% score_alles[c(1:(lin + 
                    block[1]), (lin + m_alt + 1):(lin + m_alt + 
                    n * s))]
                }
                else {
                  Inv_r[[r]] <- chol2inv(chol(F_alles[c(1:lin, 
                    (lin + sum(block[1:(r - 1)]) + 1):(lin + 
                      sum(block[1:r])), (lin + m_alt + 1):(lin + 
                      m_alt + n * s)), c(1:lin, (lin + sum(block[1:(r - 
                    1)]) + 1):(lin + sum(block[1:r])), (lin + 
                    m_alt + 1):(lin + m_alt + n * s))]))
                  Delta_r[[r]] <- Inv_r[[r]] %*% score_alles[c(1:lin, 
                    (lin + sum(block[1:(r - 1)]) + 1):(lin + 
                      sum(block[1:r])), (lin + m_alt + 1):(lin + 
                      m_alt + n * s))]
                }
                Eta_r[, r] <- Eta + Z_r_gross[[r]] %*% Delta_r[[r]]
                Mu_r[, r] <- as.vector(family$linkinv(Eta_r[, 
                  r]))
                M1[[r]] <- (Z_r_gross[[r]] * sqrt(Sigma * D * 
                  1/Sigma * D)) %*% Inv_r[[r]] %*% t(Z_r_gross[[r]] * 
                  sqrt(D * 1/Sigma * D * 1/Sigma))
                diag_H1[r] <- sum(diag(E - ((E - M1[[r]]) %*% 
                  Uu)))
                if (family$family == "poisson") 
                  IC[l, r] <- sum(y * log(Mu_r[, r]) - Mu_r[, 
                    r])
                if (family$family == "binomial") 
                  IC[l, r] <- sum(y * log(Mu_r[, r]) + (1 - y) * 
                    log(1 - Mu_r[, r]))
                if (family$family == "gaussian") 
                  IC[l, r] <- sum(y * Mu_r[, r] - 0.5 * (Mu_r[, 
                    r]^2))
                if (control$sel.method == "aic") 
                  IC[l, r] <- (-IC[l, r]) + diag_H1[r]
                if (control$sel.method == "bic") 
                  IC[l, r] <- (-IC[l, r]) + 0.5 * diag_H1[r] * 
                    log(n)
            }
            mi_ic <- match(min(IC[l, ]), IC[l, ])
            if (mi_ic == 1) {
                g <- c(rep(T, lin), is.element(grp, 1:blocksum[1]), 
                  rep(T, s * n))
            }
            else {
                g <- c(rep(T, lin), is.element(grp, (blocksum[mi_ic - 
                  1] + 1):blocksum[mi_ic]), rep(T, s * n))
            }
            Delta[l, ] <- Delta[l - 1, ]
            Delta_r[[mi_ic]][(lin + 1):(lin + block[mi_ic])] <- control$nue * 
                Delta_r[[mi_ic]][(lin + 1):(lin + block[mi_ic])]
            Delta[l, g] <- Delta[l, g] + t(Delta_r[[mi_ic]])
            komp[l] <- mi_ic
            Nue_ma <- c(rep(1, lin), rep(control$nue, block[mi_ic]), 
                rep(1, s * n))
            Uu <- (E - ((cbind(Z_r[[mi_ic]], W) * sqrt(Sigma * 
                D * 1/Sigma * D)) %*% (Inv_r[[mi_ic]] * Nue_ma) %*% 
                t(cbind(Z_r[[mi_ic]], W) * sqrt(D * 1/Sigma * 
                  D * 1/Sigma)))) %*% Uu
            FinalHat[[l]] <- E - Uu
            defre <- sum(diag(FinalHat[[l]]))
            if (mi_ic == 1) {
                X1 <- cbind(X, as.matrix(U[, 1:block[1]]))
            }
            else {
                X1 <- cbind(X, as.matrix(U[, (blocksum[mi_ic - 
                  1] + 1):blocksum[mi_ic]]))
            }
            if (s == 1) {
                PP <- rep((Q1^(-1)), n * s)
                PP <- diag(PP)
            }
            else {
                PP <- matrix(0, n * s, n * s)
                for (jf in 1:n) PP[((jf - 1) * s + 1):(jf * s), 
                  ((jf - 1) * s + 1):(jf * s)] <- chol2inv(chol(Q1))
            }
            R1 <- chol2inv(chol((t(X1) %*% (X1 * D * 1/Sigma * 
                D)) - (t(X1) %*% (W * D * 1/Sigma * D) %*% chol2inv(chol(t(W) %*% 
                (W * D * 1/Sigma * D) + PP)) %*% t(t(X1) %*% 
                (W * D * 1/Sigma * D)))))
            R2 <- (t(X1 * D * 1/Sigma) - (t(X1) %*% (W * D * 
                1/Sigma * D) %*% chol2inv(chol(t(W) %*% (W * 
                D * 1/Sigma * D) + PP))) %*% t(W * D * 1/Sigma))
            Q_mi <- control$nue * R1[(lin + block[mi_ic]):(lin + 
                block[mi_ic]), ] %*% R2 %*% (E - FinalHat[[l - 
                1]])
            if (is.matrix(Standard_ma[[mi_ic]])) {
                Standard_ma[[mi_ic]] <- Standard_ma[[mi_ic]] + 
                  Q_mi
            }
            else {
                Standard_ma[[mi_ic]] <- Q_mi
            }
            for (jh in as.numeric(levels(as.factor(komp[1:l])))[as.numeric(levels(as.factor(komp[1:l]))) != 
                0]) {
                if (jh == 1) {
                  Standard[l, (lin + 1):(lin + blocksum[1])] <- sqrt(diag(Standard_ma[[1]] %*% 
                    t(Standard_ma[[1]] * Sigma)))
                }
                else {
                  Standard[l, (lin + blocksum[jh - 1] + 1):(lin + 
                    blocksum[jh])] <- sqrt(diag(Standard_ma[[jh]] %*% 
                    t(Standard_ma[[jh]] * Sigma)))
                }
            }
            Eta <- Z_alles %*% Delta[l, ]
            Mu <- as.vector(family$linkinv(Eta))
            Sigma <- as.vector(family$variance(Mu))
            D <- as.vector(family$mu.eta(Eta))
            if (control$method == "EM") {
                aktuell_vec <- !is.element(Delta[l, ], 0)
                X_aktuell <- Z_alles[, aktuell_vec]
                m_aktuell <- sum(block[is.element(1:m, komp)])
                if (s == 1) {
                  P_gross <- c(rep(0, lin + m_aktuell), rep((Q1^(-1)), 
                    n * s))
                  P_gross <- diag(P_gross)
                }
                else {
                  P_gross <- matrix(0, lin + m_aktuell + n * 
                    s, lin + m_aktuell + n * s)
                  for (jf in 1:n) P_gross[(lin + m_aktuell + 
                    (jf - 1) * s + 1):(lin + m_aktuell + jf * 
                    s), (lin + m_aktuell + (jf - 1) * s + 1):(lin + 
                    m_aktuell + jf * s)] <- chol2inv(chol(Q1))
                }
                F_gross <- t(X_aktuell) %*% (X_aktuell * D * 
                  1/Sigma * D) + P_gross
                InvFisher <- chol2inv(chol(F_gross))
                Q1 <- InvFisher[(lin + m_aktuell + 1):(lin + 
                  m_aktuell + s), (lin + m_aktuell + 1):(lin + 
                  m_aktuell + s)] + Delta[l, (lin + m + 1):(lin + 
                  m + s)] %*% t(Delta[l, (lin + m + 1):(lin + 
                  m + s)])
                for (i in 2:n) Q1 <- Q1 + InvFisher[(lin + m_aktuell + 
                  (i - 1) * s + 1):(lin + m_aktuell + i * s), 
                  (lin + m_aktuell + (i - 1) * s + 1):(lin + 
                    m_aktuell + i * s)] + Delta[l, (lin + m + 
                  (i - 1) * s + 1):(lin + m + i * s)] %*% t(Delta[l, 
                  (lin + m + (i - 1) * s + 1):(lin + m + i * 
                    s)])
                Q1 <- 1/n * Q1
            }
            else {
                Eta_tilde <- Eta + (y - Mu) * 1/D
                Betadach <- Delta[l, 1:(lin + m_alt)]
                aktuell_vec <- !is.element(Delta[l, 1:(lin + 
                  m_alt)], 0)
                X_aktuell <- Z_fastalles[, aktuell_vec]
                if (s == 1) {
                  optim.obj <- try(nlminb(sqrt(Q1), likelihood_bobyqa, 
                    D = D, Sigma = Sigma, X = Z_fastalles, X_aktuell = X_aktuell, 
                    Eta_tilde = Eta_tilde, n = n, Betadach = Betadach, 
                    W = W, lower = 1e-12, upper = 20))
                  Q1 <- as.matrix(optim.obj$par)^2
                }
                else {
                  Q1_vec <- c(diag(Q1), Q1[lower.tri(Q1)])
                  optim.obj <- try(bobyqa(Q1_vec, likelihood, 
                    D = D, Sigma = Sigma, X = Z_fastalles, X_aktuell = X_aktuell, 
                    Eta_tilde = Eta_tilde, n = n, s = s, k = k, 
                    Betadach = Betadach, W = W, lower = low, 
                    upper = upp))
                  Q1 <- matrix(0, s, s)
                  Q1[lower.tri(Q1)] <- optim.obj$par[(s + 1):(s * 
                    (s + 1) * 0.5)]
                  Q1 <- Q1 + t(Q1)
                  diag(Q1) <- (optim.obj$par[1:s])
                  for (ttt in 0:100) {
                    Q1[lower.tri(Q1)] <- ((0.5)^ttt) * Q1[lower.tri(Q1)]
                    Q1[upper.tri(Q1)] <- ((0.5)^ttt) * Q1[upper.tri(Q1)]
                    Q_solvetest <- try(solve(Q1))
                    if (all(eigen(Q1)$values > 0) & class(Q_solvetest) != 
                      "try-error") 
                      break
                  }
                }
            }
            if (family$family == "poisson") 
                IC_stop[l] <- sum(y * log(Mu) - Mu)
            if (family$family == "binomial") 
                IC_stop[l] <- sum(y * log(Mu) + (1 - y) * log(1 - 
                  Mu))
            if (family$family == "gaussian") 
                IC_stop[l] <- sum(y * Mu - 0.5 * (Mu^2))
            if (control$sel.method == "aic") 
                IC_stop[l] <- (-IC_stop[l]) + defre
            if (control$sel.method == "bic") 
                IC_stop[l] <- (-IC_stop[l]) + 0.5 * defre * log(n)
            if (control$overdispersion == T) {
                phi[l] <- (sum((y - Mu)^2/Mu))/(N - defre)
                Sigma <- Sigma * phi[l]
            }
            Q[, , l + 1] <- Q1
            kritval <- sqrt(sum((Delta[l - 1, ] - Delta[l, ])^2))/sqrt(sum(Delta[l, 
                ]^2))
            if (kritval < 1e-08) 
                break
            kritval2 <- abs(IC_stop[l - 1] - IC_stop[l])/abs(IC_stop[l])
            if (kritval2 < 1e-08) 
                break
        }
    }
    if (l < control$steps) {
        for (wj in (l + 1):control$steps) {
            Delta[wj, ] <- Delta[l, ]
            Q[, , wj + 1] <- Q1
            IC[wj, ] <- IC[l, ]
            IC_stop[wj] <- IC_stop[l]
        }
    }
    opt <- match(min(IC_stop), IC_stop)
    if (control$OPT == TRUE) {
        Delta_neu <- Delta[opt, ]
    }
    else {
        Delta_neu <- Delta[l, ]
    }
    Eta_opt <- Z_alles %*% Delta_neu
    Mu_opt <- as.vector(family$linkinv(Eta_opt))
    Sigma_opt <- as.vector(family$variance(Mu_opt))
    D_opt <- as.vector(family$mu.eta(Eta_opt))
    aktuell_vec <- !is.element(Delta_neu, 0)
    X_aktuell <- Z_alles[, aktuell_vec]
    m_aktuell <- sum(block[is.element(1:m, komp[1:opt])])
    if (s == 1) {
        P_opt <- c(rep(0, lin + m_aktuell), rep((Q[, , opt + 
            1]^(-1)), n * s))
        P_opt <- diag(P_opt)
    }
    else {
        P_opt <- matrix(0, lin + m_aktuell + n * s, lin + m_aktuell + 
            n * s)
        for (jf in 1:n) P_opt[(lin + m_aktuell + (jf - 1) * s + 
            1):(lin + m_aktuell + jf * s), (lin + m_aktuell + 
            (jf - 1) * s + 1):(lin + m_aktuell + jf * s)] <- chol2inv(chol(Q[, 
            , opt + 1]))
    }
    F_opt <- t(X_aktuell) %*% (X_aktuell * D_opt * 1/Sigma_opt * 
        D_opt) + P_opt
    Inv_F_opt <- chol2inv(chol(F_opt))
    Lin_errors <- sqrt(diag(Inv_F_opt))
    Standard_errors <- Standard[opt, ]
    g3 <- c(rep(T, lin), rep(F, m_alt), rep(T, s * n))
    g4 <- c(rep(T, lin), rep(F, m_aktuell), rep(T, s * n))
    Standard_errors[g3] <- Lin_errors[g4]
    y_dach <- as.vector(family$linkinv(Eta_opt))
    if (s == 1) {
        Qfinal <- sqrt(Q[, , opt + 1])
    }
    else {
        Qfinal <- Q[, , opt + 1]
    }
    if (!is.matrix(Qfinal)) 
        Qfinal <- as.matrix(Qfinal)
    colnames(Qfinal) <- random.labels
    rownames(Qfinal) <- random.labels
    names(Delta_neu)[1:lin] <- colnames(X)[1:lin]
    names(Standard_errors)[1:lin] <- colnames(X)[1:lin]
    names(Delta_neu)[(lin + 1):(lin + m_alt)] <- colnames(U)
    names(Standard_errors)[(lin + 1):(lin + m_alt)] <- colnames(U)
    permut <- numeric()
    for (ys in 1:length(old.names)) permut <- c(permut, match(old.names[ys], 
        names(Delta_neu)))
    Delta_neu[1:length(old.names)] <- Delta_neu[permut]
    Standard_errors[1:length(old.names)] <- Standard_errors[permut]
    names(Delta_neu)[1:length(old.names)] <- old.names
    names(Standard_errors)[1:length(old.names)] <- old.names
    Delta[, 1:length(old.names)] <- Delta[, permut]
    names(Delta_neu)[(lin + m_alt + 1):(lin + m_alt + n * s)] <- colnames(W)
    names(Standard_errors)[(lin + m_alt + 1):(lin + m_alt + n * 
        s)] <- colnames(W)
    colnames(Delta) <- c(old.names, colnames(W))
    colnames(IC) <- very.old.names
    comp <- character()
    for (oi in 1:length(komp)) comp <- c(comp, very.old.names[komp[oi]])
    if (warn.ind) {
        if (lin == 1) {
            cat("Warning:\n")
            cat("At least one unpenalized component has to be incorporated - \n      intercept is added!\n")
        }
        else {
            cat("Warning:\n")
            cat("Intercept is added to the formula as it is specified in the control argument!\n")
        }
    }
    if (add.icept) {
        cat("Warning:\n")
        cat("Intercept is added to unpenalized component!\n")
    }
    ret.obj = list()
    ret.obj$IC <- IC
    ret.obj$IC_sel <- IC_stop
    ret.obj$opt <- opt
    ret.obj$Deltamatrix <- Delta
    ret.obj$ranef <- Delta_neu[(lin + m_alt + 1):(lin + m_alt + 
        n * s)]
    ret.obj$coefficients <- Delta_neu[1:(lin + m_alt)]
    ret.obj$fixerror <- Standard_errors[1:(lin + m_alt)]
    ret.obj$ranerror <- Standard_errors[(lin + m_alt + 1):(lin + 
        m_alt + n * s)]
    ret.obj$components <- comp
    ret.obj$Q_long <- Q
    ret.obj$Q <- Qfinal
    ret.obj$y_hat <- Mu_opt
    ret.obj$HatMatrix <- FinalHat[[opt]]
    ret.obj$phi <- phi[opt]
    ret.obj$family <- family
    ret.obj$fix <- fix
    ret.obj$newrndfrml <- newrndfrml
    ret.obj$subject <- names(rnd)[1]
    ret.obj$data <- data
    return(ret.obj)
}
<bytecode: 0x17c4a38>
<environment: namespace:GMMBoost>
 --- function search by body ---
Function est.bGLMM in namespace GMMBoost has this body.
 ----------- END OF FAILURE REPORT -------------- 
Fatal error: the condition has length > 1