R/glmmLasso_RE.R

Defines functions est.glmmLasso.RE

est.glmmLasso.RE<-function(fix,rnd,data,lambda,family,final.re,switch.NR,control)
{  
  control<-do.call(glmmLassoControl, control)
  
  ## Print stuff.
  if(is.null(control$flushit))
    control$flushit <- TRUE
  
  ia <- if(control$flushit) interactive() else FALSE
   
  fix.old<-fix
  
  if(!is.null(fix))
  {  
    if(grepl("\\*", fix[3]))
      stop("Usage of '*' not allowed in formula! Please specify the corresponding variables separately.")  
    
    ic.dummy<-attr(terms(fix),"intercept")  
    y <- model.response(model.frame(fix, data))
    very.old.names<-attr(terms(fix, keep.order = TRUE),"term.labels")
    
    if(ic.dummy!=1 && sum(substr(very.old.names,1,9)=="as.factor")>0){
      fix.help <- update.formula(fix,~ .+1)
      orig.names <- colnames(model.matrix(terms(fix.help, keep.order = TRUE), data))[-1]
    }else{
      orig.names <- colnames(model.matrix(terms(fix, keep.order = TRUE), data))
    }
    
    if(!is.null(control$index))
    {
      order.vec<-order(control$index)
      very.old.names<-very.old.names[order.vec]
      control$index<-control$index[order.vec]
    }else{
      control$index<-1:length(very.old.names)
      order.vec <-NULL
    }
    
    if(length(control$index)!=length(very.old.names))
      stop("Length of vector defining the grouping of the variables doesn't match with 
           the formula!")
    
    attr(control$index,"names")<-very.old.names
    
    fix <- formula(paste("y~-1+",paste(very.old.names,collapse="+"))) 

    if(ic.dummy==1)
    {
      fix<-update.formula(fix,~ .+1) 
      control$index<-c(NA,control$index)
      names(control$index)[1]<-"(Intercept)"
    }
    
    index.new<-c()
    fac.variab<-logical()
    for(i in 1:length(control$index))
    {
      if(!grepl("as.factor",names(control$index)[i]))
      {
        index.new<-c(index.new,control$index[i]) 
        fac.variab<-c(fac.variab,F)
      }else{
        if(!grepl("\\:", names(control$index)[i]))
        {  
          fac.name<-strsplit(strsplit(names(control$index)[i],"\\(")[[1]][2],"\\)")[[1]][1]
        }else{
          fac.name<-unlist(strsplit(unlist(strsplit(unlist(strsplit(names(control$index)[i],"\\(")),"\\)")),"\\:"))
          fac.name<-paste(fac.name[2],":",fac.name[length(fac.name)],sep="")
        }
        if(!grepl("\\:", fac.name))
        {
          length.fac<-length(levels(as.factor(data[,fac.name])))-1
        }else{
          length.fac<-(length(levels(data[,strsplit(fac.name,":")[[1]][1]]))-1)*(length(levels(data[,strsplit(fac.name,":")[[1]][2]]))-1)
        }
        index.new<-c(index.new,rep(control$index[i],length.fac))
        fac.variab<-c(fac.variab,rep(T,length.fac))
      }
    }
    
  #  browser()
    
    
    
    if(ic.dummy!=1 && sum(substr(very.old.names,1,9)=="as.factor")>0){
      fix.help <- update.formula(fix,~ .+1)
      X <- model.matrix(terms(fix.help, keep.order = TRUE), data)[,-1]
      colnames(X) <- colnames(model.matrix(terms(fix.help, keep.order = TRUE), data)[,-1])
    }else{
      X <- model.matrix(terms(fix, keep.order = TRUE), data)  
      colnames(X) <- colnames(model.matrix(terms(fix, keep.order = TRUE), data))
    }
    
    K <- NULL
    
    if(!is.null(family$multivariate)){
      y.fac <- as.factor(y)
      K <- length(levels(y.fac))-1
      if(family$family=="acat"){
        y <- c(t(model.matrix(~0+y.fac, contrasts = list(y.fac = "contr.treatment"))[,-length(levels(y.fac))]))
      }
      if(family$family=="cumulative"){
        get.resp <- function(x){as.numeric(x <= 1:K)}
        y <- c(sapply(y,get.resp))
      }
      yhelp <- as.matrix(model.matrix(~0+y.fac, contrasts = list(y.fac = "contr.treatment")))
    }
    
     if(!is.null(family$multivariate)){
      if(all(X[,1]==1)){
        X <- X[,-1]
      }
      names.x <- colnames(X)
      theta <- matrix(rep(diag(1,K),nrow(X)),ncol=K,byrow=TRUE)
      X <- cbind(theta, matrix(rep(X,each=K),ncol=ncol(X)))
      colnames(X) <- c(paste0("theta",1:K),names.x)
      if(orig.names[1]=="(Intercept)")
      {  
        orig.names <- orig.names[-1]
        index.new <- c(rep(NA,K),index.new[-1]) 
      }else{
        index.new <- c(rep(NA,K),index.new) 
      }  
      orig.names <- c(paste0("theta",1:K),orig.names)
    }

    transf.names <- colnames(X)
    center <- control$center
    standardize <- control$standardize
    ####### Center & Standardization
    
    
    ## Which are the non-penalized parameters?
    any.notpen    <- any(is.na(index.new))
    inotpen.which <- which(is.na(index.new))
    nrnotpen      <- length(inotpen.which)
    
    intercept.which <- which(apply(X == 1, 2, all))
    has.intercept   <- length(intercept.which)
    
    ## Index vector of the penalized parameter groups
    if(any.notpen){
      ipen <- index.new[-inotpen.which]
      ipen.which <- split((1:ncol(X))[-inotpen.which], ipen)
    }else{
      if(has.intercept)
        warning("All groups are penalized, including the intercept.")
      ipen <- index.new
      ipen.which <- split((1:ncol(X)), ipen)
    }
    
    if(center){
      if(!has.intercept & is.null(family$multivariate)) ## could be removed; already handled above
        stop("Need intercept term when using center = TRUE")
      
      mu.x                 <- apply(as.matrix(X[,-intercept.which]), 2, mean)
      X[,-intercept.which] <- sweep(as.matrix(X[,-intercept.which]), 2, mu.x)
    }
    
    ## Standardize the design matrix -> blockwise orthonormalization
    if(standardize){
      ##warning("...Using standardized design matrix.\n")
      stand        <- blockstand(X, ipen.which, inotpen.which)
      X            <- stand$x
      scale.pen    <- stand$scale.pen
      scale.notpen <- stand$scale.notpen
    }
    
    ##############
    
    if(ncol(X)==1)
    {
      if(colnames(X)=="(Intercept)")
        stop("No terms to select! Use glmer, glmmPQL or glmmML!")
    }
    
    #old.names<-attr(X,"dimnames")[[2]]
    
    very.old.names<-very.old.names[!is.na(control$index)]
    ############
    ############
    ############
  }else{
    y <- control$y
    X <- control$X
    K <- control$K
    index.new <- control$index.new
    W.index <- control$W.index
    orig.names <- transf.names <- colnames(X)
    standardize <- center <- FALSE
  }
  
  
  if(control$print.iter)
    #     message()
  {
    cat(if(ia) "\r" else NULL)
    cat("Iteration  1")
    if(.Platform$OS.type != "unix" & ia) flush.console()
  }
  
  random.factor.help <- FALSE
  if(is.list(rnd))
  {
    rnd.len<-length(rnd)
    
    if(rnd.len==1)
    {
      rndformula <- as.character(rnd)
      
      trmsrnd <- terms(rnd[[1]])
      
      
      if(!is.factor(data[,names(rnd)[1]]))
      {
        data[,names(rnd)[1]] <- as.factor(data[,names(rnd)[1]])
        warning("Cluster variable should be specified as a factor variable!")  
      }
      
      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="+") 
      }
      
      data.new <- data
      data.new[,colnames(X)] <- X
      
      W_start <- model.matrix(formula(newrndfrml), data.new)
      
      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
      }
    }else{     
      rndformula <- list()
      newrndfrml <- list()
      n <- numeric()
      s <- numeric()
      k<-NULL
      random.labels<-list()
      W<-NULL
      
      for (zu in 1:rnd.len)
      {
        rndformula[[zu]] <- as.character(rnd[[zu]])
        
        trmsrnd <- terms(rnd[[zu]])
        
        #####
        ran.lab.vec <- attr(trmsrnd,"term.labels") 
        random.factor.name <- character()
        
        if(length(seq_along(ran.lab.vec))>0){
          for(er in seq_along(ran.lab.vec))
          {  
            if(ran.lab.vec[er] %in% colnames(data))
            {
              if(is.factor(data[,ran.lab.vec[er]]))
              {
                random.factor.help <- TRUE
                random.factor.name <- c(random.factor.name,colnames(model.matrix(formula(paste0("~",ran.lab.vec[er])), data[1,]))[-1])
              }else{
                random.factor.name <- c(random.factor.name,ran.lab.vec[er])
              }
            }else if(substr(ran.lab.vec[er],1,9)=="as.factor")
            {
              random.factor.help <- TRUE
              name.help <- strsplit(strsplit(ran.lab.vec[er],"\\(")[[1]][2],"\\)")[[1]][1]
              data[,name.help] <- as.factor(data[,name.help])
              random.factor.name <- c(random.factor.name,colnames(model.matrix(formula(paste0("~",name.help)), data[1,]))[-1])
            }}}
        random.factor.name <- c(names(rnd)[zu],paste(names(rnd)[zu],random.factor.name,sep=":"))
        #####
        if(!is.factor(data[,names(rnd)[zu]]))
        {
          data[,names(rnd)[zu]] <- as.factor(data[,names(rnd)[zu]])
          warning("Cluster variable should be specified as a factor variable!")  
        }
        
        newrndfrml[[zu]] <- "~ -1"
        newrndfrml[[zu]] <- paste(newrndfrml[[zu]],  if(attr(trmsrnd, "intercept")) names(rnd)[zu] else "", sep=" + ")
        
        if(length(attr(trmsrnd, "variables"))>1)
        {
          newrndfrml[[zu]] <- paste(newrndfrml[[zu]],  
                                    paste(sapply(attr(trmsrnd,"term.labels"), function(lbl){
                                      paste(lbl, names(rnd)[zu], sep=":")
                                    }), collapse=" + "), sep="+") }
        W_start <- model.matrix(formula(newrndfrml[[zu]]), data)
        
        rnlabels<-terms(formula(newrndfrml[[zu]]))
        if(random.factor.help) random.labels[[zu]]<-random.factor.name else random.labels[[zu]]<-attr(rnlabels,"term.labels")
        k1<-table(data[,colnames(data)==(names(rnd)[zu])])   
        n[zu]<-length(k1)
        s[zu]<-dim(W_start)[2]/n[zu]
        
        if(s[zu]>1)
        {
          W2<-W_start[,seq(from=1,to=1+(s[zu]-1)*n[zu],by=n[zu])]
          for (i in 2:n[zu])
            W2<-cbind(W2,W_start[,seq(from=i,to=i+(s[zu]-1)*n[zu],by=n[zu])])
        }else{
          W2<-W_start
        }
        W<-cbind(W,W2)
        k<-c(k,k1)
      }
    }
    subject.names<-names(rnd)
    
    if(!is.null(family$multivariate))
    {
      names.of.W <- colnames(W)
      W <- matrix(rep(W,each=K),ncol=ncol(W))
      colnames(W) <- names.of.W
    }
    
  }else{
    W<-rnd
    if(!is.null(control$W.index))
    {
      s<-length(unique(control$W.index))
      n<-ncol(W)/s
    }else{
      n<-ncol(W)
      s<-1
    }
    
    newrndfrml<-NULL  
    random.labels<-attr(W,"W.name")   
    k<-NULL
    rnd.len <- 1
    attr(rnd,"names")<-colnames(rnd)
    subject.names<-colnames(rnd)  
  }

  block<-as.numeric(table(index.new[!is.na(index.new)]))
  
  BLOCK<-FALSE
  if(!all(block==1))
    BLOCK<-TRUE
  
  lin<-ncol(X)
  
  if(is.null(control$start))
  {
    control$start<-c(rep(0,lin+n%*%s))
    if(family$family=="gaussian")
       control$start[1] <- mean(y)
    if(family$family=="poisson")
       control$start[1] <- log(mean(y))
    if(family$family=="binomial")
       control$start[1] <- log(mean(y)/(1-mean(y)))
  }  
  
  if(family$family=="cumulative" & all(control$start==0))
    control$start[1:K] <- ((1:K)-mean(1:K))/K
  
  if(is.null(control$q_start))
  {
    control$q_start<-rep(0.1,sum(s))
    if(sum(s)>1)
      control$q_start<-diag(control$q_start)
  }
  
  q_start<-control$q_start
  
  N<-length(y)
  
  beta_null<-control$start[1:lin]
  # if(is.null(attr(beta_null,"names")))
  #   attr(beta_null,"names")<-orig.names
  # beta_null<-beta_null[colnames(X)]
  if(!is.null(order.vec))
    beta_null<-beta_null[c(1,order.vec+1)]
  attr(beta_null,"names") <- colnames(X)
  
  ranef_null<-control$start[(lin+1):(lin+n%*%s)]
  
  #Z_fastalles<-X
  
  if(!control$overdispersion && family$family=="gaussian")
    control$overdispersion<-T
  
  phi <- 1
  #######################################################################  
  ######################## allow switch to Newton Raphson ###############
  #######################################################################  
  if(switch.NR)
  {
    #######################################################################  
    ###########################  1. No Smooth #############################  
    #######################################################################  
    if(is.null(control$smooth))
    {  
      if(lin>1)
      {
        Eta_start<-X%*%beta_null+W%*%ranef_null
      }else{
        Eta_start<-rep(beta_null,N)+W%*%ranef_null
      }
      
      if(is.null(family$multivariate)){
        Mu<-family$linkinv(Eta_start)
        D<-family$mu.eta(Eta_start)
        SigmaInv <- 1/family$variance(Mu)
      }else{
        D <- family$deriv.mat(Eta_start, K)
        Mu<-family$linkinv(Eta_start, K)
        SigmaInv <- family$SigmaInv(Mu, K)
      }
      
      if(rnd.len==1)
      {
        lin0<-sum(beta_null!=0)
          Q_start<-q_start
      }else{
        lin0<-sum(beta_null!=0)
        if(all(s==1))
        {
          Q_start<-diag(diag(q_start),sum(s))
        }else{
          Q_start<-matrix(0,sum(s),sum(s))
          Q_start[1:s[1],1:s[1]]<-q_start[1:s[1],1:s[1]]
          
          for (zu in 2:rnd.len)
          {
            Q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
          }
        }
      }
      
      U<-X[,!is.na(index.new),drop=FALSE]
      X<-X[,is.na(index.new),drop=FALSE]
      
      final.names<-c(colnames(X),colnames(U))
      
      q<-ncol(X)

      Z_fastalles<-cbind(X,U)
      Z_alles<-cbind(Z_fastalles,W)
      ########################################################## some definitions ################################################
      Delta<-matrix(0,control$steps,(lin+n%*%s))
      Delta[1,1:lin]<-beta_null[final.names]
      Delta[1,(lin+1):(lin+n%*%s)]<-t(ranef_null)
      active_old<-c(rep(T,q),!is.element(Delta[1,(q+1):lin],0),rep(T,n%*%s))
      
      Eta.ma<-matrix(0,control$steps+1,N)
      Eta.ma[1,]<-Eta_start
      
      logLik.vec<-c()
      
      control$epsilon<-control$epsilon*sqrt(dim(Delta)[2])
      
      Delta_start<-Delta[1,]
      
      Q_inv<-NULL
      Q_inv.old.temp<-NULL
      Q_inv.start<-NULL

      Q<-list()
      Q[[1]]<-Q_start

      l=1
      
      if(rnd.len==1)
      {
        if(s==1)
        {
          P1<-c(rep(0,lin),rep(1/Q_start,n*s))
          P1<-diag(P1)
        }else{
          Q_inv.start<-chol2inv(chol(Q_start))
          P1<-matrix(0,lin+n%*%s,lin+n%*%s)
          for(j in 1:n)
            P1[(lin+(j-1)*s+1):(lin+j*s),(lin+(j-1)*s+1):(lin+j*s)]<-Q_inv.start
        }
      }else{
        if(all(s==1))
        {
          P1<-c(rep(0,lin),rep(diag(Q_start)^(-1),n))
          P1<-diag(P1)
        }else{
          Q_inv.start<-list()
          Q_inv.start[[1]]<-chol2inv(chol(Q_start[1:s[1],1:s[1]]))
          P1<-matrix(0,lin+n%*%s,lin+n%*%s)
          for(jf in 1:n[1])
            P1[(lin+(jf-1)*s[1]+1):(lin+jf*s[1]),(lin+(jf-1)*s[1]+1):(lin+jf*s[1])]<-Q_inv.start[[1]]
          
          for (zu in 2:rnd.len)
          {
            Q_inv.start[[zu]]<-chol2inv(chol(Q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
            for(jf in 1:n[zu])
              P1[(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                 (lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.start[[zu]]
          }
        }
      }
      
      
      if(is.null(family$multivariate)){
      score_vec<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[1,]
      }else{
      score_vec<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[1,]
      }  
      lambda.max<-max(abs(score_vec[(q+1):lin]))

      if (BLOCK)
      {
        grad.1<-gradient.lasso.block(score.beta=score_vec[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda,block=block)
      }else{
        grad.1<-gradient.lasso(score.beta=score_vec[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda)
      }
      
      score_vec[(q+1):lin] <- grad.1
      
      crit.obj<-t_change(grad=score_vec[(q+1):lin],b=Delta[1,(q+1):lin])
      t_edge<-crit.obj$min.rate
      
      optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta_start[1:lin],ranef=Delta_start[(lin+1):(lin+n%*%s)],
                                         Grad=score_vec,family=family,P=diag(P1[(lin+1):(lin+n%*%s)]), 
                                         lower = 0, upper = Inf,K=K))
      
      t_opt<-optim.obj$par
      
      nue<-control$nue
      
      half.index<-0
      solve.test<-FALSE
      ######### big while loop for testing if the update leads to Fisher matrix which can be inverted
      while(!solve.test)
      {  
        
        solve.test2<-FALSE  
        while(!solve.test2)
        {  
          
          if(half.index>50)
          {
            stop("Fisher matrix not invertible")
          }
          
          Delta[1,]<-Delta_start+min(t_opt,t_edge)*nue*(0.5^half.index)*score_vec
          if(t_opt>t_edge & half.index==0)
            Delta[1,crit.obj$whichmin+q]<-0  
          Eta<-Z_alles%*%Delta[1,]

          if(is.null(family$multivariate)){
            D<-family$mu.eta(Eta)
            Mu<-family$linkinv(Eta)
            SigmaInv <- 1/family$variance(Mu)
          }else{
            D <- family$deriv.mat(Eta, K)
            Mu<-family$linkinv(Eta, K)
            SigmaInv <- family$SigmaInv(Mu, K)
          }
          
          ranef.logLik<- -0.5*t(Delta[1,(lin+1):(lin+n%*%s)])%*%(P1[(lin+1):(lin+n%*%s),(lin+1):(lin+n%*%s)]%*%Delta[1,(lin+1):(lin+n%*%s)])
          
          logLik.vec[1]<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu,ranef.logLik=ranef.logLik,family=family,penal=T,K=K, phi = phi)
          
          active<-c(rep(T,q),!is.element(Delta[1,(q+1):lin],0),rep(T,n%*%s))
          Z_aktuell<-Z_alles[,active]
          lin_akt<-q+sum(!is.element(Delta[1,(q+1):lin],0))
          
          if (control$method=="EM")
          {  
            if(rnd.len==1)
            {
              if(s==1)
              {
                P_akt<-diag(c(rep(0,lin_akt),rep(1/Q_start,n*s)))
                }else{
                P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
                for(jf in 1:n)
                  P_akt[(lin_akt+(jf-1)*s+1):(lin_akt+jf*s),(lin_akt+(jf-1)*s+1):(lin_akt+jf*s)]<-Q_inv.start
              }
            }else{
              if(all(s==1))
              {
                P_akt<-diag(c(rep(0,lin_akt),rep(diag(Q_start)^(-1),n)))
              }else{
                P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
                for(jf in 1:n[1])
                  P_akt[(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1]),(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1])]<-Q_inv.start[[1]]
                
                for (zu in 2:rnd.len)
                {
                  for(jf in 1:n[zu])
                    P_akt[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                          (lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.start[[zu]]
                }
              }
            }
            if(is.null(family$multivariate)){
              D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
              F_gross<-t(Z_aktuell)%*%(Z_aktuell*D*SigmaInv*D)+P_akt
            }else{
              W_opt <- RcppEigenProd2(D, SigmaInv)
              F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
            }
            
            InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
            if(inherits(InvFisher2, "try-error"))
              InvFisher2<-try(solve(F_gross),silent=T)
            if(inherits(InvFisher2, "try-error"))
            {
              #stop("Fisher matrix not invertible")  
              half.index<-half.index+1  
            }else{
              solve.test2<-TRUE 
            }}else{
              solve.test2<-TRUE
            }
        }
        
        betaact<-Delta[1,active]
        
        if (control$method=="EM")
        {   
          ############################# Q update ################
          if(rnd.len==1)
          {
            Q1<-InvFisher2[(lin_akt+1):(lin_akt+s),(lin_akt+1):(lin_akt+s)]+Delta[1,(lin+1):(lin+s)]%*%t(Delta[1,(lin+1):(lin+s)])
            for (i in 2:n)
              Q1<-Q1+InvFisher2[(lin_akt+(i-1)*s+1):(lin_akt+i*s),(lin_akt+(i-1)*s+1):(lin_akt+i*s)]+Delta[1,(lin+(i-1)*s+1):(lin+i*s)]%*%t(Delta[1,(lin+(i-1)*s+1):(lin+i*s)])
            Q1<-1/n*Q1
          }else{
            Q1<-matrix(0,sum(s),sum(s))
            Q1[1:s[1],1:s[1]]<-InvFisher2[(lin_akt+1):(lin_akt+s[1]),(lin_akt+1):(lin_akt+s[1])]+Delta[1,(lin+1):(lin+s[1])]%*%t(Delta[1,(lin+1):(lin+s[1])])
            for (i in 2:n[1])
              Q1[1:s[1],1:s[1]]<-Q1[1:s[1],1:s[1]]+InvFisher2[(lin_akt+(i-1)*s[1]+1):(lin_akt+i*s[1]),(lin_akt+(i-1)*s[1]+1):(lin_akt+i*s[1])]+Delta[1,(lin+(i-1)*s[1]+1):(lin+i*s[1])]%*%t(Delta[1,(lin+(i-1)*s[1]+1):(lin+i*s[1])])
            Q1[1:s[1],1:s[1]]<-1/n[1]*Q1[1:s[1],1:s[1]]
            
            for (zu in 2:rnd.len)
            {
              Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-InvFisher2[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu]),(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]+Delta[1,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]%*%t(Delta[1,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])])
              for (i in 2:n[zu])
                Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]+InvFisher2[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu]),(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]+Delta[1,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]%*%t(Delta[1,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])])
              Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-1/n[zu]*Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
            }
          }
        }else{
          if(is.null(family$multivariate)){
          Eta_tilde<-Eta+(y-Mu)*1/D
          }else{
            Eta_tilde<-Eta+solve(D)%*%(y-Mu)
          }
          Betadach<-Delta[1,1:(lin)]     
          aktuell_vec<-!is.element(Delta[1,1:(lin)],0)
          X_aktuell<-Z_fastalles[,aktuell_vec]
 
          if(rnd.len==1)
          {
            
            if(s==1)
            {
              upp<-min(20,50*Q_start)
              low<-1e-14
              optim.obj<-nlminb(sqrt(Q_start),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=Z_fastalles,X_aktuell=X_aktuell,Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp)
              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)))
              optim.obj<-try(bobyqa(q_start_vec,likelihood,D=D,SigmaInv=SigmaInv,family=family,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])
              
              #### Check for positive definitness ########
              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) & !inherits(Q_solvetest, "try-error"))
                  break
              }
            }   
          }else{
            if(all(s==1))
            {
              q_start_vec<-diag(q_start)
              upp<-rep(min(20,50*diag(q_start)),sum(s))
              low<-rep(1e-14,sum(s))
              optim.obj<-try(bobyqa(sqrt(q_start_vec),likelihood_diag,D=D,SigmaInv=SigmaInv,family=family,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,rnd.len=rnd.len))
              Q1<-diag(optim.obj$par)^2
            }else{
              q_start_vec<-c(diag(q_start)[1:s[1]],q_start[1:s[1],1:s[1]][lower.tri(q_start[1:s[1],1:s[1]])])
              up1<-min(20,50*max(q_start_vec))
              low<-c(rep(0,s[1]),rep(-up1,0.5*(s[1]^2-s[1])))
              
              for (zu in 2:rnd.len)
              {
                q_start_vec<-c(q_start_vec,c(diag(q_start)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])],q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]))
                up1<-min(20,50*max(q_start_vec))
                low<-c(low,c(rep(0,s[zu]),rep(-up1,0.5*(s[zu]^2-s[zu]))))
              }
              upp<-rep(up1,length(q_start_vec))
              optim.obj<-try(bobyqa(q_start_vec,likelihood_block,D=D,SigmaInv=SigmaInv,family=family,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,rnd.len=rnd.len))
              optim.vec<-optim.obj$par
              
              
              Q1<-matrix(0,sum(s),sum(s))
              diag(Q1)[1:s[1]]<-optim.vec[1:s[1]]
              if(s[1]>1)
                Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])]<-optim.vec[(s[1]+1):(s[1]*(s[1]+1)*0.5)]
              optim.vec<-optim.vec[-c(1:(s[1]*(s[1]+1)*0.5))]
              
              for (zu in 2:rnd.len)
              {
                diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-optim.vec[1:s[zu]]
                if(s[zu]>1)
                  Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]<-optim.vec[(s[zu]+1):(s[zu]*(s[zu]+1)*0.5)]
                optim.vec<-optim.vec[-c(1:(s[zu]*(s[zu]+1)*0.5))]
              }
              
              #### Check for positive definitness ########
              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) & !inherits(Q_solvetest, "try-error"))
                  break
              }
            }
          }}
        
        Q[[2]]<-Q1
        
        NRstep<-F
        
        if(rnd.len==1)
        {
          if(s==1)
          {
            P1<-c(rep(0,lin),rep(1/Q1,n*s))
            P1<-diag(P1)
          }else{
            Q_inv<-solve(Q1)
            P1<-matrix(0,lin+n%*%s,lin+n%*%s)
            for(j in 1:n)
              P1[(lin+(j-1)*s+1):(lin+j*s),(lin+(j-1)*s+1):(lin+j*s)]<-Q_inv
          }
        }else{
          if(all(s==1))
          {
            P1<-c(rep(0,lin),rep(diag(Q1)^(-1),n))
            P1<-diag(P1)
          }else{
            Q_inv<-list()
            Q_inv[[1]]<-chol2inv(chol(Q1[1:s[1],1:s[1]]))
            P1<-matrix(0,lin+n%*%s,lin+n%*%s)
            for(jf in 1:n[1])
              P1[(lin+(jf-1)*s[1]+1):(lin+jf*s[1]),(lin+(jf-1)*s[1]+1):(lin+jf*s[1])]<-Q_inv[[1]]
            
            for (zu in 2:rnd.len)
            {
              Q_inv[[zu]]<-chol2inv(chol(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
              for(jf in 1:n[zu])
                P1[(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                   (lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
            }
          }
        }
        
        if(is.null(family$multivariate)){
          score_old2<-score_vec2<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[1,]
        }else{
          score_old2<-score_vec2<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[1,]
        }  
        
        lambda.max<-max(abs(score_vec2[(q+1):lin]))
        
        if(BLOCK)
        {
          grad.1<-gradient.lasso.block(score.beta=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda,block=block)
        }else{
          grad.1<-gradient.lasso(score.beta=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda)
        }
        
        score_vec2[(q+1):lin] <- grad.1
        
        crit.obj<-t_change(grad=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin])
        t_edge<-crit.obj$min.rate
        
        optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta[1,1:lin],ranef=Delta[1,(lin+1):(lin+n%*%s)],
                                           Grad=score_vec2,family=family,P=diag(P1[(lin+1):(lin+n%*%s)]), 
                                           lower = 0, upper = Inf,K=K))
        
        t_opt<-optim.obj$par
        
        if(!NRstep && !(all(active_old==active)))
           NRstep <- T
           
           tryNR <- (t_opt<t_edge)  && NRstep  
        
        vorz <- T
        
        if(tryNR) 
        {
          lin_akt<-q+sum(!is.element(Delta[1,(q+1):lin],0))
          if(rnd.len==1)
          {
            if(s==1)
            {
              P_akt<-diag(c(rep(0,lin_akt),rep(1/Q1,n*s)))
            }else{
              P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
              for(jf in 1:n)
                P_akt[(lin_akt+(jf-1)*s+1):(lin_akt+jf*s),(lin_akt+(jf-1)*s+1):(lin_akt+jf*s)]<-Q_inv
            }
          }else{
            if(all(s==1))
            {
              P_akt<-diag(c(rep(0,lin_akt),rep(diag(Q1)^(-1),n)))
            }else{
              P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
              for(jf in 1:n[1])
                P_akt[(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1]),(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1])]<-Q_inv[[1]]
              
              for (zu in 2:rnd.len)
              {
                for(jf in 1:n[zu])
                  P_akt[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                        (lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
              }
            }
          }
          if(is.null(family$multivariate)){
            D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
            F_gross<-t(Z_aktuell)%*%(Z_aktuell*D*SigmaInv*D)+P_akt
          }else{
            W_opt <- RcppEigenProd2(D, SigmaInv)
            F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
          }

          InvFisher<-try(chol2inv(chol(F_gross)),silent=T)
          if(inherits(InvFisher, "try-error"))
            InvFisher<-try(solve(F_gross),silent=T)
          if(inherits(InvFisher, "try-error"))
          {
            half.index<-half.index+1  
          }else{
            solve.test<-TRUE 
            Delta.test<-Delta[1,active]+nue*InvFisher%*%score_old2[active]
            if(lin_akt>q)
            {
              vorz<-all(sign(Delta.test[(q+1):lin_akt])==sign(betaact[(q+1):lin_akt]))  
            }else{
              vorz<-T  
            }
          }
        }else{
          solve.test<-TRUE  
        }
      }   
      
      Eta.ma[2,]<-Eta
      
      score_old<-score_old2
      score_vec<-score_vec2
      Q1.old<-Q1
      Q_inv.old<-Q_inv    
      
      Q1.very.old<-Q_start
      Q_inv.very.old<-Q_inv.start
      
      active_old<-active
      ###############################################################################################################################################
      ################################################################### Main Iteration ###################################################################
      
      if(control$steps!=1)
      {
        for (l in 2:control$steps)
        {
          if(control$print.iter)
          #  message("Iteration ",l)
          {
            cat(if(ia) "\r" else if(l > 1) "\n" else NULL)
            cat(paste("Iteration ",l))
            if(.Platform$OS.type != "unix" & ia) flush.console()
          }
          
          
          if(!vorz)
            tryNR<-F

          half.index<-0
          solve.test<-FALSE
          ######### big while loop for testing if the update leads to Fisher matrix which can be inverted
          while(!solve.test)
          {  
            
            solve.test2<-FALSE  
            while(!solve.test2)
            {  
              if(half.index>50)
              {
                half.index<-Inf;Q1.old<-Q1.very.old;Q_inv.old<-Q_inv.very.old
              }
              
              if(tryNR)
              {
                Delta[l,active_old]<-Delta[l-1,active_old]+nue*(0.5^half.index)*InvFisher%*%score_old[active_old]
                NRstep<-T
#                print("NR")
              }else{
                Delta[l,]<-Delta[l-1,]+min(t_opt,t_edge)*nue*(0.5^half.index)*nue*score_vec
                if(t_opt>t_edge & half.index==0)
                  Delta[l,crit.obj$whichmin+q]<-0  
                NRstep<-F
              }
              
              Eta<-Z_alles%*%Delta[l,]
               if(is.null(family$multivariate)){
                D<-family$mu.eta(Eta)
                Mu<-family$linkinv(Eta)
                SigmaInv <- 1/family$variance(Mu)
              }else{
                D <- family$deriv.mat(Eta, K)
                Mu<-family$linkinv(Eta, K)
                SigmaInv <- family$SigmaInv(Mu, K)
              }
              
              ranef.logLik<- -0.5*t(Delta[l,(lin+1):(lin+n%*%s)])%*%(P1[(lin+1):(lin+n%*%s),(lin+1):(lin+n%*%s)]%*%Delta[l,(lin+1):(lin+n%*%s)])
              
              logLik.vec[l]<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu,ranef.logLik=ranef.logLik,family=family,penal=T,K=K, phi = phi)
              
              active<-c(rep(T,q),!is.element(Delta[l,(q+1):lin],0),rep(T,n%*%s))
              Z_aktuell<-Z_alles[,active]
              lin_akt<-q+sum(!is.element(Delta[l,(q+1):lin],0))
              
              if(control$method=="EM")
              {  
                if(rnd.len==1)
                {
                  if(s==1)
                  {
                    P_akt<-diag(c(rep(0,lin_akt),rep(1/Q1.old,n*s)))
                  }else{
                    P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
                    for(jf in 1:n)
                      P_akt[(lin_akt+(jf-1)*s+1):(lin_akt+jf*s),(lin_akt+(jf-1)*s+1):(lin_akt+jf*s)]<-Q_inv.old
                  }
                }else{
                  if(all(s==1))
                  {
                    P_akt<-diag(c(rep(0,lin_akt),rep(diag(Q1.old)^(-1),n)))
                  }else{
                    P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
                    for(jf in 1:n[1])
                      P_akt[(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1]),(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1])]<-Q_inv.old[[1]]
                    
                    for (zu in 2:rnd.len)
                    {
                      for(jf in 1:n[zu])
                        P_akt[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                              (lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.old[[zu]]
                    }
                  }
                }
                if(is.null(family$multivariate)){
                  D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
                  W_opt <- D*SigmaInv*D
                  F_gross<-t(Z_aktuell)%*%(Z_aktuell*W_opt)+P_akt
                }else{
                  W_opt <- RcppEigenProd2(D, SigmaInv)
                  F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
                }
                
                InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
                if(inherits(InvFisher2, "try-error"))
                  InvFisher2<-try(solve(F_gross),silent=T)
                if(inherits(InvFisher2, "try-error"))
                {
                  half.index<-half.index+1  
                }else{
                  solve.test2<-TRUE 
                }}else{
                  solve.test2<-TRUE 
                }
            }
            
            betaact<-Delta[l,active]
            
            if(control$method=="EM")
            {        
              ############################# Q update ################
              if(rnd.len==1)
              {
                Q1<-InvFisher2[(lin_akt+1):(lin_akt+s),(lin_akt+1):(lin_akt+s)]+Delta[l,(lin+1):(lin+s)]%*%t(Delta[l,(lin+1):(lin+s)])
                for (i in 2:n)
                  Q1<-Q1+InvFisher2[(lin_akt+(i-1)*s+1):(lin_akt+i*s),(lin_akt+(i-1)*s+1):(lin_akt+i*s)]+Delta[l,(lin+(i-1)*s+1):(lin+i*s)]%*%t(Delta[l,(lin+(i-1)*s+1):(lin+i*s)])
                Q1<-1/n*Q1
              }else{
                Q1<-matrix(0,sum(s),sum(s))
                Q1[1:s[1],1:s[1]]<-InvFisher2[(lin_akt+1):(lin_akt+s[1]),(lin_akt+1):(lin_akt+s[1])]+Delta[l,(lin+1):(lin+s[1])]%*%t(Delta[l,(lin+1):(lin+s[1])])
                for (i in 2:n[1])
                  Q1[1:s[1],1:s[1]]<-Q1[1:s[1],1:s[1]]+InvFisher2[(lin_akt+(i-1)*s[1]+1):(lin_akt+i*s[1]),(lin_akt+(i-1)*s[1]+1):(lin_akt+i*s[1])]+Delta[l,(lin+(i-1)*s[1]+1):(lin+i*s[1])]%*%t(Delta[l,(lin+(i-1)*s[1]+1):(lin+i*s[1])])
                Q1[1:s[1],1:s[1]]<-1/n[1]*Q1[1:s[1],1:s[1]]
                
                for (zu in 2:rnd.len)
                {
                  Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-InvFisher2[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu]),(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]+Delta[l,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]%*%t(Delta[l,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])])
                  for (i in 2:n[zu])
                    Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]+InvFisher2[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu]),(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]+Delta[l,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]%*%t(Delta[l,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])])
                  Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-1/n[zu]*Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
                }
              }  
            }else{
              if(is.null(family$multivariate)){
                Eta_tilde<-Eta+(y-Mu)*1/D
              }else{
                Eta_tilde<-Eta+solve(D)%*%(y-Mu)
              }
              Betadach<-Delta[l,1:(lin)]
              
              aktuell_vec<-!is.element(Delta[l,1:(lin)],0)
              X_aktuell<-Z_fastalles[,aktuell_vec]
              
              if(rnd.len==1)
              {
                
                if(s==1)
                {
                  if(Q1<1e-14)
                    low<-0
                  
                  optim.obj<-nlminb(sqrt(Q1),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=Z_fastalles,X_aktuell=X_aktuell,Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp)
                  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,SigmaInv=SigmaInv,family=family,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])
                  
                  #### Check for positiv definitness ########
                  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) & !inherits(Q_solvetest, "try-error"))
                      break
                  }
                }   
              }else{
                if(all(s==1))
                {
                  Q1_vec<-diag(Q1)
                  optim.obj<-try(bobyqa(sqrt(Q1_vec),likelihood_diag,D=D,SigmaInv=SigmaInv,family=family,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,rnd.len=rnd.len))
                  Q1<-diag(optim.obj$par)^2
                }else{
                  Q1_vec<-c(diag(Q1)[1:s[1]],Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])])
                  
                  for (zu in 2:rnd.len)
                    Q1_vec<-c(Q1_vec,c(diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])],Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]))
                  
                  optim.obj<-try(bobyqa(Q1_vec,likelihood_block,D=D,SigmaInv=SigmaInv,family=family,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,rnd.len=rnd.len))
                  optim.vec<-optim.obj$par
                  
                  Q1<-matrix(0,sum(s),sum(s))
                  diag(Q1)[1:s[1]]<-optim.vec[1:s[1]]
                  if(s[1]>1)
                    Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])]<-optim.vec[(s[1]+1):(s[1]*(s[1]+1)*0.5)]
                  optim.vec<-optim.vec[-c(1:(s[1]*(s[1]+1)*0.5))]
                  
                  for (zu in 2:rnd.len)
                  {
                    diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-optim.vec[1:s[zu]]
                    if(s[zu]>1)
                      Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]<-optim.vec[(s[zu]+1):(s[zu]*(s[zu]+1)*0.5)]
                    optim.vec<-optim.vec[-c(1:(s[zu]*(s[zu]+1)*0.5))]
                  }
                  
                  #### Check for positive definitness ########
                  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) & !inherits(Q_solvetest, "try-error"))
                      break
                  }
                }
              }}
            
            Q[[l+1]]<-Q1
            
            if(rnd.len==1)
            {
              if(s==1)
              {
                P1<-c(rep(0,lin),rep(1/Q1,n*s))
                P1<-diag(P1)
              }else{
                Q_inv<-solve(Q1)
                P1<-matrix(0,lin+n%*%s,lin+n%*%s)
                for(j in 1:n)
                  P1[(lin+(j-1)*s+1):(lin+j*s),(lin+(j-1)*s+1):(lin+j*s)]<-Q_inv
              }
            }else{
              if(all(s==1))
              {
                P1<-c(rep(0,lin),rep(diag(Q1)^(-1),n))
                P1<-diag(P1)
              }else{
                Q_inv<-list()
                Q_inv[[1]]<-chol2inv(chol(Q1[1:s[1],1:s[1]]))
                P1<-matrix(0,lin+n%*%s,lin+n%*%s)
                for(jf in 1:n[1])
                  P1[(lin+(jf-1)*s[1]+1):(lin+jf*s[1]),(lin+(jf-1)*s[1]+1):(lin+jf*s[1])]<-Q_inv[[1]]
                
                for (zu in 2:rnd.len)
                {
                  Q_inv[[zu]]<-chol2inv(chol(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
                  for(jf in 1:n[zu])
                    P1[(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                       (lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
                }
              }
            }
            
            if(is.null(family$multivariate)){
              score_old2<-score_vec2<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[l,]
            }else{
              score_old2<-score_vec2<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[l,]
            }  
            lambda.max<-max(abs(score_vec2[(q+1):lin]))
            
            if (BLOCK)
            {
              grad.1<-gradient.lasso.block(score.beta=score_vec2[(q+1):lin],b=Delta[l,(q+1):lin],lambda.b=lambda,block=block)
            }else{
              grad.1<-gradient.lasso(score.beta=score_vec2[(q+1):lin],b=Delta[l,(q+1):lin],lambda.b=lambda)
            }
            score_vec2[(q+1):lin] <- grad.1
            
            crit.obj<-t_change(grad=score_vec2[(q+1):lin],b=Delta[l,(q+1):lin])
            t_edge<-crit.obj$min.rate
            
            optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta[l,1:lin],ranef=Delta[l,(lin+1):(lin+n%*%s)],
                                               Grad=score_vec2,family=family,P=diag(P1[(lin+1):(lin+n%*%s)]), 
                                               lower = 0, upper = Inf,K=K))
            
            t_opt<-optim.obj$par
            
            if(!NRstep && !(all(active_old==active)))
              NRstep <- T
            
            tryNR <- (t_opt<t_edge)  && NRstep  
            
            if(tryNR) 
            {
              lin_akt<-q+sum(!is.element(Delta[l,(q+1):lin],0))
              if(rnd.len==1)
              {
                if(s==1)
                {
                  P_akt<-diag(c(rep(0,lin_akt),rep((Q1^(-1)),n*s)))
                }else{
                  P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
                  for(jf in 1:n)
                    P_akt[(lin_akt+(jf-1)*s+1):(lin_akt+jf*s),(lin_akt+(jf-1)*s+1):(lin_akt+jf*s)]<-Q_inv
                }
              }else{
                if(all(s==1))
                {
                  P_akt<-diag(c(rep(0,lin_akt),rep(diag(Q1)^(-1),n)))
                }else{
                  P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
                  for(jf in 1:n[1])
                    P_akt[(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1]),(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1])]<-Q_inv[[1]]
                  
                  for (zu in 2:rnd.len)
                  {
                    for(jf in 1:n[zu])
                      P_akt[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                            (lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
                  }
                }
              }
              if(is.null(family$multivariate)){
                D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
                F_gross<-t(Z_aktuell)%*%(Z_aktuell*D*SigmaInv*D)+P_akt
              }else{
                W_opt <- RcppEigenProd2(D, SigmaInv)
                F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
              }

              InvFisher3<-try(chol2inv(chol(F_gross)),silent=T)
              if(inherits(InvFisher3, "try-error"))
                InvFisher3<-try(solve(F_gross),silent=T)
              if(inherits(InvFisher3, "try-error"))
              {
                half.index<-half.index+1  
              }else{
                
                solve.test<-TRUE 
                Delta.test<-Delta[l,active]+nue*InvFisher3%*%score_old2[active]
                if(lin_akt>q)
                {
                  vorz<-all(sign(Delta.test[(q+1):lin_akt])==sign(betaact[(q+1):lin_akt]))  
                }else{
                  vorz<-T  
                }
              }
            }else{
              solve.test<-TRUE  
            }
          }  
          
          if(tryNR) 
            InvFisher<-InvFisher3
          
          score_old<-score_old2
          score_vec<-score_vec2
          Q1.very.old<-Q1.old
          Q_inv.very.old<-Q_inv.old
          Q1.old<-Q1
          Q_inv.old<-Q_inv    
          
          Eta.ma[l+1,]<-Eta
          active_old<-active

          finish<-(sqrt(sum((Eta.ma[l,]-Eta.ma[l+1,])^2))/sqrt(sum((Eta.ma[l,])^2))<control$epsilon)
          finish2<-(sqrt(sum((Eta.ma[l-1,]-Eta.ma[l+1,])^2))/sqrt(sum((Eta.ma[l-1,])^2))<control$epsilon)
          if(finish ||  finish2) #|| (all(grad.1 == 0) ))
            break
          Eta.old<-Eta
        }}
      ############################################################################

        if(rnd.len==1)
        {
          if(s==1)
          {
            P_akt<-diag(c(rep(0,lin_akt),rep(1/Q1.old,n*s)))
          }else{
            P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
            for(jf in 1:n)
              P_akt[(lin_akt+(jf-1)*s+1):(lin_akt+jf*s),(lin_akt+(jf-1)*s+1):(lin_akt+jf*s)]<-Q_inv.old
          }
        }else{
          if(all(s==1))
          {
            P_akt<-diag(c(rep(0,lin_akt),rep(diag(Q1.old)^(-1),n)))
          }else{
            P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
            for(jf in 1:n[1])
              P_akt[(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1]),(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1])]<-Q_inv.old[[1]]
            
            for (zu in 2:rnd.len)
            {
              for(jf in 1:n[zu])
                P_akt[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                      (lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.old[[zu]]
            }
          }
        }

        if(is.null(family$multivariate)){
          if(control$method!="EM")
          {
            D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
            W_opt <- D*SigmaInv*D
          }  
          F_gross <- t(Z_aktuell)%*%(Z_aktuell*W_opt)+P_akt
         }else{
           if(control$method!="EM")
             W_opt <- RcppEigenProd2(D, SigmaInv)
          W_inv_t <- RcppEigenSpChol(W_opt)
          F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
        }

        InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
        if(inherits(InvFisher2, "try-error"))
          InvFisher2<-try(solve(F_gross),silent=T)

        if(is.null(family$multivariate)){
          Z_temp <- Z_aktuell*sqrt(W_opt)
          FinalHat.df<-Z_temp%*%InvFisher2%*%t(Z_temp)
        }else{
          FinalHat.df <- RcppEigenProd3(W_inv_t, Z_aktuell, InvFisher2)
        }
        df<-sum(diag(FinalHat.df))
      if(control$overdispersion)
        phi<-(sum((y-Mu)^2/family$variance(Mu)))/(N-df)

      conv.step<-l
      phi.med<-phi
      
      if(conv.step==control$steps)
      {
        cat("Warning:\n")
        cat("Algorithm did not converge!\n")
      }
      
      Delta_neu<-Delta[l,]
      Mu_opt<-Mu
      Qfinal<-Q[[l+1]]
      
      aaa<-!is.element(Delta_neu[1:(lin)],0)
      aaa <- correct.cat(aaa,c(rep(1,sum(is.na(index.new))),block))
      
      if(final.re)
      {    
        ############ final re-estimation
        
        if(rnd.len==1 && s==1)
        {  
          Q.max<-max(sqrt(unlist(Q)))
          Q.min<-min(sqrt(unlist(Q)))
        }else{
          Q.max<-max(Qfinal)+1
          Q.min<-min(Qfinal)-1e-10
        }
        if(rnd.len==1)
        {
          glmm_fin<-try(glmm_final(y,Z_fastalles[,aaa],W,k,n,q_start=Qfinal,K=K,
                                   Delta_start=Delta_neu[c(aaa,rep(T,n%*%s))],s,steps=control$maxIter,
                                   family=family,method=control$method.final,overdispersion=control$overdispersion,
                                   phi=phi,print.iter.final=control$print.iter.final,flushit=control$flushit,eps.final=control$eps.final,
                                   Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
          
          if(inherits(glmm_fin, "try-error") || glmm_fin$opt>control$maxIter-10)
          {  
            glmm_fin2<-try(glmm_final(y,Z_fastalles[,aaa],W,k,n,q_start=q_start,K=K,
                                      Delta_start=Delta_start[c(aaa,rep(T,n%*%s))],s,steps=control$maxIter,
                                      family=family,method=control$method.final,overdispersion=control$overdispersion,
                                      phi=control$phi,print.iter.final=control$print.iter.final,flushit=control$flushit,
                                      eps.final=control$eps.final,Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
            if(!inherits(glmm_fin2, "try-error"))
            {    
              if(inherits(glmm_fin, "try-error") || (glmm_fin$opt>control$maxIter-10 && glmm_fin2$opt<control$maxIter-10))
                glmm_fin<-glmm_fin2 
            }
          }
        }else{
          glmm_fin<-try(glmm_final_multi_random(y,Z_fastalles[,aaa],W,k,q_start=Qfinal,K=K,
                                                Delta_start=Delta_neu[c(aaa,rep(T,n%*%s))],s,n,steps=control$maxIter,
                                                family=family,method=control$method.final,overdispersion=control$overdispersion,
                                                phi=phi,rnd.len=rnd.len,print.iter.final=control$print.iter.final,flushit=control$flushit,
                                                eps.final=control$eps.final,Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
          if(inherits(glmm_fin, "try-error")|| glmm_fin$opt>control$maxIter-10)
          {  
            glmm_fin2<-try(glmm_final_multi_random(y,Z_fastalles[,aaa],W,k,q_start=q_start,K=K,
                                                   Delta_start=Delta_start[c(aaa,rep(T,n%*%s))],s,n,steps=control$maxIter,
                                                   family=family,method=control$method.final,overdispersion=control$overdispersion,
                                                   phi=control$phi,rnd.len=rnd.len,print.iter.final=control$print.iter.final,flushit=control$flushit,
                                                   eps.final=control$eps.final,Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)    
            if(!inherits(glmm_fin2, "try-error"))
            {    
              if(inherits(glmm_fin, "try-error") || (glmm_fin$opt>control$maxIter-10 && glmm_fin2$opt<control$maxIter-10))
                glmm_fin<-glmm_fin2 
            }
          }
        }
        
        
        if(inherits(glmm_fin, "try-error") || glmm_fin$opt>control$maxIter-10)
        {
          cat("\nWarning:\n")
          cat("Final Fisher scoring reestimation did not converge!\n")
        }
        #######
      }else{
        glmm_fin<-NA  
        class(glmm_fin)<-"try-error"
      }  
      Standard_errors<-matrix(NA,length(Delta_neu),length(Delta_neu))
      
      if(!inherits(glmm_fin, "try-error"))
      {
        Delta_neu2<-Delta_neu
        Delta_neu2[c(aaa,rep(T,n%*%s))]<-glmm_fin$Delta
        Standard_errors[c(aaa,rep(T,n%*%s)),c(aaa,rep(T,n%*%s))]<-glmm_fin$Standard_errors
        Qfinal<-glmm_fin$Q
        phi<-glmm_fin$phi
        complexity<-glmm_fin$complexity
        Delta_neu<-Delta_neu2
        Eta_opt<-Z_alles%*%Delta_neu
        if(is.null(family$multivariate)){
          Mu_opt<-family$linkinv(Eta_opt)
        }else{
          Mu_opt <- family$linkinv(Eta_opt, K)
        }
      }else{
        glmm_fin<-list()
        glmm_fin$ranef.logLik<-ranef.logLik
        complexity<-df
      }
      
      if(rnd.len==1)
      {
        if(s==1)
          Qfinal<-sqrt(Qfinal)
        
        if(!is.matrix(Qfinal))
          Qfinal<-as.matrix(Qfinal)
        colnames(Qfinal)<-random.labels
        rownames(Qfinal)<-random.labels
      }else{
        Qfinal_old<-Qfinal
        Qfinal<-list()
        Qfinal[[1]]<-as.matrix(Qfinal_old[1:s[1],1:s[1]])
        colnames(Qfinal[[1]])<-random.labels[[1]]
        rownames(Qfinal[[1]])<-random.labels[[1]]
        
        if(s[1]==1)
          Qfinal[[1]]<-sqrt(Qfinal[[1]])
        
        for (zu in 2:rnd.len)
        {
          Qfinal[[zu]]<-as.matrix(Qfinal_old[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])
          colnames(Qfinal[[zu]])<-random.labels[[zu]]
          rownames(Qfinal[[zu]])<-random.labels[[zu]]
          
          if(s[zu]==1)
            Qfinal[[zu]]<-sqrt(Qfinal[[zu]])
          
        }
      }

      colnames(Standard_errors) <- rownames(Standard_errors) <- paste0("help.",1:nrow(Standard_errors))
      
      if(dim(X)[2]>0)
      {  
        names(Delta_neu)[1:dim(X)[2]]<-colnames(X)
        colnames(Standard_errors)[1:dim(X)[2]]<-rownames(Standard_errors)[1:dim(X)[2]]<-colnames(X)
      }
      
      if(lin>1)
      {
        names(Delta_neu)[(dim(X)[2]+1):lin]<-colnames(U)
        colnames(Standard_errors)[(dim(X)[2]+1):lin]<-rownames(Standard_errors)[(dim(X)[2]+1):lin]<-colnames(U)
      }
      names(Delta_neu)[(lin+1):(lin+n%*%s)]<-colnames(W)
      colnames(Standard_errors)[(lin+1):(lin+n%*%s)]<-rownames(Standard_errors)[(lin+1):(lin+n%*%s)]<-colnames(W)
      colnames(Delta)<-c(final.names,colnames(W))

      Delta_neu[1:lin] <- Delta_neu[1:lin][transf.names]
      Standard_errors[1:lin,1:lin] <-Standard_errors[1:lin,1:lin][transf.names,transf.names]
      names(Delta_neu)[1:lin] <- transf.names
      colnames(Standard_errors)[1:lin] <- rownames(Standard_errors)[1:lin] <- transf.names
      
      ## Transform the coefficients back to the original scale if the design
      ## matrix was standardized
      if(standardize){
        if(any.notpen)
        {
          Delta_neu[inotpen.which] <- (1 / scale.notpen) * Delta_neu[inotpen.which]
          if(!inherits(glmm_fin, "try-error"))
            Standard_errors[inotpen.which] <- (1 / scale.notpen) * Standard_errors[inotpen.which]
        }
        ## For df > 1 we have to use a matrix inversion to go back to the
        ## original scale
        for(j in 1:length(ipen.which)){
          ind <- ipen.which[[j]]
          Delta_neu[ind] <- solve(scale.pen[[j]], Delta_neu[ind,drop = FALSE])
          if(!inherits(glmm_fin, "try-error"))
          {
            Sc.help <- solve(scale.pen[[j]])
            Standard_errors[ind,ind] <- t(Sc.help)%*%(Standard_errors[ind,ind]%*%Sc.help)
          }
        }
        if(any(s>1))
          warning("Random slopes are not standardized back!")
      }
      
      Standard_errors <- sqrt(diag(Standard_errors))
      ## Need to adjust intercept if we have performed centering
      if(center){
        Delta_neu[intercept.which] <- Delta_neu[intercept.which] -
          sum(Delta_neu[1:lin][-intercept.which,drop = FALSE] * mu.x)   
      }
      
      aic<-NaN
      bic<-NaN
      if(is.element(family$family,c("gaussian", "binomial", "poisson","acat","cumulative"))) 
      {
        
        loglik<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu_opt,ranef.logLik=glmm_fin$ranef.logLik,family=family,penal=T,K=K, phi = phi)
        
        if(control$complexity!="hat.matrix")  
        {  
          if(rnd.len==1)
          {
            complexity<-0.5*(s*(s+1))
          }else{
            complexity<-0.5*(s[1]*(s[1]+1))
            for(zu in 2:rnd.len)
              complexity<-complexity+0.5*(s[zu]*(s[zu]+1))
          }
          complexity<-complexity+sum(Delta_neu[1:(lin)]!=0)
        }    
        aic<--2*loglik+2*complexity
        bic<--2*loglik+log(N)*complexity
      }else{
        warning("For the specified family (so far) no AIC and BIC are available!")  
      }
      
      ret.obj=list()
      ret.obj$aic<-aic
      ret.obj$bic<-bic
      ret.obj$Deltamatrix<-Delta
      ret.obj$ranef<-Delta_neu[(lin+1):(lin+n%*%s)]
      ret.obj$coefficients<-Delta_neu[1:(lin)]
      ret.obj$fixerror<-Standard_errors[1:(lin)]
      ret.obj$ranerror<-Standard_errors[(lin+1):(lin+n%*%s)]
      ret.obj$Q_long<-Q
      ret.obj$Q<-Qfinal
      ret.obj$y_hat<-Mu_opt
      ret.obj$phi<-phi
      ret.obj$family<-family
      ret.obj$fix<-fix.old
      ret.obj$conv.step<-conv.step
      ret.obj$newrndfrml<-newrndfrml
      ret.obj$subject<-subject.names
      ret.obj$data<-data
      ret.obj$rnd.len<-rnd.len
      ret.obj$phi.med<-phi.med
      ret.obj$y <- y
      ret.obj$X <- cbind(X,U)
      ret.obj$W <- W
      ret.obj$K <- K
      ret.obj$df<-df
      ret.obj$loglik<-loglik
      ret.obj$lambda.max<-lambda.max
      ret.obj$logLik.vec<-logLik.vec
      return(ret.obj)
      ##############################################################  
      ######################## 2. Smooth ###########################  
      ##############################################################  
    }else{
      
      smooth<-control$smooth
      
      if(attr(terms(smooth$formula), "intercept")==0)
      {
        variables<-attr(terms(smooth$formula),"term.labels")
        smooth$formula<- "~ +1"
        for (ir in 1:length(variables))
          smooth$formula <- paste(smooth$formula, variables[ir], sep="+")
        smooth$formula <-formula(smooth$formula)
      }
      
      B <- model.matrix(smooth$formula, data)
      B.names<-attr(B,"dimnames")[[2]]
      
      B<-as.matrix(B[,-1])
      attr(B,"dimnames")[[2]]<-B.names[-1]
      
      nbasis<-smooth$nbasis
      diff.ord<-smooth$diff.ord
      spline.degree<-smooth$spline.degree
      knots.no<-nbasis-1
      if(spline.degree<3 && (spline.degree-diff.ord)<2)
        knots.no<-knots.no+1  
      penal<-smooth$penal
      
      if(!(diff.ord<spline.degree))
        stop("Order of differences must be lower than degree of B-spline polynomials!")
      
      m<-dim(B)[2]
      
      Phi<-numeric()
      
      for (r in 1:m)
      {
        Basis<-bs.design(B[,r],diff.ord=diff.ord,spline.degree=spline.degree,knots.no=knots.no)
        Phi_temp<-cbind(Basis$X[,-1],Basis$Z)
        colnames(Phi_temp)<-paste(colnames(B)[r],rep(1:dim(Phi_temp)[2],each=1), sep=".")
        Phi<-cbind(Phi,Phi_temp)
      }
      
      dim.smooth<-dim(Phi)[2]
      
      if(is.null(smooth$start))
        smooth$start<-rep(0,dim.smooth)  
      
      smooth_null<-smooth$start
      
      if(lin>1)
      {
        Eta_start<-X%*%beta_null+Phi%*%smooth_null+W%*%ranef_null
      }else{
        Eta_start<-rep(beta_null,N)+Phi%*%smooth_null+W%*%ranef_null
      }
      
      if(is.null(family$multivariate)){
        D<-family$mu.eta(Eta_start)
        Mu<-family$linkinv(Eta_start)
        SigmaInv <- 1/family$variance(Mu)
      }else{
        D <- family$deriv.mat(Eta_start, K)
        Mu<-family$linkinv(Eta_start, K)
        SigmaInv <- family$SigmaInv(Mu, K)
      }
      
      if(rnd.len==1)
      {
        if(s==1)
          Q_start<-q_start
      }else{
        if(all(s==1))
        {
          Q_start<-diag(diag(q_start),sum(s))
        }else{
          Q_start<-matrix(0,sum(s),sum(s))
          Q_start[1:s[1],1:s[1]]<-q_start[1:s[1],1:s[1]]
          
          for (zu in 2:rnd.len)
          {
            Q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
          }
        }
      }
      
      U<-X[,!is.na(index.new),drop=FALSE]
      X<-X[,is.na(index.new),drop=FALSE]
      
      final.names<-c(colnames(X),colnames(U))
      
      q<-dim(X)[2]
      
      Z_fastalles<-cbind(X,U)
      Z_alles<-cbind(Z_fastalles,Phi,W)
      
      ########################################################## some definitions ################################################
      Delta<-matrix(0,control$steps,(lin+dim.smooth+n%*%s))
      Delta[1,1:lin]<-beta_null[final.names]
      Delta[1,(lin+1):(lin+dim.smooth)]<-smooth_null
      Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]<-t(ranef_null)
      active_old<-c(rep(T,q),!is.element(Delta[1,(q+1):lin],0),rep(T,dim.smooth+n%*%s))
      
      control$epsilon<-control$epsilon*sqrt(dim(Delta)[2])
      
      Delta_start<-Delta[1,]
      
      logLik.vec<-c()
      
      Eta.ma<-matrix(0,control$steps+1,N)
      Eta.ma[1,]<-Eta_start
      
      Q_inv<-NULL
      Q_inv.old.temp<-NULL
      Q_inv.start<-NULL
      
      Q<-list()
      Q[[1]]<-Q_start
      
      l=1
      
      if(diff.ord>1)
      {
        k2<-c(rep(0,diff.ord-1),rep(1,nbasis-diff.ord+1))
      }else{
        k2<-rep(1,nbasis)
      }
      k22<-rep(k2,m)
      penal.vec<-penal*k22
      
      Q_inv<-NULL
      Q_inv.old.temp<-NULL
      Q_inv.start<-NULL
      
      if(rnd.len==1)
      {
        if(s==1)
        {
          P1<-c(rep(0,lin),penal.vec,rep(1/Q_start,n*s))
          P1<-diag(P1)
        }else{
          Q_inv.start<-chol2inv(chol(Q_start))
          P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
          diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
          for(j in 1:n)
            P1[(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s),(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s)]<-Q_inv
        }
      }else{
        if(all(s==1))
        {
          P1<-c(rep(0,lin),penal.vec,rep(diag(Q_start)^(-1),n))
          P1<-diag(P1)
        }else{
          Q_inv.start<-list()
          Q_inv.start[[1]]<-chol2inv(chol(Q_start[1:s[1],1:s[1]]))
          P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
          diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
          for(jf in 1:n[1])
            P1[(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1]),(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1])]<-Q_inv[[1]]
          
          for (zu in 2:rnd.len)
          {
            Q_inv.start[[zu]]<-chol2inv(chol(Q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
            for(jf in 1:n[zu])
              P1[(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                 (lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
          }
        }
      }
      
      if(is.null(family$multivariate)){
        score_vec<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[1,]
      }else{
        score_vec<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[1,]
      }  
      lambda.max<-max(abs(score_vec[(q+1):lin]))
      
      if (BLOCK)
      {
        grad.1<-gradient.lasso.block(score.beta=score_vec[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda,block=block)
      }else{
        grad.1<-gradient.lasso(score.beta=score_vec[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda)
      }
      
      score_vec[(q+1):lin] <- grad.1
      
      crit.obj<-t_change(grad=score_vec[(q+1):lin],b=Delta[1,(q+1):lin])
      t_edge<-crit.obj$min.rate
      
      optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta_start[1:(lin+dim.smooth)],ranef=Delta_start[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)],
                                         Grad=score_vec,family=family,P=diag(P1[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]), 
                                         lower = 0, upper = Inf,K=K))
      
      t_opt<-optim.obj$par
      
      nue<-control$nue
      
      half.index<-0
      solve.test<-FALSE
      ######### big while loop for testing if the update leads to Fisher matrix which can be inverted
      while(!solve.test)
      {  
        
        solve.test2<-FALSE  
        while(!solve.test2)
        {  
          
          if(half.index>50)
          {
            stop("Fisher matrix not invertible")
          }
          
          Delta[1,]<-Delta_start+min(t_opt,t_edge)*nue*(0.5^half.index)*score_vec
          if(t_opt>t_edge & half.index==0)
            Delta[1,crit.obj$whichmin+q]<-0  
          Eta<-Z_alles%*%Delta[1,]

          if(is.null(family$multivariate)){
            D<-family$mu.eta(Eta)
            Mu<-family$linkinv(Eta)
            SigmaInv <- 1/family$variance(Mu)
          }else{
            D <- family$deriv.mat(Eta, K)
            Mu<-family$linkinv(Eta, K)
            SigmaInv <- family$SigmaInv(Mu, K)
          }          
          ranef.logLik<- -0.5*t(Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)])%*%(P1[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s),(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]%*%Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)])
          
          logLik.vec[1]<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu,ranef.logLik=ranef.logLik,family=family,penal=T,K=K, phi = phi)
          
          active<-c(rep(T,q),!is.element(Delta[1,(q+1):lin],0),rep(T,dim.smooth+n%*%s))
          Z_aktuell<-Z_alles[,active]
          lin_akt<-q+sum(!is.element(Delta[1,(q+1):lin],0))
          
          if (control$method=="EM")
          {  
            if(rnd.len==1)
            {
              if(s==1)
              {
                P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(1/Q_start,n*s)))
              }else{
                P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
                diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
                for(jf in 1:n)
                  P_akt[(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s),(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s)]<-Q_inv.start
              }
            }else{
              if(all(s==1))
              {
                P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(diag(Q_start)^(-1),n)))
              }else{
                P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
                diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
                for(jf in 1:n[1])
                  P_akt[(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1]),(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1])]<-Q_inv.start[[1]]
                
                for (zu in 2:rnd.len)
                {
                  for(jf in 1:n[zu])
                    P_akt[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                          (lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.start[[zu]]
                }
              }
            }
            if(is.null(family$multivariate)){
              D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
              F_gross<-t(Z_aktuell)%*%(Z_aktuell*D*SigmaInv*D)+P_akt
            }else{
              W_opt <- RcppEigenProd2(D, SigmaInv)
              F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
            }
            InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
            if(inherits(InvFisher2, "try-error"))
              InvFisher2<-try(solve(F_gross),silent=T)
            if(inherits(InvFisher2, "try-error"))
            {
              #stop("Fisher matrix not invertible")  
              half.index<-half.index+1  
            }else{
              solve.test2<-TRUE 
            }}else{
              solve.test2<-TRUE
            }
        }
        
        betaact<-Delta[1,active]
        
        if (control$method=="EM")
        {   
          ############################# Q update ################
          if(rnd.len==1)
          {
            Q1<-InvFisher2[(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s),(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s)]+Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+s)]%*%t(Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+s)])
            for (i in 2:n)
              Q1<-Q1+InvFisher2[(lin_akt+dim.smooth+(i-1)*s+1):(lin_akt+dim.smooth+i*s),(lin_akt+dim.smooth+(i-1)*s+1):(lin_akt+dim.smooth+i*s)]+Delta[1,(lin+dim.smooth+(i-1)*s+1):(lin+dim.smooth+i*s)]%*%t(Delta[1,(lin+dim.smooth+(i-1)*s+1):(lin+dim.smooth+i*s)])
            Q1<-1/n*Q1
          }else{
            Q1<-matrix(0,sum(s),sum(s))
            Q1[1:s[1],1:s[1]]<-InvFisher2[(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s[1]),(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s[1])]+Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+s[1])]%*%t(Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+s[1])])
            for (i in 2:n[1])
              Q1[1:s[1],1:s[1]]<-Q1[1:s[1],1:s[1]]+InvFisher2[(lin_akt+dim.smooth+(i-1)*s[1]+1):(lin_akt+dim.smooth+i*s[1]),(lin_akt+dim.smooth+(i-1)*s[1]+1):(lin_akt+dim.smooth+i*s[1])]+Delta[1,(lin+dim.smooth+(i-1)*s[1]+1):(lin+dim.smooth+i*s[1])]%*%t(Delta[1,(lin+dim.smooth+(i-1)*s[1]+1):(lin+dim.smooth+i*s[1])])
            Q1[1:s[1],1:s[1]]<-1/n[1]*Q1[1:s[1],1:s[1]]
            
            for (zu in 2:rnd.len)
            {
              Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-InvFisher2[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu]),(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]+Delta[1,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]%*%t(Delta[1,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])])
              for (i in 2:n[zu])
                Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]+InvFisher2[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu]),(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]+Delta[1,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]%*%t(Delta[1,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])])
              Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-1/n[zu]*Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
            }
          }
        }else{
          if(is.null(family$multivariate)){
            Eta_tilde<-Eta+(y-Mu)*1/D
          }else{
            Eta_tilde<-Eta+solve(D)%*%(y-Mu)
          }
          Betadach<-Delta[1,1:(lin+dim.smooth)]     
          aktuell_vec<-!is.element(Delta[1,1:(lin)],0)
          X_aktuell<-Z_fastalles[,aktuell_vec]
          
          if(rnd.len==1)
          {
            if(s==1)
            {
              upp<-min(20,50*Q_start)
              low<-1e-14
              optim.obj<-try(nlminb(sqrt(Q_start),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),
                                    Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp))
              if(inherits(optim.obj, "try-error"))
                optim.obj<-try(bobyqa(sqrt(Q_start),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),
                                      Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp))
              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))#[(s+1):(s*(s+1)*0.5)]))
              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,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),
                                    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])
              
              #### Check for positive definitness ########
              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) & !inherits(Q_solvetest, "try-error"))
                  break
              }
            }   
          }else{
            if(all(s==1))
            {
              q_start_vec<-diag(q_start)
              upp<-rep(min(20,50*diag(q_start)),sum(s))
              low<-rep(1e-14,sum(s))
              optim.obj<-try(bobyqa(sqrt(q_start_vec),likelihood_diag,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,s=s,k=k,Betadach=Betadach,W=W, lower=low,upper=upp,rnd.len=rnd.len))
              Q1<-diag(optim.obj$par)^2
            }else{
              q_start_vec<-c(diag(q_start)[1:s[1]],q_start[1:s[1],1:s[1]][lower.tri(q_start[1:s[1],1:s[1]])])
              up1<-min(20,50*max(q_start_vec))
              low<-c(rep(0,s[1]),rep(-up1,0.5*(s[1]^2-s[1])))
              
              for (zu in 2:rnd.len)
              {
                q_start_vec<-c(q_start_vec,c(diag(q_start)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])],q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]))
                up1<-min(20,50*max(q_start_vec))
                low<-c(low,c(rep(0,s[zu]),rep(-up1,0.5*(s[zu]^2-s[zu]))))
              }
              upp<-rep(up1,length(q_start_vec))
              optim.obj<-try(bobyqa(q_start_vec,likelihood_block,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,s=s,k=k,Betadach=Betadach,W=W, lower=low,upper=upp,rnd.len=rnd.len))
              optim.vec<-optim.obj$par
              
              Q1<-matrix(0,sum(s),sum(s))
              diag(Q1)[1:s[1]]<-optim.vec[1:s[1]]
              if(s[1]>1)
                Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])]<-optim.vec[(s[1]+1):(s[1]*(s[1]+1)*0.5)]
              optim.vec<-optim.vec[-c(1:(s[1]*(s[1]+1)*0.5))]
              
              for (zu in 2:rnd.len)
              {
                diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-optim.vec[1:s[zu]]
                if(s[zu]>1)
                  Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]<-optim.vec[(s[zu]+1):(s[zu]*(s[zu]+1)*0.5)]
                optim.vec<-optim.vec[-c(1:(s[zu]*(s[zu]+1)*0.5))]
              }
              
              #### Check for positive definitness ########
              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) & !inherits(Q_solvetest, "try-error"))
                  break
              }
            }
          }}
        
        Q[[2]]<-Q1

        NRstep<-F
        
        if(rnd.len==1)
        {
          if(s==1)
          {
            P1<-c(rep(0,lin),penal.vec,rep(1/Q1,n*s))
            P1<-diag(P1)
          }else{
            Q_inv<-chol2inv(chol(Q1))
            Q_inv.old.temp<-Q_inv
            P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
            diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
            for(j in 1:n)
              P1[(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s),(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s)]<-Q_inv
          }
        }else{
          if(all(s==1))
          {
            P1<-c(rep(0,lin),penal.vec,rep(diag(Q1)^(-1),n))
            P1<-diag(P1)
          }else{
            Q_inv<-list()
            Q_inv[[1]]<-chol2inv(chol(Q1[1:s[1],1:s[1]]))
            P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
            diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
            for(jf in 1:n[1])
              P1[(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1]),(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1])]<-Q_inv[[1]]
            
            for (zu in 2:rnd.len)
            {
              Q_inv[[zu]]<-chol2inv(chol(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
              for(jf in 1:n[zu])
                P1[(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                   (lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
            }
          }
        }
        
        if(is.null(family$multivariate)){
          score_old2<-score_vec2<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[1,]
        }else{
          score_old2<-score_vec2<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[1,]
        }  
        lambda.max<-max(abs(score_vec2[(q+1):lin]))
        
        if (BLOCK)
        {
          grad.1<-gradient.lasso.block(score.beta=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda,block=block)
        }else{
          grad.1<-gradient.lasso(score.beta=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda)
        }
        score_vec2[(q+1):lin] <- grad.1
        
        crit.obj<-t_change(grad=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin])
        t_edge<-crit.obj$min.rate
        
        optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta[1,1:(lin+dim.smooth)],ranef=Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)],
                                           Grad=score_vec2,family=family,P=diag(P1[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]), 
                                           lower = 0, upper = Inf,K=K))
        
        t_opt<-optim.obj$par
        
        if(!NRstep && !(all(active_old==active)))
          NRstep <- T
        
        tryNR <- (t_opt<t_edge)  && NRstep  
        
        vorz <- T
        
        if(tryNR) 
        {
          lin_akt<-q+sum(!is.element(Delta[1,(q+1):lin],0))
          if(rnd.len==1)
          {
            if(s==1)
            {
              P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(1/Q1,n*s)))
            }else{
              P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
              diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
              for(jf in 1:n)
                P_akt[(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s),
                      (lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s)]<-Q_inv
            }
          }else{
            if(all(s==1))
            {
              P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(diag(Q1)^(-1),n)))
            }else{
              P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
              diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
              for(jf in 1:n[1])
                P_akt[(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1]),
                      (lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1])]<-Q_inv[[1]]
              
              for (zu in 2:rnd.len)
              {
                for(jf in 1:n[zu])
                  P_akt[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                        (lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
              }
            }
          }
          if(is.null(family$multivariate)){
            D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
            W_opt <- D*SigmaInv*D
            F_gross<-t(Z_aktuell)%*%(Z_aktuell*W_opt)+P_akt
          }else{
            W_opt <- RcppEigenProd2(D, SigmaInv)
            F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
          }
          InvFisher<-try(chol2inv(chol(F_gross)),silent=T)
          if(inherits(InvFisher, "try-error"))
            InvFisher<-try(solve(F_gross),silent=T)
          if(inherits(InvFisher, "try-error"))
          {
            half.index<-half.index+1  
          }else{
            solve.test<-TRUE 
            Delta.test<-Delta[1,active]+nue*InvFisher%*%score_vec2[active]
            if(lin_akt>q)
            {
              vorz<-all(sign(Delta.test[(q+1):lin_akt])==sign(betaact[(q+1):lin_akt]))  
            }else{
              vorz<-T  
            }
          }
        }else{
          solve.test<-TRUE  
        }
      }   
      
      Eta.ma[2,]<-Eta
      
      score_old<-score_old2
      score_vec<-score_vec2
      Q1.old<-Q1
      Q_inv.old<-Q_inv    
      
      Q1.very.old<-Q_start
      Q_inv.very.old<-Q_inv.start

      active_old<-active
      
      ###############################################################################################################################################
      ################################################################### Main Iteration ###################################################################
      if(control$steps!=1)
      {
        for (l in 2:control$steps)
        {

          if(control$print.iter)
#            message("Iteration ",l)
          {
            cat(if(ia) "\r" else if(l > 1) "\n" else NULL)
            cat(paste("Iteration ",l))
            if(.Platform$OS.type != "unix" & ia) flush.console()
          }

          if(!vorz)
            tryNR<-F
          
          half.index<-0
          solve.test<-FALSE
          ######### big while loop for testing if the update leads to Fisher matrix which can be inverted
          while(!solve.test)
          {  
            
            solve.test2<-FALSE  
            while(!solve.test2)
            {  
              
              if(half.index>50)
              {
                half.index<-Inf;Q1.old<-Q1.very.old;Q_inv.old<-Q_inv.very.old;
              }
              
              if(tryNR)
              {
                Delta[l,active_old]<- Delta[l-1,active_old]+nue*(0.5^half.index)*InvFisher%*%score_old[active_old]
                NRstep<-T
#                print("NR")
              }else{
                Delta[l,]<-Delta[l-1,]+min(t_opt,t_edge)*nue*(0.5^half.index)*score_vec
                if(t_opt>t_edge & half.index==0)
                  Delta[l,crit.obj$whichmin+q]<-0  
                NRstep<-F
              }
              
              Eta<-Z_alles%*%Delta[l,]
              if(is.null(family$multivariate)){
                D<-family$mu.eta(Eta)
                Mu<-family$linkinv(Eta)
                SigmaInv <- 1/family$variance(Mu)
              }else{
                D <- family$deriv.mat(Eta, K)
                Mu<-family$linkinv(Eta, K)
                SigmaInv <- family$SigmaInv(Mu, K)
              }
              
              ranef.logLik<- -0.5*t(Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)])%*%(P1[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s),(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]%*%Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)])
              
              logLik.vec[l]<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu,ranef.logLik=ranef.logLik,family=family,penal=T,K=K, phi = phi)
              
              active<-c(rep(T,q),!is.element(Delta[l,(q+1):lin],0),rep(T,dim.smooth+n%*%s))
              Z_aktuell<-Z_alles[,active]
              lin_akt<-q+sum(!is.element(Delta[l,(q+1):lin],0))
              
              if (control$method=="EM")
              {  
                if(rnd.len==1)
                {
                  if(s==1)
                  {
                    P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(1/Q1.old,n*s)))
                  }else{
                    P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
                    diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
                    for(jf in 1:n)
                      P_akt[(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s),
                            (lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s)]<-Q_inv.old
                  }
                }else{
                  if(all(s==1))
                  {
                    P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(diag(Q1.old)^(-1),n)))
                  }else{
                    P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
                    diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
                    for(jf in 1:n[1])
                      P_akt[(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1]),
                            (lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1])]<-Q_inv.old[[1]]
                    
                    for (zu in 2:rnd.len)
                    {
                      for(jf in 1:n[zu])
                        P_akt[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                              (lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.old[[zu]]
                    }
                  }
                 }
                if(is.null(family$multivariate)){
                  D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
                  W_opt <- D*SigmaInv*D
                  F_gross<-t(Z_aktuell)%*%(Z_aktuell*W_opt)+P_akt
                }else{
                  W_opt <- RcppEigenProd2(D, SigmaInv)
                  F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
                }
                InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
                if(inherits(InvFisher2, "try-error"))
                  InvFisher2<-try(solve(F_gross),silent=T)
                if(inherits(InvFisher2, "try-error"))
                {
                  half.index<-half.index+1  
                }else{
                  solve.test2<-TRUE 
                }}else{
                  solve.test2<-TRUE 
                }
            }
            
            betaact<-Delta[l,active]
            
            if (control$method=="EM")
            {        
              ############################# Q update ################
              if(rnd.len==1)
              {
                Q1<-InvFisher2[(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s),(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s)]+Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+s)]%*%t(Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+s)])
                for (i in 2:n)
                  Q1<-Q1+InvFisher2[(lin_akt+dim.smooth+(i-1)*s+1):(lin_akt+dim.smooth+i*s),(lin_akt+dim.smooth+(i-1)*s+1):(lin_akt+dim.smooth+i*s)]+Delta[l,(lin+dim.smooth+(i-1)*s+1):(lin+dim.smooth+i*s)]%*%t(Delta[l,(lin+dim.smooth+(i-1)*s+1):(lin+dim.smooth+i*s)])
                Q1<-1/n*Q1
              }else{
                Q1<-matrix(0,sum(s),sum(s))
                Q1[1:s[1],1:s[1]]<-InvFisher2[(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s[1]),(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s[1])]+Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+s[1])]%*%t(Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+s[1])])
                for (i in 2:n[1])
                  Q1[1:s[1],1:s[1]]<-Q1[1:s[1],1:s[1]]+InvFisher2[(lin_akt+dim.smooth+(i-1)*s[1]+1):(lin_akt+dim.smooth+i*s[1]),(lin_akt+dim.smooth+(i-1)*s[1]+1):(lin_akt+dim.smooth+i*s[1])]+Delta[l,(lin+dim.smooth+(i-1)*s[1]+1):(lin+dim.smooth+i*s[1])]%*%t(Delta[l,(lin+dim.smooth+(i-1)*s[1]+1):(lin+dim.smooth+i*s[1])])
                Q1[1:s[1],1:s[1]]<-1/n[1]*Q1[1:s[1],1:s[1]]
                
                for (zu in 2:rnd.len)
                {
                  Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-InvFisher2[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu]),(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]+Delta[l,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]%*%t(Delta[l,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])])
                  for (i in 2:n[zu])
                    Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]+InvFisher2[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu]),(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]+Delta[l,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]%*%t(Delta[l,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])])
                  Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-1/n[zu]*Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
                }
              }  
            }else{
              if(is.null(family$multivariate)){
                Eta_tilde<-Eta+(y-Mu)*1/D
              }else{
                Eta_tilde<-Eta+solve(D)%*%(y-Mu)
              }
              
              Betadach<-Delta[l,1:(lin+dim.smooth)]
              
              aktuell_vec<-!is.element(Delta[l,1:(lin)],0)
              X_aktuell<-Z_fastalles[,aktuell_vec]
              
              if(rnd.len==1)
              {
                
                if(s==1)
                {
                  if(Q1<1e-14)
                    low<-0
                  
                  optim.obj<-try(nlminb(sqrt(Q1),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp))
                  if(inherits(optim.obj, "try-error"))
                    optim.obj<-bobyqa(sqrt(Q1),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp)
                  
                  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,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),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])
                  
                  #### Check for positiv definitness ########
                  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) & !inherits(Q_solvetest, "try-error"))
                      break
                  }
                }   
              }else{
                if(all(s==1))
                {
                  Q1_vec<-diag(Q1)
                  optim.obj<-try(bobyqa(sqrt(Q1_vec),likelihood_diag,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,s=s,k=k,Betadach=Betadach,W=W, lower=low,upper=upp,rnd.len=rnd.len))
                  Q1<-diag(optim.obj$par)^2
                }else{
                  Q1_vec<-c(diag(Q1)[1:s[1]],Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])])
                  
                  for (zu in 2:rnd.len)
                    Q1_vec<-c(Q1_vec,c(diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])],Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]))
                  
                  optim.obj<-try(bobyqa(Q1_vec,likelihood_block,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,s=s,k=k,Betadach=Betadach,W=W, lower=low,upper=upp,rnd.len=rnd.len))
                  optim.vec<-optim.obj$par
                  
                  Q1<-matrix(0,sum(s),sum(s))
                  diag(Q1)[1:s[1]]<-optim.vec[1:s[1]]
                  if(s[1]>1)
                    Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])]<-optim.vec[(s[1]+1):(s[1]*(s[1]+1)*0.5)]
                  optim.vec<-optim.vec[-c(1:(s[1]*(s[1]+1)*0.5))]
                  
                  for (zu in 2:rnd.len)
                  {
                    diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-optim.vec[1:s[zu]]
                    if(s[zu]>1)
                      Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]<-optim.vec[(s[zu]+1):(s[zu]*(s[zu]+1)*0.5)]
                    optim.vec<-optim.vec[-c(1:(s[zu]*(s[zu]+1)*0.5))]
                  }
                  
                  #### Check for positive definitness ########
                  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) & !inherits(Q_solvetest, "try-error"))
                      break
                  }
                }
              }}

            Q[[l+1]]<-Q1
            
            if(rnd.len==1)
            {
              if(s==1)
              {
                P1<-c(rep(0,lin),penal.vec,rep(1/Q1,n*s))
                P1<-diag(P1)
              }else{
                Q_inv<-chol2inv(chol(Q1))
                P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
                diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
                for(j in 1:n)
                  P1[(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s),(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s)]<-Q_inv
              }
            }else{
              if(all(s==1))
              {
                P1<-c(rep(0,lin),penal.vec,rep(diag(Q1)^(-1),n))
                P1<-diag(P1)
              }else{
                Q_inv<-list()
                Q_inv[[1]]<-chol2inv(chol(Q1[1:s[1],1:s[1]]))
                P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
                diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
                for(jf in 1:n[1])
                  P1[(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1]),(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1])]<-Q_inv[[1]]
                
                for (zu in 2:rnd.len)
                {
                  Q_inv[[zu]]<-chol2inv(chol(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
                  for(jf in 1:n[zu])
                    P1[(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                       (lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
                }
              }
            }
            
            if(is.null(family$multivariate)){
              score_old2<-score_vec2<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[l,]
            }else{
              score_old2<-score_vec2<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[l,]
            }  
            lambda.max<-max(abs(score_vec2[(q+1):lin]))
            
            if (BLOCK)
            {
              grad.1<-gradient.lasso.block(score.beta=score_vec2[(q+1):lin],b=Delta[l,(q+1):lin],lambda.b=lambda,block=block)
            }else{
              grad.1<-gradient.lasso(score.beta=score_vec2[(q+1):lin],b=Delta[l,(q+1):lin],lambda.b=lambda)
            }
            score_vec2[(q+1):lin] <- grad.1
            
            crit.obj<-t_change(grad=score_vec2[(q+1):lin],b=Delta[l,(q+1):lin])
            t_edge<-crit.obj$min.rate
            
            optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta[l,1:(lin+dim.smooth)],ranef=Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)],
                                               Grad=score_vec2,family=family,P=diag(P1[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]), 
                                               lower = 0, upper = Inf,K=K))
            
            t_opt<-optim.obj$par
            
            if(!NRstep && !(all(active_old==active)))
              NRstep <- T
            
            tryNR <- (t_opt<t_edge)  && NRstep  
            
            if(tryNR) 
            {
              lin_akt<-q+sum(!is.element(Delta[l,(q+1):lin],0))
              if(rnd.len==1)
              {
                if(s==1)
                {
                  P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(1/Q1,n*s)))
                }else{
                  P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
                  diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
                  for(jf in 1:n)
                    P_akt[(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s),
                          (lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s)]<-Q_inv
                }
              }else{
                if(all(s==1))
                {
                  P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(diag(Q1)^(-1),n)))
                }else{
                  P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
                  diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
                  for(jf in 1:n[1])
                    P_akt[(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1]),
                          (lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1])]<-Q_inv[[1]]
                  
                  for (zu in 2:rnd.len)
                  {
                    for(jf in 1:n[zu])
                      P_akt[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                            (lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
                  }
                }
              }
              if(is.null(family$multivariate)){
                D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
                F_gross<-t(Z_aktuell)%*%(Z_aktuell*D*SigmaInv*D)+P_akt
              }else{
                W_opt <- RcppEigenProd2(D, SigmaInv)
                F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
              }
              InvFisher3<-try(chol2inv(chol(F_gross)),silent=T)
              if(inherits(InvFisher3, "try-error"))
                InvFisher3<-try(solve(F_gross),silent=T)
              if(inherits(InvFisher3, "try-error"))
              {
                half.index<-half.index+1  
              }else{
                
                
                solve.test<-TRUE 
                Delta.test<-Delta[l,active]+nue*InvFisher3%*%score_old2[active]
                if(lin_akt>q)
                {
                  vorz<-all(sign(Delta.test[(q+1):lin_akt])==sign(betaact[(q+1):lin_akt]))  
                }else{
                  vorz<-T  
                }
              }
            }else{
              solve.test<-TRUE  
            }
            
          }  
          
          if(tryNR) 
            InvFisher<-InvFisher3
          
          score_old<-score_old2
          score_vec<-score_vec2
          Q1.very.old<-Q1.old
          Q_inv.very.old<-Q_inv.old
          Q1.old<-Q1
          Q_inv.old<-Q_inv    
          
          Eta.ma[l+1,]<-Eta
          active_old<-active
          
          
          finish<-(sqrt(sum((Eta.ma[l,]-Eta.ma[l+1,])^2))/sqrt(sum((Eta.ma[l,])^2))<control$epsilon)
          finish2<-(sqrt(sum((Eta.ma[l-1,]-Eta.ma[l+1,])^2))/sqrt(sum((Eta.ma[l-1,])^2))<control$epsilon)
          if(finish ||  finish2) #|| (all(grad.1 == 0) ))
            break
          Eta.old<-Eta
        }}
##############################################################################      

        if(rnd.len==1)
        {
          if(s==1)
          {
            P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(1/Q1.old,n*s)))
          }else{
            P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
            diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
            for(jf in 1:n)
              P_akt[(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s),
                    (lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s)]<-Q_inv.old
          }
        }else{
          if(all(s==1))
          {
            P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(diag(Q1.old)^(-1),n)))
          }else{
            P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
            diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
            for(jf in 1:n[1])
              P_akt[(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1]),
                    (lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1])]<-Q_inv.old[[1]]
            
            for (zu in 2:rnd.len)
            {
              for(jf in 1:n[zu])
                P_akt[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                      (lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.old[[zu]]
            }
          }
        }
        if(is.null(family$multivariate)){
          if(control$method!="EM")
          {
            D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
            W_opt <- D*SigmaInv*D
          }  
          F_gross <- t(Z_aktuell)%*%(Z_aktuell*W_opt)+P_akt
        }else{
          if(control$method!="EM")
            W_opt <- RcppEigenProd2(D, SigmaInv)
          W_inv_t <- RcppEigenSpChol(W_opt)
          F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
        }
        InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
        if(inherits(InvFisher2, "try-error"))
          InvFisher2<-try(solve(F_gross),silent=T)
      
      if(is.null(family$multivariate)){
        FinalHat.df<-(Z_aktuell*sqrt(W_opt))%*%InvFisher2%*%t(Z_aktuell*sqrt(W_opt))
      }else{
        FinalHat.df <- RcppEigenProd3(W_inv_t, Z_aktuell, InvFisher2)
      }
      df<-sum(diag(FinalHat.df))
      
      if(control$overdispersion)
        phi<-(sum((y-Mu)^2/family$variance(Mu)))/(N-df)
      
      ######## Final calculation
      conv.step<-l
      phi.med<-phi
      
      if(conv.step==control$steps)
      {
        cat("Warning:\n")
        cat("Algorithm did not converge!\n")
      }
      
      Delta_neu<-Delta[l,]
      Mu_opt<-Mu
      Qfinal<-Q[[l+1]]
      
      aaa<-!is.element(Delta_neu[1:(lin)],0)
      aaa <- correct.cat(aaa,c(rep(1,sum(is.na(index.new))),block))
      
      if(final.re)
      {    
        ############ final re-estimation
        
        if(s==1)
        {  
          Q.max<-max(sqrt(unlist(Q)))
          Q.min<-min(sqrt(unlist(Q)))
        }else{
          Q.max<-max(Qfinal)+1
          Q.min<-min(Qfinal)-1e-10
        }
        
        if(rnd.len==1)
        {
          glmm_fin<-try(glmm_final_smooth(y,Z_fastalles[,aaa],Phi,W,k,n,penal.vec,q_start=Qfinal,K=K,
                                          Delta_start=Delta_neu[c(aaa,rep(T,dim.smooth+n%*%s))],
                                          s,steps=control$maxIter,family=family,method=control$method.final,overdispersion=control$overdispersion,
                                          phi=phi,print.iter.final=control$print.iter.final,flushit=control$flushit,eps.final=control$eps.final,
                                          Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
          if(inherits(glmm_fin, "try-error") || glmm_fin$opt>control$maxIter-10)
          {  
            glmm_fin2<-try(glmm_final_smooth(y,Z_fastalles[,aaa],Phi,W,k,n,penal.vec,q_start=q_start,K=K,
                                             Delta_start=Delta_start[c(aaa,rep(T,dim.smooth+n%*%s))],
                                             s,steps=control$maxIter,family=family,method=control$method.final,overdispersion=control$overdispersion,
                                             phi=control$phi,print.iter.final=control$print.iter.final,flushit=control$flushit,eps.final=control$eps.final,
                                             Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
            if(!inherits(glmm_fin2, "try-error"))
            {    
              if(inherits(glmm_fin, "try-error") || (glmm_fin$opt>control$maxIter-10 && glmm_fin2$opt<control$maxIter-10))
                glmm_fin<-glmm_fin2 
            }
          }
        }else{
          glmm_fin<-try(glmm_final_multi_random_smooth(y,Z_fastalles[,aaa],Phi,W,k,penal.vec,q_start=Qfinal,K=K,
                                                       Delta_start=Delta_neu[c(aaa,rep(T,dim.smooth+n%*%s))],
                                                       s,n,steps=control$maxIter,family=family,method=control$method.final,overdispersion=control$overdispersion,
                                                       phi=phi,rnd.len=rnd.len,print.iter.final=control$print.iter.final,flushit=control$flushit,eps.final=control$eps.final,
                                                       Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
          if(inherits(glmm_fin, "try-error") || glmm_fin$opt>control$maxIter-10)
          {  
            glmm_fin2<-try(glmm_final_multi_random_smooth(y,Z_fastalles[,aaa],Phi,W,k,penal.vec,q_start=q_start,K=K,
                                                          Delta_start=Delta_start[c(aaa,rep(T,dim.smooth+n%*%s))],
                                                          s,n,steps=control$maxIter,family=family,method=control$method.final,overdispersion=control$overdispersion,
                                                          phi=control$phi,rnd.len=rnd.len,print.iter.final=control$print.iter.final,flushit=control$flushit,
                                                          eps.final=control$eps.final,Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
            if(!inherits(glmm_fin2, "try-error"))
            {    
              if(inherits(glmm_fin, "try-error") || (glmm_fin$opt>control$maxIter-10 && glmm_fin2$opt<control$maxIter-10))
                glmm_fin<-glmm_fin2 
            }
          }
        }
        
        if(inherits(glmm_fin, "try-error") || glmm_fin$opt>control$maxIter-10)
        {
          cat("\nWarning:\n")
          cat("Final Fisher scoring reestimation did not converge!\n")
        }
        
        #######
      }else{
        glmm_fin<-NA  
        class(glmm_fin)<-"try-error"
      }  
      Standard_errors<-matrix(NA,length(Delta_neu),length(Delta_neu))
      
      if(!inherits(glmm_fin, "try-error"))
      {
        EDF.matrix<-glmm_fin$EDF.matrix
        complexity.smooth<-sum(diag(EDF.matrix)[c(T,rep(F,sum(aaa)-1),rep(T,dim.smooth),rep(F,n%*%s))])  
        Delta_neu2<-Delta_neu
        Delta_neu2[c(aaa,rep(T,dim.smooth+n%*%s))]<-glmm_fin$Delta
        Standard_errors[c(aaa,rep(T,dim.smooth+n%*%s)),c(aaa,rep(T,dim.smooth+n%*%s))]<-glmm_fin$Standard_errors
        Qfinal<-glmm_fin$Q
        phi<-glmm_fin$phi
        complexity<-glmm_fin$complexity
        Delta_neu<-Delta_neu2
        Eta_opt<-Z_alles%*%Delta_neu
        if(is.null(family$multivariate)){
          Mu_opt<-family$linkinv(Eta_opt)
        }else{
          Mu_opt <- family$linkinv(Eta_opt, K)
        }
      }else{
        glmm_fin<-list()
        glmm_fin$ranef.logLik<-ranef.logLik
        complexity<-df
        
        lin_akt<-q+sum(!is.element(Delta_neu[(q+1):lin],0))
        if(rnd.len==1)
        {
          if(s==1)
          {
            P1a<-diag(c(rep(0,lin_akt+dim.smooth),rep(1/Q1,n*s)))
          }else{
            P1a<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
            for(jf in 1:n)
            {
              P1a[(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s),
                  (lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s)]<-Q_inv
            }  
          }
        }else{
          if(all(s==1))
          {
            P1a<-diag(c(rep(0,lin_akt+dim.smooth),rep(diag(Q1)^(-1),n)))
          }else{
            P1a<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
            for(jf in 1:n[1])
            {
              P1a[(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1]),
                  (lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1])]<-Q_inv[[1]]
            }
            for (zu in 2:rnd.len)
            {
              for(jf in 1:n[zu])
              {
                P1a[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                    (lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
              }}
          }
        }
        if(inherits(InvFisher2, "try-error"))
        {
          warning("No EDF's for smooth functions available, as Fisher matrix not invertible!")
          complexity.smooth<-dim.smooth
        }else{  
          ###### EDF of spline; compare Wood's Book on page 167
          if(is.null(family$multivariate)){
            D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
            W_opt <- D*SigmaInv*D
            EDF.matrix<-InvFisher2%*%(t(Z_aktuell)%*%(Z_aktuell*W_opt)+P1a)
          }else{
            W_opt <- RcppEigenProd2(D, SigmaInv)
            F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P1a
            EDF.matrix<-InvFisher2%*%F_gross
          }
          complexity.smooth<-sum(diag(EDF.matrix)[c(T,rep(F,sum(aaa)-1),rep(T,dim.smooth),rep(F,n%*%s))])
        }
      }  
      
      if(!(complexity.smooth>=1 && complexity.smooth<=dim.smooth))
        complexity.smooth<-dim.smooth
      
      if(rnd.len==1)
      {
        if(s==1)
          Qfinal<-sqrt(Qfinal)
        
        if(!is.matrix(Qfinal))
          Qfinal<-as.matrix(Qfinal)
        colnames(Qfinal)<-random.labels
        rownames(Qfinal)<-random.labels
      }else{
        Qfinal_old<-Qfinal
        Qfinal<-list()
        Qfinal[[1]]<-as.matrix(Qfinal_old[1:s[1],1:s[1]])
        colnames(Qfinal[[1]])<-random.labels[[1]]
        rownames(Qfinal[[1]])<-random.labels[[1]]
        
        for (zu in 2:rnd.len)
        {
          Qfinal[[zu]]<-as.matrix(Qfinal_old[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])
          colnames(Qfinal[[zu]])<-random.labels[[zu]]
          rownames(Qfinal[[zu]])<-random.labels[[zu]]
        }
      }
      
      colnames(Standard_errors) <- rownames(Standard_errors) <- paste0("help.",1:nrow(Standard_errors))
      
      if(dim(X)[2]>0)
      {  
        names(Delta_neu)[1:dim(X)[2]]<-colnames(X)
        colnames(Standard_errors)[1:dim(X)[2]]<-rownames(Standard_errors)[1:dim(X)[2]]<-colnames(X)
      }
      
      if(lin>1)
      {
        names(Delta_neu)[(dim(X)[2]+1):lin]<-colnames(U)
        colnames(Standard_errors)[(dim(X)[2]+1):lin]<-rownames(Standard_errors)[(dim(X)[2]+1):lin]<-colnames(U)
      }
      
      names(Delta_neu)[(lin+1):(lin+dim.smooth)]<-colnames(Phi)
      colnames(Standard_errors)[(lin+1):(lin+dim.smooth)]<-rownames(Standard_errors)[(lin+1):(lin+dim.smooth)]<-colnames(Phi)
      names(Delta_neu)[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]<-colnames(W)
      colnames(Standard_errors)[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]<-rownames(Standard_errors)[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]<-colnames(W)
      colnames(Delta)<-c(final.names,colnames(Phi),colnames(W))

      Delta_neu[1:lin] <- Delta_neu[1:lin][transf.names]
      Standard_errors[1:lin,1:lin] <-Standard_errors[1:lin,1:lin][transf.names,transf.names]
      names(Delta_neu)[1:lin] <- transf.names
      colnames(Standard_errors)[1:lin] <- rownames(Standard_errors)[1:lin] <- transf.names
      ## Transform the coefficients back to the original scale if the design
      ## matrix was standardized
      if(standardize){
        if(any.notpen)
        {
          Delta_neu[inotpen.which] <- (1 / scale.notpen) * Delta_neu[inotpen.which]
          if(!inherits(glmm_fin, "try-error"))
            Standard_errors[inotpen.which] <- (1 / scale.notpen) * Standard_errors[inotpen.which]
        }
        ## For df > 1 we have to use a matrix inversion to go back to the
        ## original scale
        for(j in 1:length(ipen.which)){
          ind <- ipen.which[[j]]
          Delta_neu[ind] <- solve(scale.pen[[j]], Delta_neu[ind,drop = FALSE])
          if(!inherits(glmm_fin, "try-error"))
          {
            Sc.help <- solve(scale.pen[[j]])
            Standard_errors[ind,ind] <- t(Sc.help)%*%(Standard_errors[ind,ind]%*%Sc.help)
          }
        }
        if(any(s>1))
          warning("Random slopes are not standardized back!")
      }
      
      Standard_errors <- sqrt(diag(Standard_errors))
      ## Need to adjust intercept if we have performed centering
      if(center){
        Delta_neu[intercept.which] <- Delta_neu[intercept.which] -
          sum(Delta_neu[1:lin][-intercept.which,drop = FALSE] * mu.x)   
      }
      
      aic<-NaN
      bic<-NaN
      
      if(is.element(family$family,c("gaussian", "binomial", "poisson","acat","cumulative"))) 
      {
        loglik<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu_opt,ranef.logLik=glmm_fin$ranef.logLik,family=family,penal=T,K=K, phi = phi)
        
        if(control$complexity!="hat.matrix")  
        {  
          if(rnd.len==1)
          {
            complexity<-0.5*(s*(s+1))
          }else{
            complexity<-0.5*(s[1]*(s[1]+1))
            for(zu in 2:rnd.len)
              complexity<-complexity+0.5*(s[zu]*(s[zu]+1))
          }
          complexity<-complexity+sum(Delta_neu[1:(lin)]!=0)+complexity.smooth
        }
        aic<--2*loglik+2*complexity
        bic<--2*loglik+log(N)*complexity
      }else{
        warning("For the specified family (so far) no AIC and BIC are available!")  
      }
      
      ret.obj=list()
      ret.obj$aic<-aic
      ret.obj$bic<-bic
      ret.obj$Deltamatrix<-Delta
      ret.obj$smooth<-Delta_neu[(lin+1):(lin+dim.smooth)]
      ret.obj$ranef<-Delta_neu[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]
      ret.obj$coefficients<-Delta_neu[1:(lin)]
      ret.obj$fixerror<-Standard_errors[1:(lin)]
      ret.obj$ranerror<-Standard_errors[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]
      ret.obj$Q_long<-Q
      ret.obj$Q<-Qfinal
      ret.obj$y_hat<-Mu_opt
      ret.obj$phi<-phi
      ret.obj$family<-family
      ret.obj$fix<-fix.old
      ret.obj$newrndfrml<-newrndfrml
      ret.obj$subject<-names(rnd)
      ret.obj$data<-data
      ret.obj$rnd.len<-rnd.len
      ret.obj$B<-B
      ret.obj$nbasis<-nbasis
      ret.obj$spline.degree<-spline.degree
      ret.obj$diff.ord<-diff.ord
      ret.obj$knots.no<-knots.no
      ret.obj$conv.step<-conv.step
      ret.obj$phi.med<-phi.med
      ret.obj$complexity.smooth<-complexity.smooth
      ret.obj$y <- y
      ret.obj$X <- cbind(X,U)
      ret.obj$W <- W
      ret.obj$K <- K
      ret.obj$df<-df
      ret.obj$loglik<-loglik
      ret.obj$lambda.max<-lambda.max
      ret.obj$logLik.vec<-logLik.vec
      return(ret.obj)
    }  
    
    #######################################################################  
    ######################## no switch to Newton Raphson ###############
    #######################################################################  
  }else{
    #######################################################################  
    ###########################  1. No Smooth #############################  
    #######################################################################  
    if(is.null(control$smooth))
    {  
      if(lin>1)
      {
        Eta_start<-X%*%beta_null+W%*%ranef_null
      }else{
        Eta_start<-rep(beta_null,N)+W%*%ranef_null
      }
      
      if(is.null(family$multivariate)){
        D<-family$mu.eta(Eta_start)
        Mu<-family$linkinv(Eta_start)
        SigmaInv <- 1/family$variance(Mu)
      }else{
        D <- family$deriv.mat(Eta_start, K)
        Mu<-family$linkinv(Eta_start, K)
        SigmaInv <- family$SigmaInv(Mu, K)
      }
      
      if(rnd.len==1)
      {
        lin0<-sum(beta_null!=0)
        if(s==1)
        {
          Q_start<-diag(q_start,s)
        }else{
          Q_start<-q_start
        }
      }else{
        lin0<-sum(beta_null!=0)
        if(all(s==1))
        {
          Q_start<-diag(diag(q_start),sum(s))
        }else{
          Q_start<-matrix(0,sum(s),sum(s))
          Q_start[1:s[1],1:s[1]]<-q_start[1:s[1],1:s[1]]
          
          for (zu in 2:rnd.len)
          {
            Q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
          }
        }
      }
      
      U<-X[,!is.na(index.new),drop=FALSE]
      X<-X[,is.na(index.new),drop=FALSE]
      
      final.names<-c(colnames(X),colnames(U))
      
      q<-dim(X)[2]
      
      Z_fastalles<-cbind(X,U)
      Z_alles<-cbind(Z_fastalles,W)
      
      ########################################################## some definitions ################################################
      Delta<-matrix(0,control$steps,(lin+n%*%s))
      Delta[1,1:lin]<-beta_null[final.names]
      Delta[1,(lin+1):(lin+n%*%s)]<-t(ranef_null)

      Eta.ma<-matrix(0,control$steps+1,N)
      Eta.ma[1,]<-Eta_start
      
      logLik.vec<-c()

      control$epsilon<-control$epsilon*sqrt(dim(Delta)[2])
      
      Delta_start<-Delta[1,]
      
      Q_inv<-NULL
      Q_inv.old.temp<-NULL
      Q_inv.start<-NULL
      
      Q<-list()
      Q[[1]]<-Q_start
      
      l=1
      
      if(rnd.len==1)
      {
        if(s==1)
        {
          P1<-c(rep(0,lin),rep((Q_start^(-1)),n*s))
          P1<-diag(P1)
        }else{
          Q_inv.start<-chol2inv(chol(Q_start))
          P1<-matrix(0,lin+n%*%s,lin+n%*%s)
          for(j in 1:n)
            P1[(lin+(j-1)*s+1):(lin+j*s),(lin+(j-1)*s+1):(lin+j*s)]<-Q_inv.start
        }
      }else{
        if(all(s==1))
        {
          P1<-c(rep(0,lin),rep(diag(Q_start)^(-1),n))
          P1<-diag(P1)
        }else{
          Q_inv.start<-list()
          Q_inv.start[[1]]<-chol2inv(chol(Q_start[1:s[1],1:s[1]]))
          P1<-matrix(0,lin+n%*%s,lin+n%*%s)
          for(jf in 1:n[1])
            P1[(lin+(jf-1)*s[1]+1):(lin+jf*s[1]),(lin+(jf-1)*s[1]+1):(lin+jf*s[1])]<-Q_inv.start[[1]]
          
          for (zu in 2:rnd.len)
          {
            Q_inv.start[[zu]]<-chol2inv(chol(Q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
            for(jf in 1:n[zu])
              P1[(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                 (lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.start[[zu]]
          }
        }
      }
      
      if(is.null(family$multivariate)){
        score_vec<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[1,]
      }else{
        score_vec<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[1,]
      }  
      lambda.max<-max(abs(score_vec[(q+1):lin]))
      
      if (BLOCK)
      {
        grad.1<-gradient.lasso.block(score.beta=score_vec[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda,block=block)
      }else{
        grad.1<-gradient.lasso(score.beta=score_vec[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda)
      }
      
      score_vec[(q+1):lin] <- grad.1
      
      crit.obj<-t_change(grad=score_vec[(q+1):lin],b=Delta[1,(q+1):lin])
      t_edge<-crit.obj$min.rate
      
      optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta_start[1:lin],ranef=Delta_start[(lin+1):(lin+n%*%s)],
                                         Grad=score_vec,family=family,P=diag(P1[(lin+1):(lin+n%*%s)]), 
                                         lower = 0, upper = Inf,K=K))
      
      t_opt<-optim.obj$par
      
      nue<-control$nue
      
      half.index<-0
      
      solve.test2<-FALSE  
      while(!solve.test2)
      {  
        
        if(half.index>50)
        {
          stop("Fisher matrix not invertible")
        }
        
        Delta[1,]<-Delta_start+min(t_opt,t_edge)*nue*(0.5^half.index)*score_vec
        if(t_opt>t_edge & half.index==0)
          Delta[1,crit.obj$whichmin+q]<-0  
        Eta<-Z_alles%*%Delta[1,]

        if(is.null(family$multivariate)){
          D<-family$mu.eta(Eta)
          Mu<-family$linkinv(Eta)
          SigmaInv <- 1/family$variance(Mu)
        }else{
          D <- family$deriv.mat(Eta, K)
          Mu<-family$linkinv(Eta, K)
          SigmaInv <- family$SigmaInv(Mu, K)
        }
        
        ranef.logLik<- -0.5*t(Delta[1,(lin+1):(lin+n%*%s)])%*%(P1[(lin+1):(lin+n%*%s),(lin+1):(lin+n%*%s)]%*%Delta[1,(lin+1):(lin+n%*%s)])
        
        logLik.vec[1]<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu,ranef.logLik=ranef.logLik,family=family,penal=T,K=K, phi = phi)
        
        active<-c(rep(T,q),!is.element(Delta[1,(q+1):lin],0),rep(T,n%*%s))
        Z_aktuell<-Z_alles[,active]
        lin_akt<-q+sum(!is.element(Delta[1,(q+1):lin],0))
        
        if (control$method=="EM")
        {  
          if(rnd.len==1)
          {
            if(s==1)
            {
              P_akt<-diag(c(rep(0,lin_akt),rep(1/Q_start,n*s)))
            }else{
              P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
              for(jf in 1:n)
                P_akt[(lin_akt+(jf-1)*s+1):(lin_akt+jf*s),(lin_akt+(jf-1)*s+1):(lin_akt+jf*s)]<-Q_inv.start
            }
          }else{
            if(all(s==1))
            {
              P_akt<-diag(c(rep(0,lin_akt),rep(diag(Q_start)^(-1),n)))
            }else{
              P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
              for(jf in 1:n[1])
                P_akt[(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1]),(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1])]<-Q_inv.start[[1]]
              
              for (zu in 2:rnd.len)
              {
                for(jf in 1:n[zu])
                  P_akt[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                        (lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.start[[zu]]
              }
            }
          }
          if(is.null(family$multivariate)){
            D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
            W_opt <- D*SigmaInv*D
            F_gross<-t(Z_aktuell)%*%(Z_aktuell*W_opt)+P_akt
          }else{
            W_opt <- RcppEigenProd2(D, SigmaInv)
            F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
          }
          InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
          if(inherits(InvFisher2, "try-error"))
            InvFisher2<-try(solve(F_gross),silent=T)
          if(inherits(InvFisher2, "try-error"))
          {
            #stop("Fisher matrix not invertible")  
            half.index<-half.index+1  
          }else{
            solve.test2<-TRUE 
          }}else{
            solve.test2<-TRUE
          }
      }
      
      betaact<-Delta[1,active]
      
      if (control$method=="EM")
      {   
        ############################# Q update ################
        if(rnd.len==1)
        {
          Q1<-InvFisher2[(lin_akt+1):(lin_akt+s),(lin_akt+1):(lin_akt+s)]+Delta[1,(lin+1):(lin+s)]%*%t(Delta[1,(lin+1):(lin+s)])
          for (i in 2:n)
            Q1<-Q1+InvFisher2[(lin_akt+(i-1)*s+1):(lin_akt+i*s),(lin_akt+(i-1)*s+1):(lin_akt+i*s)]+Delta[1,(lin+(i-1)*s+1):(lin+i*s)]%*%t(Delta[1,(lin+(i-1)*s+1):(lin+i*s)])
          Q1<-1/n*Q1
        }else{
          Q1<-matrix(0,sum(s),sum(s))
          Q1[1:s[1],1:s[1]]<-InvFisher2[(lin_akt+1):(lin_akt+s[1]),(lin_akt+1):(lin_akt+s[1])]+Delta[1,(lin+1):(lin+s[1])]%*%t(Delta[1,(lin+1):(lin+s[1])])
          for (i in 2:n[1])
            Q1[1:s[1],1:s[1]]<-Q1[1:s[1],1:s[1]]+InvFisher2[(lin_akt+(i-1)*s[1]+1):(lin_akt+i*s[1]),(lin_akt+(i-1)*s[1]+1):(lin_akt+i*s[1])]+Delta[1,(lin+(i-1)*s[1]+1):(lin+i*s[1])]%*%t(Delta[1,(lin+(i-1)*s[1]+1):(lin+i*s[1])])
          Q1[1:s[1],1:s[1]]<-1/n[1]*Q1[1:s[1],1:s[1]]
          
          for (zu in 2:rnd.len)
          {
            Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-InvFisher2[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu]),(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]+Delta[1,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]%*%t(Delta[1,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])])
            for (i in 2:n[zu])
              Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]+InvFisher2[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu]),(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]+Delta[1,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]%*%t(Delta[1,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])])
            Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-1/n[zu]*Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
          }
        }
      }else{
        if(is.null(family$multivariate)){
          Eta_tilde<-Eta+(y-Mu)*1/D
        }else{
          Eta_tilde<-Eta+solve(D)%*%(y-Mu)
        }
        Betadach<-Delta[1,1:(lin)]     
        aktuell_vec<-!is.element(Delta[1,1:(lin)],0)
        X_aktuell<-Z_fastalles[,aktuell_vec]
        
        if(rnd.len==1)
        {
          
          if(s==1)
          {
            upp<-min(20,50*Q_start)
            low<-1e-14
            optim.obj<-nlminb(sqrt(Q_start),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=Z_fastalles,X_aktuell=X_aktuell,Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp)
            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))#[(s+1):(s*(s+1)*0.5)]))
            upp<-rep(up1,length(q_start_vec))
            low<-c(rep(0,s),rep(-up1,0.5*(s^2-s)))
            optim.obj<-try(bobyqa(q_start_vec,likelihood,D=D,SigmaInv=SigmaInv,family=family,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])
            
            #### Check for positive definitness ########
            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) & !inherits(Q_solvetest, "try-error"))
                break
            }
          }   
        }else{
          if(all(s==1))
          {
            q_start_vec<-diag(q_start)
            upp<-rep(min(20,50*diag(q_start)),sum(s))
            low<-rep(1e-14,sum(s))
            optim.obj<-try(bobyqa(sqrt(q_start_vec),likelihood_diag,D=D,SigmaInv=SigmaInv,family=family,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,rnd.len=rnd.len))
            Q1<-diag(optim.obj$par)^2
          }else{
            q_start_vec<-c(diag(q_start)[1:s[1]],q_start[1:s[1],1:s[1]][lower.tri(q_start[1:s[1],1:s[1]])])
            up1<-min(20,50*max(q_start_vec))
            low<-c(rep(0,s[1]),rep(-up1,0.5*(s[1]^2-s[1])))
            
            for (zu in 2:rnd.len)
            {
              q_start_vec<-c(q_start_vec,c(diag(q_start)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])],q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]))
              up1<-min(20,50*max(q_start_vec))
              low<-c(low,c(rep(0,s[zu]),rep(-up1,0.5*(s[zu]^2-s[zu]))))
            }
            upp<-rep(up1,length(q_start_vec))
            optim.obj<-try(bobyqa(q_start_vec,likelihood_block,D=D,SigmaInv=SigmaInv,family=family,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,rnd.len=rnd.len))
            optim.vec<-optim.obj$par
            
            
            Q1<-matrix(0,sum(s),sum(s))
            diag(Q1)[1:s[1]]<-optim.vec[1:s[1]]
            if(s[1]>1)
              Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])]<-optim.vec[(s[1]+1):(s[1]*(s[1]+1)*0.5)]
            optim.vec<-optim.vec[-c(1:(s[1]*(s[1]+1)*0.5))]
            
            for (zu in 2:rnd.len)
            {
              diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-optim.vec[1:s[zu]]
              if(s[zu]>1)
                Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]<-optim.vec[(s[zu]+1):(s[zu]*(s[zu]+1)*0.5)]
              optim.vec<-optim.vec[-c(1:(s[zu]*(s[zu]+1)*0.5))]
            }
            
            #### Check for positive definitness ########
            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) & !inherits(Q_solvetest, "try-error"))
                break
            }
          }
        }}
       Q[[2]]<-Q1
      
      if(rnd.len==1)
      {
        if(s==1)
        {
          P1<-c(rep(0,lin),rep((Q1^(-1)),n*s))
          P1<-diag(P1)
        }else{
          Q_inv<-solve(Q1)
          P1<-matrix(0,lin+n%*%s,lin+n%*%s)
          for(j in 1:n)
            P1[(lin+(j-1)*s+1):(lin+j*s),(lin+(j-1)*s+1):(lin+j*s)]<-Q_inv
        }
      }else{
        if(all(s==1))
        {
          P1<-c(rep(0,lin),rep(diag(Q1)^(-1),n))
          P1<-diag(P1)
        }else{
          Q_inv<-list()
          Q_inv[[1]]<-chol2inv(chol(Q1[1:s[1],1:s[1]]))
          P1<-matrix(0,lin+n%*%s,lin+n%*%s)
          for(jf in 1:n[1])
            P1[(lin+(jf-1)*s[1]+1):(lin+jf*s[1]),(lin+(jf-1)*s[1]+1):(lin+jf*s[1])]<-Q_inv[[1]]
          
          for (zu in 2:rnd.len)
          {
            Q_inv[[zu]]<-chol2inv(chol(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
            for(jf in 1:n[zu])
              P1[(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                 (lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
          }
        }
      }
      
      if(is.null(family$multivariate)){
        score_vec2<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[1,]
      }else{
        score_vec2<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[1,]
      }  
      lambda.max<-max(abs(score_vec2[(q+1):lin]))
      
      if (BLOCK)
      {
        grad.1<-gradient.lasso.block(score.beta=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda,block=block)
      }else{
        grad.1<-gradient.lasso(score.beta=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda)
      }
      score_vec2[(q+1):lin] <- grad.1

      crit.obj<-t_change(grad=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin])
      t_edge<-crit.obj$min.rate
      
      optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta[1,1:lin],ranef=Delta[1,(lin+1):(lin+n%*%s)],
                                         Grad=score_vec2,family=family,P=diag(P1[(lin+1):(lin+n%*%s)]), 
                                         lower = 0, upper = Inf,K=K))
      
      t_opt<-optim.obj$par
      
      Eta.ma[2,]<-Eta
      
      score_vec<-score_vec2
      Q1.old<-Q1
      Q_inv.old<-Q_inv    
      
      Q1.very.old<-Q_start
      Q_inv.very.old<-Q_inv.start
      
      ###############################################################################################################################################
      ################################################################### Main Iteration ###################################################################
      if(control$steps!=1)
      {
        for (l in 2:control$steps)
        {
          if(control$print.iter)
          #  message("Iteration ",l)
          {
            cat(if(ia) "\r" else if(l > 1) "\n" else NULL)
            cat(paste("Iteration ",l))
            if(.Platform$OS.type != "unix" & ia) flush.console()
          }
          
          half.index<-0
          
          solve.test2<-FALSE  
          while(!solve.test2)
          {  
            
            if(half.index>50)
            {
              half.index<-Inf;Q1.old<-Q1.very.old;Q_inv.old<-Q_inv.very.old
            }
            
            Delta[l,]<-Delta[l-1,]+min(t_opt,t_edge)*nue*(0.5^half.index)*score_vec
            if(t_opt>t_edge & half.index==0)
              Delta[l,crit.obj$whichmin+q]<-0  
            
            Eta<-Z_alles%*%Delta[l,]
            if(is.null(family$multivariate)){
              D<-family$mu.eta(Eta)
              Mu<-family$linkinv(Eta)
              SigmaInv <- 1/family$variance(Mu)
            }else{
              D <- family$deriv.mat(Eta, K)
              Mu<-family$linkinv(Eta, K)
              SigmaInv <- family$SigmaInv(Mu, K)
            }

            ranef.logLik<- -0.5*t(Delta[l,(lin+1):(lin+n%*%s)])%*%(P1[(lin+1):(lin+n%*%s),(lin+1):(lin+n%*%s)]%*%Delta[l,(lin+1):(lin+n%*%s)])
            
            logLik.vec[l]<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu,ranef.logLik=ranef.logLik,family=family,penal=T,K=K, phi = phi)
            
            active<-c(rep(T,q),!is.element(Delta[l,(q+1):lin],0),rep(T,n%*%s))
            Z_aktuell<-Z_alles[,active]
            lin_akt<-q+sum(!is.element(Delta[l,(q+1):lin],0))
            
            if(control$method=="EM")
            {  
              if(rnd.len==1)
              {
                if(s==1)
                {
                  P_akt<-diag(c(rep(0,lin_akt),rep(1/Q1.old,n*s)))
                }else{
                  P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
                  for(jf in 1:n)
                    P_akt[(lin_akt+(jf-1)*s+1):(lin_akt+jf*s),(lin_akt+(jf-1)*s+1):(lin_akt+jf*s)]<-Q_inv.old
                }
              }else{
                if(all(s==1))
                {
                  P_akt<-diag(c(rep(0,lin_akt),rep(diag(Q1.old)^(-1),n)))
                }else{
                  P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
                  for(jf in 1:n[1])
                    P_akt[(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1]),(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1])]<-Q_inv.old[[1]]
                  
                  for (zu in 2:rnd.len)
                  {
                    for(jf in 1:n[zu])
                      P_akt[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                            (lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.old[[zu]]
                  }
                }
              }
              if(is.null(family$multivariate)){
                D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
                W_opt <- D*SigmaInv*D
                F_gross<-t(Z_aktuell)%*%(Z_aktuell*W_opt)+P_akt
              }else{
                W_opt <- RcppEigenProd2(D, SigmaInv)
                F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
              }
              InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
              if(inherits(InvFisher2, "try-error"))
                InvFisher2<-try(solve(F_gross),silent=T)
              if(inherits(InvFisher2, "try-error"))
              {
                half.index<-half.index+1  
              }else{
                solve.test2<-TRUE 
              }}else{
                solve.test2<-TRUE 
              }
          }
          
          betaact<-Delta[l,active]
          
          if (control$method=="EM")
          {        
            ############################# Q update ################
            if(rnd.len==1)
            {
              Q1<-InvFisher2[(lin_akt+1):(lin_akt+s),(lin_akt+1):(lin_akt+s)]+Delta[l,(lin+1):(lin+s)]%*%t(Delta[l,(lin+1):(lin+s)])
              for (i in 2:n)
                Q1<-Q1+InvFisher2[(lin_akt+(i-1)*s+1):(lin_akt+i*s),(lin_akt+(i-1)*s+1):(lin_akt+i*s)]+Delta[l,(lin+(i-1)*s+1):(lin+i*s)]%*%t(Delta[l,(lin+(i-1)*s+1):(lin+i*s)])
              Q1<-1/n*Q1
            }else{
              Q1<-matrix(0,sum(s),sum(s))
              Q1[1:s[1],1:s[1]]<-InvFisher2[(lin_akt+1):(lin_akt+s[1]),(lin_akt+1):(lin_akt+s[1])]+Delta[l,(lin+1):(lin+s[1])]%*%t(Delta[l,(lin+1):(lin+s[1])])
              for (i in 2:n[1])
                Q1[1:s[1],1:s[1]]<-Q1[1:s[1],1:s[1]]+InvFisher2[(lin_akt+(i-1)*s[1]+1):(lin_akt+i*s[1]),(lin_akt+(i-1)*s[1]+1):(lin_akt+i*s[1])]+Delta[l,(lin+(i-1)*s[1]+1):(lin+i*s[1])]%*%t(Delta[l,(lin+(i-1)*s[1]+1):(lin+i*s[1])])
              Q1[1:s[1],1:s[1]]<-1/n[1]*Q1[1:s[1],1:s[1]]
              
              for (zu in 2:rnd.len)
              {
                Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-InvFisher2[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu]),(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]+Delta[l,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]%*%t(Delta[l,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])])
                for (i in 2:n[zu])
                  Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]+InvFisher2[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu]),(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]+Delta[l,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]%*%t(Delta[l,(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])])
                Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-1/n[zu]*Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
              }
            }  
          }else{
            if(is.null(family$multivariate)){
              Eta_tilde<-Eta+(y-Mu)*1/D
            }else{
              Eta_tilde<-Eta+solve(D)%*%(y-Mu)
            }
            
            Betadach<-Delta[l,1:(lin)]
            
            aktuell_vec<-!is.element(Delta[l,1:(lin)],0)
            X_aktuell<-Z_fastalles[,aktuell_vec]
            
            if(rnd.len==1)
            {
              if(s==1)
              {
                if(Q1<1e-14)
                  low<-0
                
                optim.obj<-nlminb(sqrt(Q1),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=Z_fastalles,X_aktuell=X_aktuell,Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp)
                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,SigmaInv=SigmaInv,family=family,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])
                
                #### Check for positiv definitness ########
                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) & !inherits(Q_solvetest, "try-error"))
                    break
                }
              }   
            }else{
              if(all(s==1))
              {
                Q1_vec<-diag(Q1)
                optim.obj<-try(bobyqa(sqrt(Q1_vec),likelihood_diag,D=D,SigmaInv=SigmaInv,family=family,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,rnd.len=rnd.len))
                Q1<-diag(optim.obj$par)^2
              }else{
                Q1_vec<-c(diag(Q1)[1:s[1]],Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])])
                
                for (zu in 2:rnd.len)
                  Q1_vec<-c(Q1_vec,c(diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])],Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]))
                
                optim.obj<-try(bobyqa(Q1_vec,likelihood_block,D=D,SigmaInv=SigmaInv,family=family,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,rnd.len=rnd.len))
                optim.vec<-optim.obj$par
                
                Q1<-matrix(0,sum(s),sum(s))
                diag(Q1)[1:s[1]]<-optim.vec[1:s[1]]
                if(s[1]>1)
                  Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])]<-optim.vec[(s[1]+1):(s[1]*(s[1]+1)*0.5)]
                optim.vec<-optim.vec[-c(1:(s[1]*(s[1]+1)*0.5))]
                
                for (zu in 2:rnd.len)
                {
                  diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-optim.vec[1:s[zu]]
                  if(s[zu]>1)
                    Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]<-optim.vec[(s[zu]+1):(s[zu]*(s[zu]+1)*0.5)]
                  optim.vec<-optim.vec[-c(1:(s[zu]*(s[zu]+1)*0.5))]
                }
                
                #### Check for positive definitness ########
                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) & !inherits(Q_solvetest, "try-error"))
                    break
                }
              }
            }}
          
          Q[[l+1]]<-Q1
          
          if(rnd.len==1)
          {
            if(s==1)
            {
              P1<-c(rep(0,lin),rep(1/Q1,n*s))
              P1<-diag(P1)
            }else{
              Q_inv<-solve(Q1)
              P1<-matrix(0,lin+n%*%s,lin+n%*%s)
              for(j in 1:n)
                P1[(lin+(j-1)*s+1):(lin+j*s),(lin+(j-1)*s+1):(lin+j*s)]<-Q_inv
            }
          }else{
            if(all(s==1))
            {
              P1<-c(rep(0,lin),rep(diag(Q1)^(-1),n))
              P1<-diag(P1)
            }else{
              Q_inv<-list()
              Q_inv[[1]]<-chol2inv(chol(Q1[1:s[1],1:s[1]]))
              P1<-matrix(0,lin+n%*%s,lin+n%*%s)
              for(jf in 1:n[1])
                P1[(lin+(jf-1)*s[1]+1):(lin+jf*s[1]),(lin+(jf-1)*s[1]+1):(lin+jf*s[1])]<-Q_inv[[1]]
              
              for (zu in 2:rnd.len)
              {
                Q_inv[[zu]]<-chol2inv(chol(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
                for(jf in 1:n[zu])
                  P1[(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                     (lin+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
              }
            }
          }
          
          if(is.null(family$multivariate)){
            score_old2<-score_vec2<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[l,]
          }else{
            score_old2<-score_vec2<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[l,]
          }  
          lambda.max<-max(abs(score_vec2[(q+1):lin]))
          
          if (BLOCK)
          {
            grad.1<-gradient.lasso.block(score.beta=score_vec2[(q+1):lin],b=Delta[l,(q+1):lin],lambda.b=lambda,block=block)
          }else{
            grad.1<-gradient.lasso(score.beta=score_vec2[(q+1):lin],b=Delta[l,(q+1):lin],lambda.b=lambda)
          }
          
          score_vec2[(q+1):lin] <- grad.1
          
          t_edge<-crit.obj$min.rate

          optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta[l,1:lin],ranef=Delta[l,(lin+1):(lin+n%*%s)],
                                             Grad=score_vec2,family=family,P=diag(P1[(lin+1):(lin+n%*%s)]), 
                                             lower = 0, upper = Inf,K=K))
          
          t_opt<-optim.obj$par
          
          score_vec<-score_vec2
          Q1.very.old<-Q1.old
          Q_inv.very.old<-Q_inv.old
          Q1.old<-Q1
          Q_inv.old<-Q_inv    
          
          Eta.ma[l+1,]<-Eta
          finish<-(sqrt(sum((Eta.ma[l,]-Eta.ma[l+1,])^2))/sqrt(sum((Eta.ma[l,])^2))<control$epsilon)
          finish2<-(sqrt(sum((Eta.ma[l-1,]-Eta.ma[l+1,])^2))/sqrt(sum((Eta.ma[l-1,])^2))<control$epsilon)
          if(finish ||  finish2) #|| (all(grad.1 == 0) ))
            break
          Eta.old<-Eta
        }}
##############################################################      
        if(rnd.len==1)
        {
          if(s==1)
          {
            P_akt<-diag(c(rep(0,lin_akt),rep(1/Q1.old,n*s)))
          }else{
            P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
            for(jf in 1:n)
              P_akt[(lin_akt+(jf-1)*s+1):(lin_akt+jf*s),(lin_akt+(jf-1)*s+1):(lin_akt+jf*s)]<-Q_inv.old
          }
        }else{
          if(all(s==1))
          {
            P_akt<-diag(c(rep(0,lin_akt),rep(diag(Q1.old)^(-1),n)))
          }else{
            P_akt<-matrix(0,lin_akt+n%*%s,lin_akt+n%*%s)
            for(jf in 1:n[1])
              P_akt[(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1]),(lin_akt+(jf-1)*s[1]+1):(lin_akt+jf*s[1])]<-Q_inv.old[[1]]
            
            for (zu in 2:rnd.len)
            {
              for(jf in 1:n[zu])
                P_akt[(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                      (lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.old[[zu]]
            }
          }
        }
        
        if(is.null(family$multivariate)){
          if(control$method!="EM")
          {
          D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
          W_opt <- D*SigmaInv*D
          }  
          F_gross <- t(Z_aktuell)%*%(Z_aktuell*W_opt)+P_akt
        }else{
          if(control$method!="EM")
            W_opt <- RcppEigenProd2(D, SigmaInv)
          W_inv_t <- RcppEigenSpChol(W_opt)
          F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
        }
        
        InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
        if(inherits(InvFisher2, "try-error"))
          InvFisher2<-try(solve(F_gross),silent=T)
      
      
      if(is.null(family$multivariate)){
        FinalHat.df<-(Z_aktuell*sqrt(W_opt))%*%InvFisher2%*%t(Z_aktuell*sqrt(W_opt))
      }else{
        FinalHat.df <- RcppEigenProd3(W_inv_t, Z_aktuell, InvFisher2)
      }
      df<-sum(diag(FinalHat.df))
      if(abs(df)>1e+10)
      df<-Inf
        
      if(control$overdispersion)
        phi<-(sum((y-Mu)^2/family$variance(Mu)))/(N-df)
      
      conv.step<-l
      phi.med<-phi
      
      if(conv.step==control$steps)
      {
        cat("Warning:\n")
        cat("Algorithm did not converge!\n")
      }

      Delta_neu<-Delta[l,]
      Mu_opt<-Mu
      Qfinal<-Q[[l+1]]
      
      aaa<-!is.element(Delta_neu[1:(lin)],0)
      aaa <- correct.cat(aaa,c(rep(1,sum(is.na(index.new))),block))
      
      if(final.re)
      {    
        ############ final re-estimation
 
        if(rnd.len==1 && s==1)
        {  
          Q.max<-max(sqrt(unlist(Q)))
          Q.min<-min(sqrt(unlist(Q)))
        }else{
          Q.max<-max(Qfinal)+1
          Q.min<-min(Qfinal)-1e-10
        }
        
        if(rnd.len==1)
        {
          glmm_fin<-try(glmm_final(y,Z_fastalles[,aaa],W,k,n,q_start=Qfinal,K=K,
                                   Delta_start=Delta_neu[c(aaa,rep(T,n%*%s))],s,steps=control$maxIter,
                                   family=family,method=control$method.final,overdispersion=control$overdispersion,
                                   phi=phi,print.iter.final=control$print.iter.final,flushit=control$flushit,eps.final=control$eps.final,
                                   Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
          
          if(inherits(glmm_fin, "try-error") || glmm_fin$opt>control$maxIter-10)
          {  
            glmm_fin2<-try(glmm_final(y,Z_fastalles[,aaa],W,k,n,q_start=q_start,K=K,
                                      Delta_start=Delta_start[c(aaa,rep(T,n%*%s))],s,steps=control$maxIter,
                                      family=family,method=control$method.final,overdispersion=control$overdispersion,
                                      phi=control$phi,print.iter.final=control$print.iter.final,flushit=control$flushit,
                                      eps.final=control$eps.final,Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
            if(!inherits(glmm_fin2, "try-error"))
            {    
              if(inherits(glmm_fin, "try-error") || (glmm_fin$opt>control$maxIter-10 && glmm_fin2$opt<control$maxIter-10))
                glmm_fin<-glmm_fin2 
            }
          }
        }else{
          glmm_fin<-try(glmm_final_multi_random(y,Z_fastalles[,aaa],W,k,q_start=Qfinal,K=K,
                                                Delta_start=Delta_neu[c(aaa,rep(T,n%*%s))],s,n,steps=control$maxIter,
                                                family=family,method=control$method.final,overdispersion=control$overdispersion,
                                                phi=phi,rnd.len=rnd.len,print.iter.final=control$print.iter.final,flushit=control$flushit,
                                                eps.final=control$eps.final,Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
          if(inherits(glmm_fin, "try-error") || glmm_fin$opt>control$maxIter-10)
          {  
            glmm_fin2<-try(glmm_final_multi_random(y,Z_fastalles[,aaa],W,k,q_start=q_start,K=K,
                                                   Delta_start=Delta_start[c(aaa,rep(T,n%*%s))],s,n,steps=control$maxIter,
                                                   family=family,method=control$method.final,overdispersion=control$overdispersion,
                                                   phi=control$phi,rnd.len=rnd.len,print.iter.final=control$print.iter.final,flushit=control$flushit,
                                                   eps.final=control$eps.final,Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)    
            if(!inherits(glmm_fin2, "try-error"))
            {    
              if(inherits(glmm_fin, "try-error") || (glmm_fin$opt>control$maxIter-10 && glmm_fin2$opt<control$maxIter-10))
                glmm_fin<-glmm_fin2 
            }
          }
        }
        
        if(inherits(glmm_fin, "try-error") || glmm_fin$opt>control$maxIter-10)
        {
          cat("\nWarning:\n")
          cat("Final Fisher scoring reestimation did not converge!\n")
        }
        #######
        ranef.logLik<-glmm_fin$ranef.logLik
      }else{
        glmm_fin<-NA  
        class(glmm_fin)<-"try-error"
      }  
      Standard_errors<-matrix(NA,length(Delta_neu),length(Delta_neu))
      
      if(!inherits(glmm_fin, "try-error"))
      {
        Delta_neu2<-Delta_neu
        Delta_neu2[c(aaa,rep(T,n%*%s))]<-glmm_fin$Delta
        Standard_errors[c(aaa,rep(T,n%*%s)),c(aaa,rep(T,n%*%s))]<-glmm_fin$Standard_errors
        Qfinal<-glmm_fin$Q
        phi<-glmm_fin$phi
        complexity<-glmm_fin$complexity
        Delta_neu<-Delta_neu2
        Eta_opt<-Z_alles%*%Delta_neu
        if(is.null(family$multivariate)){
          Mu_opt<-family$linkinv(Eta_opt)
        }else{
          Mu_opt <- family$linkinv(Eta_opt, K)
        }
      }else{
        glmm_fin<-list()
        glmm_fin$ranef.logLik<-ranef.logLik
        complexity<-df
      }

      if(rnd.len==1)
      {
        if(s==1)
          Qfinal<-sqrt(Qfinal)
        
        if(!is.matrix(Qfinal))
          Qfinal<-as.matrix(Qfinal)
        colnames(Qfinal)<-random.labels
        rownames(Qfinal)<-random.labels
      }else{
        Qfinal_old<-Qfinal
        Qfinal<-list()
        Qfinal[[1]]<-as.matrix(Qfinal_old[1:s[1],1:s[1]])
        colnames(Qfinal[[1]])<-random.labels[[1]]
        rownames(Qfinal[[1]])<-random.labels[[1]]
        
        if(s[1]==1)
          Qfinal[[1]]<-sqrt(Qfinal[[1]])
        
        for (zu in 2:rnd.len)
        {
          Qfinal[[zu]]<-as.matrix(Qfinal_old[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])
          colnames(Qfinal[[zu]])<-random.labels[[zu]]
          rownames(Qfinal[[zu]])<-random.labels[[zu]]
          
          if(s[zu]==1)
            Qfinal[[zu]]<-sqrt(Qfinal[[zu]])
          
        }
      }

      colnames(Standard_errors) <- rownames(Standard_errors) <- paste0("help.",1:nrow(Standard_errors))
      
      if(dim(X)[2]>0)
      {  
      names(Delta_neu)[1:dim(X)[2]]<-colnames(X)
      colnames(Standard_errors)[1:dim(X)[2]]<-rownames(Standard_errors)[1:dim(X)[2]]<-colnames(X)
      }
      
      if(lin>1)
      {
        names(Delta_neu)[(dim(X)[2]+1):lin]<-colnames(U)
        colnames(Standard_errors)[(dim(X)[2]+1):lin]<-rownames(Standard_errors)[(dim(X)[2]+1):lin]<-colnames(U)
      }
      
      names(Delta_neu)[(lin+1):(lin+n%*%s)]<-colnames(W)
      colnames(Standard_errors)[(lin+1):(lin+n%*%s)]<-rownames(Standard_errors)[(lin+1):(lin+n%*%s)]<-colnames(W)
      colnames(Delta)<-c(final.names,colnames(W))
      
      Delta_neu[1:lin] <- Delta_neu[1:lin][transf.names]
      Standard_errors[1:lin,1:lin] <-Standard_errors[1:lin,1:lin][transf.names,transf.names]
      names(Delta_neu)[1:lin] <- transf.names
      colnames(Standard_errors)[1:lin] <- rownames(Standard_errors)[1:lin] <- transf.names
      
      ## Transform the coefficients back to the original scale if the design
      ## matrix was standardized
      if(standardize){
        if(any.notpen)
        {
          Delta_neu[inotpen.which] <- (1 / scale.notpen) * Delta_neu[inotpen.which]
          if(!inherits(glmm_fin, "try-error"))
            Standard_errors[inotpen.which] <- (1 / scale.notpen) * Standard_errors[inotpen.which]
        }
        ## For df > 1 we have to use a matrix inversion to go back to the
        ## original scale
        for(j in 1:length(ipen.which)){
          ind <- ipen.which[[j]]
          Delta_neu[ind] <- solve(scale.pen[[j]], Delta_neu[ind,drop = FALSE])
          if(!inherits(glmm_fin, "try-error"))
          {
            Sc.help <- solve(scale.pen[[j]])
            Standard_errors[ind,ind] <- t(Sc.help)%*%(Standard_errors[ind,ind]%*%Sc.help)
          }
        }
        if(any(s>1))
          warning("Random slopes are not standardized back!")
      }
      
      Standard_errors <- sqrt(diag(Standard_errors))
      ## Need to adjust intercept if we have performed centering
      if(center){
        Delta_neu[intercept.which] <- Delta_neu[intercept.which] -
          sum(Delta_neu[1:lin][-intercept.which,drop = FALSE] * mu.x)   
      }

      aic<-NaN
      bic<-NaN
      
      if(is.element(family$family,c("gaussian", "binomial", "poisson","acat","cumulative"))) 
      {
        
        loglik<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu_opt,ranef.logLik=ranef.logLik,family=family,penal=T,K=K, phi = phi)
        
        if(control$complexity!="hat.matrix")  
        {  
          if(rnd.len==1)
          {
            complexity<-0.5*(s*(s+1))
          }else{
            complexity<-0.5*(s[1]*(s[1]+1))
            for(zu in 2:rnd.len)
              complexity<-complexity+0.5*(s[zu]*(s[zu]+1))
          }
          complexity<-complexity+sum(Delta_neu[1:(lin)]!=0)
        }      
        aic<--2*loglik+2*complexity
        bic<--2*loglik+log(N)*complexity
      }else{
        warning("For the specified family (so far) no AIC and BIC are available!")  
      }
      
      ret.obj=list()
      ret.obj$aic<-aic
      ret.obj$bic<-bic
      ret.obj$Deltamatrix<-Delta
      ret.obj$ranef<-Delta_neu[(lin+1):(lin+n%*%s)]
      ret.obj$coefficients<-Delta_neu[1:(lin)]
      ret.obj$fixerror<-Standard_errors[1:(lin)]
      ret.obj$ranerror<-Standard_errors[(lin+1):(lin+n%*%s)]
      ret.obj$Q_long<-Q
      ret.obj$Q<-Qfinal
      ret.obj$y_hat<-Mu_opt
      ret.obj$phi<-phi
      ret.obj$family<-family
      ret.obj$fix<-fix.old
      ret.obj$conv.step<-conv.step
      ret.obj$newrndfrml<-newrndfrml
      ret.obj$subject<-subject.names
      ret.obj$data<-data
      ret.obj$rnd.len<-rnd.len
      ret.obj$phi.med<-phi.med
      ret.obj$y <- y
      ret.obj$X <- cbind(X,U)
      ret.obj$W <- W
      ret.obj$K <- K
      ret.obj$df<-df
      ret.obj$loglik<-loglik
      ret.obj$lambda.max<-lambda.max
      ret.obj$logLik.vec<-logLik.vec
      #    ret.obj$logLik.test<-logLik.test
      return(ret.obj)
      ##############################################################  
      ######################## 2. Smooth ###########################  
      ##############################################################  
    }else{

      smooth<-control$smooth
      
      if(attr(terms(smooth$formula), "intercept")==0)
      {
        variables<-attr(terms(smooth$formula),"term.labels")
        smooth$formula<- "~ +1"
        for (ir in 1:length(variables))
          smooth$formula <- paste(smooth$formula, variables[ir], sep="+")
        smooth$formula <-formula(smooth$formula)
      }
      
      B <- model.matrix(smooth$formula, data)
      B.names<-attr(B,"dimnames")[[2]]
      
      B<-as.matrix(B[,-1])
      attr(B,"dimnames")[[2]]<-B.names[-1]
      
      nbasis<-smooth$nbasis
      diff.ord<-smooth$diff.ord
      spline.degree<-smooth$spline.degree
      knots.no<-nbasis-1
      if(spline.degree<3 && (spline.degree-diff.ord)<2)
        knots.no<-knots.no+1  
      penal<-smooth$penal
      
      if(!(diff.ord<spline.degree))
        stop("Order of differences must be lower than degree of B-spline polynomials!")
      
      m<-dim(B)[2]
      
      Phi<-numeric()
      
      for (r in 1:m)
      {
        Basis<-bs.design(B[,r],diff.ord=diff.ord,spline.degree=spline.degree,knots.no=knots.no)
        Phi_temp<-cbind(Basis$X[,-1],Basis$Z)
        colnames(Phi_temp)<-paste(colnames(B)[r],rep(1:dim(Phi_temp)[2],each=1), sep=".")
        Phi<-cbind(Phi,Phi_temp)
      }
      
      dim.smooth<-dim(Phi)[2]
      
      if(is.null(smooth$start))
        smooth$start<-rep(0,dim.smooth)  
      
      smooth_null<-smooth$start
      
      if(lin>1)
      {
        Eta_start<-X%*%beta_null+Phi%*%smooth_null+W%*%ranef_null
      }else{
        Eta_start<-rep(beta_null,N)+Phi%*%smooth_null+W%*%ranef_null
      }
      
      if(is.null(family$multivariate)){
        D<-family$mu.eta(Eta_start)
        Mu<-family$linkinv(Eta_start)
        SigmaInv <- 1/family$variance(Mu)
      }else{
        D <- family$deriv.mat(Eta_start, K)
        Mu<-family$linkinv(Eta_start, K)
        SigmaInv <- family$SigmaInv(Mu, K)
      }
      
      if(rnd.len==1)
      {
        if(s==1)
        {
          Q_start<-diag(q_start,s)
        }else{
          Q_start<-q_start
        }
      }else{
        if(all(s==1))
        {
          Q_start<-diag(diag(q_start),sum(s))
        }else{
          Q_start<-matrix(0,sum(s),sum(s))
          Q_start[1:s[1],1:s[1]]<-q_start[1:s[1],1:s[1]]
          
          for (zu in 2:rnd.len)
          {
            Q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
          }
        }
      }
      
      U<-X[,!is.na(index.new),drop=FALSE]
      X<-X[,is.na(index.new),drop=FALSE]
      
      final.names<-c(colnames(X),colnames(U))
      
      q<-dim(X)[2]
      
      Z_fastalles<-cbind(X,U)
      Z_alles<-cbind(Z_fastalles,Phi,W)

      ########################################################## some definitions ################################################
      Delta<-matrix(0,control$steps,(lin+dim.smooth+n%*%s))
      Delta[1,1:lin]<-beta_null[final.names]
      Delta[1,(lin+1):(lin+dim.smooth)]<-smooth_null
      Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]<-t(ranef_null)

      control$epsilon<-control$epsilon*sqrt(dim(Delta)[2])
      
      Delta_start<-Delta[1,]
      
      active_old<-!is.element(Delta[1,],0)
      logLik.vec<-c()
      
      Eta.ma<-matrix(0,control$steps+1,N)
      Eta.ma[1,]<-Eta_start
      
      Q_inv<-NULL
      Q_inv.old.temp<-NULL
      Q_inv.start<-NULL
      
      Q<-list()
      Q[[1]]<-Q_start
      
      l=1
      
      if(diff.ord>1)
      {
        k2<-c(rep(0,diff.ord-1),rep(1,nbasis-diff.ord+1))
      }else{
        k2<-rep(1,nbasis)
      }
      k22<-rep(k2,m)
      penal.vec<-penal*k22
      
      if(rnd.len==1)
      {
        if(s==1)
        {
          P1<-c(rep(0,lin),penal.vec,rep(1/Q_start,n*s))
          P1<-diag(P1)
        }else{
          Q_inv.start<-chol2inv(chol(Q_start))
          P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
          diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
          for(j in 1:n)
            P1[(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s),(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s)]<-Q_inv
        }
      }else{
        if(all(s==1))
        {
          P1<-c(rep(0,lin),penal.vec,rep(diag(Q_start)^(-1),n))
          P1<-diag(P1)
        }else{
          Q_inv.start<-list()
          Q_inv.start[[1]]<-chol2inv(chol(Q_start[1:s[1],1:s[1]]))
          P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
          diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
          for(jf in 1:n[1])
            P1[(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1]),(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1])]<-Q_inv[[1]]
          
          for (zu in 2:rnd.len)
          {
            Q_inv.start[[zu]]<-chol2inv(chol(Q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
            for(jf in 1:n[zu])
              P1[(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                 (lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
          }
        }
      }
      
      if(is.null(family$multivariate)){
        score_vec<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[1,]
      }else{
        score_vec<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[1,]
      }  
      lambda.max<-max(abs(score_vec[(q+1):lin]))
      
      if (BLOCK)
      {
        grad.1<-gradient.lasso.block(score.beta=score_vec[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda,block=block)
      }else{
        grad.1<-gradient.lasso(score.beta=score_vec[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda)
      }
      
      score_vec[(q+1):lin] <- grad.1
      
      crit.obj<-t_change(grad=score_vec[(q+1):lin],b=Delta[1,(q+1):lin])
      t_edge<-crit.obj$min.rate
      
      optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta_start[1:(lin+dim.smooth)],ranef=Delta_start[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)],
                                         Grad=score_vec,family=family,P=diag(P1[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]), 
                                         lower = 0, upper = Inf,K=K))
      
      t_opt<-optim.obj$par
      
      nue<-control$nue
      
      half.index<-0
      
      solve.test2<-FALSE  
      while(!solve.test2)
      {  
        
        if(half.index>50)
        {
          stop("Fisher matrix not invertible")
        }
        
        Delta[1,]<-Delta_start+min(t_opt,t_edge)*nue*(0.5^half.index)*score_vec
        if(t_opt>t_edge & half.index==0)
          Delta[1,crit.obj$whichmin+q]<-0  
        Eta<-Z_alles%*%Delta[1,]
        
        if(is.null(family$multivariate)){
          D<-family$mu.eta(Eta)
          Mu<-family$linkinv(Eta)
          SigmaInv <- 1/family$variance(Mu)
        }else{
          D <- family$deriv.mat(Eta_start)
          Mu<-family$linkinv(Eta, K)
          SigmaInv <- family$SigmaInv(Mu, K)
        }          
        ranef.logLik<- -0.5*t(Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)])%*%(P1[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s),(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]%*%Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)])
        
        logLik.vec[1]<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu,ranef.logLik=ranef.logLik,family=family,penal=T,K=K, phi = phi)
        
        active<-c(rep(T,q),!is.element(Delta[1,(q+1):lin],0),rep(T,dim.smooth+n%*%s))
        Z_aktuell<-Z_alles[,active]
        lin_akt<-q+sum(!is.element(Delta[1,(q+1):lin],0))
        
        if (control$method=="EM")
        {  
          if(rnd.len==1)
          {
            if(s==1)
            {
              P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(1/Q_start,n*s)))
            }else{
              P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
              diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
              for(jf in 1:n)
                P_akt[(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s),(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s)]<-Q_inv.start
            }
          }else{
            if(all(s==1))
            {
              P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(diag(Q_start)^(-1),n)))
            }else{
              P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
              diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
              for(jf in 1:n[1])
                P_akt[(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1]),(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1])]<-Q_inv.start[[1]]
              
              for (zu in 2:rnd.len)
              {
                for(jf in 1:n[zu])
                  P_akt[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                        (lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.start[[zu]]
              }
            }
          }
          if(is.null(family$multivariate)){
            D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
            W_opt <- D*SigmaInv*D
            F_gross<-t(Z_aktuell)%*%(Z_aktuell*W_opt)+P_akt
          }else{
            W_opt <- RcppEigenProd2(D, SigmaInv)
            F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
          }
          InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
          if(inherits(InvFisher2, "try-error"))
            InvFisher2<-try(solve(F_gross),silent=T)
          if(inherits(InvFisher2, "try-error"))
          {
            #stop("Fisher matrix not invertible")  
            half.index<-half.index+1  
          }else{
            solve.test2<-TRUE 
          }}else{
            solve.test2<-TRUE
          }
      }
      
      betaact<-Delta[1,active]
      
      if (control$method=="EM")
      {   
        ############################# Q update ################
        if(rnd.len==1)
        {
          Q1<-InvFisher2[(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s),(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s)]+Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+s)]%*%t(Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+s)])
          for (i in 2:n)
            Q1<-Q1+InvFisher2[(lin_akt+dim.smooth+(i-1)*s+1):(lin_akt+dim.smooth+i*s),(lin_akt+dim.smooth+(i-1)*s+1):(lin_akt+dim.smooth+i*s)]+Delta[1,(lin+dim.smooth+(i-1)*s+1):(lin+dim.smooth+i*s)]%*%t(Delta[1,(lin+dim.smooth+(i-1)*s+1):(lin+dim.smooth+i*s)])
          Q1<-1/n*Q1
        }else{
          Q1<-matrix(0,sum(s),sum(s))
          Q1[1:s[1],1:s[1]]<-InvFisher2[(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s[1]),(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s[1])]+Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+s[1])]%*%t(Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+s[1])])
          for (i in 2:n[1])
            Q1[1:s[1],1:s[1]]<-Q1[1:s[1],1:s[1]]+InvFisher2[(lin_akt+dim.smooth+(i-1)*s[1]+1):(lin_akt+dim.smooth+i*s[1]),(lin_akt+dim.smooth+(i-1)*s[1]+1):(lin_akt+dim.smooth+i*s[1])]+Delta[1,(lin+dim.smooth+(i-1)*s[1]+1):(lin+dim.smooth+i*s[1])]%*%t(Delta[1,(lin+dim.smooth+(i-1)*s[1]+1):(lin+dim.smooth+i*s[1])])
          Q1[1:s[1],1:s[1]]<-1/n[1]*Q1[1:s[1],1:s[1]]
          
          for (zu in 2:rnd.len)
          {
            Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-InvFisher2[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu]),(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]+Delta[1,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]%*%t(Delta[1,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])])
            for (i in 2:n[zu])
              Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]+InvFisher2[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu]),(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]+Delta[1,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]%*%t(Delta[1,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])])
            Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-1/n[zu]*Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
          }
        }
      }else{
        if(is.null(family$multivariate)){
          Eta_tilde<-Eta+(y-Mu)*1/D
        }else{
          Eta_tilde<-Eta+solve(D)%*%(y-Mu)
        }
        Betadach<-Delta[1,1:(lin+dim.smooth)]     
        aktuell_vec<-!is.element(Delta[1,1:(lin)],0)
        X_aktuell<-Z_fastalles[,aktuell_vec]
        
        if(rnd.len==1)
        {
          if(s==1)
          {
            upp<-min(20,50*Q_start)
            low<-1e-14
            optim.obj<-try(nlminb(sqrt(Q_start),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),
                                  Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp))
            if(inherits(optim.obj, "try-error"))
              optim.obj<-try(bobyqa(sqrt(Q_start),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),
                                    Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp))
            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))#[(s+1):(s*(s+1)*0.5)]))
            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,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),
                                  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])
            
            #### Check for positive definitness ########
            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) & !inherits(Q_solvetest, "try-error"))
                break
            }
          }   
        }else{
          if(all(s==1))
          {
            q_start_vec<-diag(q_start)
            upp<-rep(min(20,50*diag(q_start)),sum(s))
            low<-rep(1e-14,sum(s))
            optim.obj<-try(bobyqa(sqrt(q_start_vec),likelihood_diag,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,s=s,k=k,Betadach=Betadach,W=W, lower=low,upper=upp,rnd.len=rnd.len))
            Q1<-diag(optim.obj$par)^2
          }else{
            q_start_vec<-c(diag(q_start)[1:s[1]],q_start[1:s[1],1:s[1]][lower.tri(q_start[1:s[1],1:s[1]])])
            up1<-min(20,50*max(q_start_vec))
            low<-c(rep(0,s[1]),rep(-up1,0.5*(s[1]^2-s[1])))
            
            for (zu in 2:rnd.len)
            {
              q_start_vec<-c(q_start_vec,c(diag(q_start)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])],q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(q_start[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]))
              up1<-min(20,50*max(q_start_vec))
              low<-c(low,c(rep(0,s[zu]),rep(-up1,0.5*(s[zu]^2-s[zu]))))
            }
            upp<-rep(up1,length(q_start_vec))
            optim.obj<-try(bobyqa(q_start_vec,likelihood_block,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,s=s,k=k,Betadach=Betadach,W=W, lower=low,upper=upp,rnd.len=rnd.len))
            optim.vec<-optim.obj$par
            
            Q1<-matrix(0,sum(s),sum(s))
            diag(Q1)[1:s[1]]<-optim.vec[1:s[1]]
            if(s[1]>1)
              Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])]<-optim.vec[(s[1]+1):(s[1]*(s[1]+1)*0.5)]
            optim.vec<-optim.vec[-c(1:(s[1]*(s[1]+1)*0.5))]
            
            for (zu in 2:rnd.len)
            {
              diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-optim.vec[1:s[zu]]
              if(s[zu]>1)
                Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]<-optim.vec[(s[zu]+1):(s[zu]*(s[zu]+1)*0.5)]
              optim.vec<-optim.vec[-c(1:(s[zu]*(s[zu]+1)*0.5))]
            }
            
            #### Check for positive definitness ########
            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) & !inherits(Q_solvetest, "try-error"))
                break
            }
          }
        }}
      
      Q[[2]]<-Q1
      
      if(rnd.len==1)
      {
        if(s==1)
        {
          P1<-c(rep(0,lin),penal.vec,rep(1/Q1,n*s))
          P1<-diag(P1)
        }else{
          Q_inv<-chol2inv(chol(Q1))
          Q_inv.old.temp<-Q_inv
          P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
          diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
          for(j in 1:n)
            P1[(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s),(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s)]<-Q_inv
        }
      }else{
        if(all(s==1))
        {
          P1<-c(rep(0,lin),penal.vec,rep(diag(Q1)^(-1),n))
          P1<-diag(P1)
        }else{
          Q_inv<-list()
          Q_inv[[1]]<-chol2inv(chol(Q1[1:s[1],1:s[1]]))
          P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
          diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
          for(jf in 1:n[1])
            P1[(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1]),(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1])]<-Q_inv[[1]]
          
          for (zu in 2:rnd.len)
          {
            Q_inv[[zu]]<-chol2inv(chol(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
            for(jf in 1:n[zu])
              P1[(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                 (lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
          }
        }
      }
      
      if(is.null(family$multivariate)){
        score_vec2<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[1,]
      }else{
        score_vec2<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[1,]
      }  
      lambda.max<-max(abs(score_vec2[(q+1):lin]))
      
      if (BLOCK)
      {
        grad.1<-gradient.lasso.block(score.beta=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda,block=block)
      }else{
        grad.1<-gradient.lasso(score.beta=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin],lambda.b=lambda)
      }
      score_vec2[(q+1):lin] <- grad.1

      crit.obj<-t_change(grad=score_vec2[(q+1):lin],b=Delta[1,(q+1):lin])
      t_edge<-crit.obj$min.rate

      optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta[1,1:(lin+dim.smooth)],ranef=Delta[1,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)],
                                         Grad=score_vec2,family=family,P=diag(P1[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]), 
                                         lower = 0, upper = Inf,K=K))
      
      t_opt<-optim.obj$par
      
      Eta.ma[2,]<-Eta
      
      score_vec<-score_vec2
      Q1.old<-Q1
      Q_inv.old<-Q_inv    
      
      Q1.very.old<-Q_start
      Q_inv.very.old<-Q_inv.start
      
      ###############################################################################################################################################
      ################################################################### Main Iteration ###################################################################
      if(control$steps!=1)
      {
        for (l in 2:control$steps)
        {
          if(control$print.iter)
          #  message("Iteration ",l)
          {
            cat(if(ia) "\r" else if(l > 1) "\n" else NULL)
            cat(paste("Iteration ",l))
            if(.Platform$OS.type != "unix" & ia) flush.console()
          }
          
          half.index<-0
          
          solve.test2<-FALSE  
          while(!solve.test2)
          {  
            
            if(half.index>50)
            {
              half.index<-Inf;Q1.old<-Q1.very.old;Q_inv.old<-Q_inv.very.old;
            }
            
            Delta[l,]<-Delta[l-1,]+min(t_opt,t_edge)*nue*(0.5^half.index)*score_vec
            if(t_opt>t_edge & half.index==0)
              Delta[l,crit.obj$whichmin+q]<-0
            
            Eta<-Z_alles%*%Delta[l,]
            if(is.null(family$multivariate)){
              D<-family$mu.eta(Eta)
              Mu<-family$linkinv(Eta)
              SigmaInv <- 1/family$variance(Mu)
            }else{
              D <- family$deriv.mat(Eta, K)
              Mu<-family$linkinv(Eta, K)
              SigmaInv <- family$SigmaInv(Mu, K)
            }
            
            ranef.logLik<- -0.5*t(Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)])%*%(P1[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s),(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]%*%Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)])
            
            logLik.vec[l]<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu,ranef.logLik=ranef.logLik,family=family,penal=T,K=K, phi = phi)
            
            active_old<-active
            active<-c(rep(T,q),!is.element(Delta[l,(q+1):lin],0),rep(T,dim.smooth+n%*%s))
            Z_aktuell<-Z_alles[,active]
            lin_akt<-q+sum(!is.element(Delta[l,(q+1):lin],0))
            
            if(control$method=="EM")
            {  
              if(rnd.len==1)
              {
                if(s==1)
                {
                  P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(1/Q1.old,n*s)))
                }else{
                  P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
                  diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
                  for(jf in 1:n)
                    P_akt[(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s),
                          (lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s)]<-Q_inv.old
                }
              }else{
                if(all(s==1))
                {
                  P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(diag(Q1.old)^(-1),n)))
                }else{
                  P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
                  diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
                  for(jf in 1:n[1])
                    P_akt[(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1]),
                          (lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1])]<-Q_inv.old[[1]]
                  
                  for (zu in 2:rnd.len)
                  {
                    for(jf in 1:n[zu])
                      P_akt[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                            (lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.old[[zu]]
                  }
                }
              }
              if(is.null(family$multivariate)){
                D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
                W_opt <- D*SigmaInv*D
                F_gross<-t(Z_aktuell)%*%(Z_aktuell*W_opt)+P_akt
              }else{
                W_opt <- RcppEigenProd2(D, SigmaInv)
                F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
              }
              InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
              if(inherits(InvFisher2, "try-error"))
                InvFisher2<-try(solve(F_gross),silent=T)
              if(inherits(InvFisher2, "try-error"))
              {
                half.index<-half.index+1  
              }else{
                solve.test2<-TRUE 
              }}else{
                solve.test2<-TRUE 
              }
          }
          
          betaact<-Delta[l,active]
          
          if (control$method=="EM")
          {        
            ############################# Q update ################
            if(rnd.len==1)
            {
              Q1<-InvFisher2[(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s),(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s)]+Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+s)]%*%t(Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+s)])
              for (i in 2:n)
                Q1<-Q1+InvFisher2[(lin_akt+dim.smooth+(i-1)*s+1):(lin_akt+dim.smooth+i*s),(lin_akt+dim.smooth+(i-1)*s+1):(lin_akt+dim.smooth+i*s)]+Delta[l,(lin+dim.smooth+(i-1)*s+1):(lin+dim.smooth+i*s)]%*%t(Delta[l,(lin+dim.smooth+(i-1)*s+1):(lin+dim.smooth+i*s)])
              Q1<-1/n*Q1
            }else{
              Q1<-matrix(0,sum(s),sum(s))
              Q1[1:s[1],1:s[1]]<-InvFisher2[(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s[1]),(lin_akt+dim.smooth+1):(lin_akt+dim.smooth+s[1])]+Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+s[1])]%*%t(Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+s[1])])
              for (i in 2:n[1])
                Q1[1:s[1],1:s[1]]<-Q1[1:s[1],1:s[1]]+InvFisher2[(lin_akt+dim.smooth+(i-1)*s[1]+1):(lin_akt+dim.smooth+i*s[1]),(lin_akt+dim.smooth+(i-1)*s[1]+1):(lin_akt+dim.smooth+i*s[1])]+Delta[l,(lin+dim.smooth+(i-1)*s[1]+1):(lin+dim.smooth+i*s[1])]%*%t(Delta[l,(lin+dim.smooth+(i-1)*s[1]+1):(lin+dim.smooth+i*s[1])])
              Q1[1:s[1],1:s[1]]<-1/n[1]*Q1[1:s[1],1:s[1]]
              
              for (zu in 2:rnd.len)
              {
                Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-InvFisher2[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu]),(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]+Delta[l,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])]%*%t(Delta[l,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+s[zu])])
                for (i in 2:n[zu])
                  Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]+InvFisher2[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu]),(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]+Delta[l,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])]%*%t(Delta[l,(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(i-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+i*s[zu])])
                Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-1/n[zu]*Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]
              }
            }  
          }else{
            if(is.null(family$multivariate)){
              Eta_tilde<-Eta+(y-Mu)*1/D
            }else{
              Eta_tilde<-Eta+solve(D)%*%(y-Mu)
            }
            
            Betadach<-Delta[l,1:(lin+dim.smooth)]
            
            aktuell_vec<-!is.element(Delta[l,1:(lin)],0)
            X_aktuell<-Z_fastalles[,aktuell_vec]
            
            if(rnd.len==1)
            {
              
              if(s==1)
              {
                if(Q1<1e-14)
                  low<-0
                
                optim.obj<-try(nlminb(sqrt(Q1),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp))
                if(inherits(optim.obj, "try-error"))
                  optim.obj<-bobyqa(sqrt(Q1),likelihood_nlminb,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,Betadach=Betadach,W=W, lower = low, upper = upp)
                
                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,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),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])
                
                #### Check for positiv definitness ########
                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) & !inherits(Q_solvetest, "try-error"))
                    break
                }
              }   
            }else{
              if(all(s==1))
              {
                Q1_vec<-diag(Q1)
                optim.obj<-try(bobyqa(sqrt(Q1_vec),likelihood_diag,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,s=s,k=k,Betadach=Betadach,W=W, lower=low,upper=upp,rnd.len=rnd.len))
                Q1<-diag(optim.obj$par)^2
              }else{
                Q1_vec<-c(diag(Q1)[1:s[1]],Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])])
                
                for (zu in 2:rnd.len)
                  Q1_vec<-c(Q1_vec,c(diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])],Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]))
                
                optim.obj<-try(bobyqa(Q1_vec,likelihood_block,D=D,SigmaInv=SigmaInv,family=family,X=cbind(Z_fastalles,Phi),X_aktuell=cbind(X_aktuell,Phi),Eta_tilde=Eta_tilde,n=n,s=s,k=k,Betadach=Betadach,W=W, lower=low,upper=upp,rnd.len=rnd.len))
                optim.vec<-optim.obj$par
                
                Q1<-matrix(0,sum(s),sum(s))
                diag(Q1)[1:s[1]]<-optim.vec[1:s[1]]
                if(s[1]>1)
                  Q1[1:s[1],1:s[1]][lower.tri(Q1[1:s[1],1:s[1]])]<-optim.vec[(s[1]+1):(s[1]*(s[1]+1)*0.5)]
                optim.vec<-optim.vec[-c(1:(s[1]*(s[1]+1)*0.5))]
                
                for (zu in 2:rnd.len)
                {
                  diag(Q1)[(sum(s[1:(zu-1)])+1):sum(s[1:zu])]<-optim.vec[1:s[zu]]
                  if(s[zu]>1)
                    Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])][lower.tri(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])]<-optim.vec[(s[zu]+1):(s[zu]*(s[zu]+1)*0.5)]
                  optim.vec<-optim.vec[-c(1:(s[zu]*(s[zu]+1)*0.5))]
                }
                
                #### Check for positive definitness ########
                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) & !inherits(Q_solvetest, "try-error"))
                    break
                }
              }
            }}
          
          Q[[l+1]]<-Q1
          
          if(rnd.len==1)
          {
            if(s==1)
            {
              P1<-c(rep(0,lin),penal.vec,rep(1/Q1,n*s))
              P1<-diag(P1)
            }else{
              Q_inv<-chol2inv(chol(Q1))
              P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
              diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
              for(j in 1:n)
                P1[(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s),(lin+dim.smooth+(j-1)*s+1):(lin+dim.smooth+j*s)]<-Q_inv
            }
          }else{
            if(all(s==1))
            {
              P1<-c(rep(0,lin),penal.vec,rep(diag(Q1)^(-1),n))
              P1<-diag(P1)
            }else{
              Q_inv<-list()
              Q_inv[[1]]<-chol2inv(chol(Q1[1:s[1],1:s[1]]))
              P1<-matrix(0,lin+dim.smooth+n%*%s,lin+dim.smooth+n%*%s)
              diag(P1)[(lin+1):(lin+dim.smooth)]<-penal.vec
              for(jf in 1:n[1])
                P1[(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1]),(lin+dim.smooth+(jf-1)*s[1]+1):(lin+dim.smooth+jf*s[1])]<-Q_inv[[1]]
              
              for (zu in 2:rnd.len)
              {
                Q_inv[[zu]]<-chol2inv(chol(Q1[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])]))
                for(jf in 1:n[zu])
                  P1[(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                     (lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
              }
            }
          }
          
          if(is.null(family$multivariate)){
            score_vec2<-t(Z_alles)%*%((y-Mu)*D*SigmaInv)-P1%*%Delta[l,]
          }else{
            score_vec2<-RcppEigenProd1(Z_alles, D, SigmaInv, y, Mu)-P1%*%Delta[l,]
          }  
          lambda.max<-max(abs(score_vec2[(q+1):lin]))
          
          if (BLOCK)
          {
            grad.1<-gradient.lasso.block(score.beta=score_vec2[(q+1):lin],b=Delta[l,(q+1):lin],lambda.b=lambda,block=block)
          }else{
            grad.1<-gradient.lasso(score.beta=score_vec2[(q+1):lin],b=Delta[l,(q+1):lin],lambda.b=lambda)
          }
          score_vec2[(q+1):lin] <- grad.1
          
          crit.obj<-t_change(grad=score_vec2[(q+1):lin],b=Delta[l,(q+1):lin])
          t_edge<-crit.obj$min.rate

          optim.obj<-suppressWarnings(nlminb(1e-16,taylor.opt,y=y,yhelp=yhelp,X=Z_alles,fixef=Delta[l,1:(lin+dim.smooth)],ranef=Delta[l,(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)],
                                             Grad=score_vec2,family=family,P=diag(P1[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]), 
                                             lower = 0, upper = Inf,K=K))
          
          t_opt<-optim.obj$par
          
          score_vec<-score_vec2
          Q1.very.old<-Q1.old
          Q_inv.very.old<-Q_inv.old
          Q1.old<-Q1
          Q_inv.old<-Q_inv    
          
          Eta.ma[l+1,]<-Eta
          
          finish<-(sqrt(sum((Eta.ma[l,]-Eta.ma[l+1,])^2))/sqrt(sum((Eta.ma[l,])^2))<control$epsilon)
          finish2<-(sqrt(sum((Eta.ma[l-1,]-Eta.ma[l+1,])^2))/sqrt(sum((Eta.ma[l-1,])^2))<control$epsilon)
          if(finish ||  finish2) #|| (all(grad.1 == 0) ))
            break
        }}
      
      ######## Final calculation
      if(control$method!="EM")
      {  
        if(rnd.len==1)
        {
          if(s==1)
          {
            P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(1/Q1.old,n*s)))
          }else{
            P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
            diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
            for(jf in 1:n)
              P_akt[(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s),
                    (lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s)]<-Q_inv.old
          }
        }else{
          if(all(s==1))
          {
            P_akt<-diag(c(rep(0,lin_akt),penal.vec,rep(diag(Q1.old)^(-1),n)))
          }else{
            P_akt<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
            diag(P_akt)[(lin_akt+1):(lin_akt+dim.smooth)]<-penal.vec
            for(jf in 1:n[1])
              P_akt[(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1]),
                    (lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1])]<-Q_inv.old[[1]]
            
            for (zu in 2:rnd.len)
            {
              for(jf in 1:n[zu])
                P_akt[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                      (lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv.old[[zu]]
            }
          }
        }
        if(is.null(family$multivariate)){
          if(control$method!="EM")
          {
            D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
            W_opt <- D*SigmaInv*D
          }  
          F_gross <- t(Z_aktuell)%*%(Z_aktuell*W_opt)+P_akt
        }else{
          if(control$method!="EM")
            W_opt <- RcppEigenProd2(D, SigmaInv)
          W_inv_t <- RcppEigenSpChol(W_opt)
          F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P_akt
        }
        InvFisher2<-try(chol2inv(chol(F_gross)),silent=T)
        if(inherits(InvFisher2, "try-error"))
          InvFisher2<-try(solve(F_gross),silent=T)
      }          
      
      if(is.null(family$multivariate)){
        FinalHat.df<-(Z_aktuell*sqrt(W_opt))%*%InvFisher2%*%t(Z_aktuell*sqrt(W_opt))
      }else{
        FinalHat.df <- RcppEigenProd3(W_inv_t, Z_aktuell, InvFisher2)
      }
      df<-sum(diag(FinalHat.df))
      
      if(control$overdispersion)
        phi<-(sum((y-Mu)^2/family$variance(Mu)))/(N-df)
      
      conv.step<-l
      phi.med<-phi
      
      if(conv.step==control$steps)
      {
        cat("Warning:\n")
        cat("Algorithm did not converge!\n")
      }
      
      Delta_neu<-Delta[l,]
      Eta_opt<-Z_alles%*%Delta_neu
      Mu_opt<-Mu
      Qfinal<-Q[[l+1]]
      
      aaa<-!is.element(Delta_neu[1:(lin)],0)
      aaa <- correct.cat(aaa,c(rep(1,sum(is.na(index.new))),block))
      
      if(final.re)
      {    
        ############ final re-estimation
        
        if(s==1)
        {  
          Q.max<-max(sqrt(unlist(Q)))
          Q.min<-min(sqrt(unlist(Q)))
        }else{
          Q.max<-max(Qfinal)+1
          Q.min<-min(Qfinal)-1e-10
        }
        
        if(rnd.len==1)
        {
          glmm_fin<-try(glmm_final_smooth(y,Z_fastalles[,aaa],Phi,W,k,n,penal.vec,q_start=Qfinal,K=K,
                                          Delta_start=Delta_neu[c(aaa,rep(T,dim.smooth+n%*%s))],
                                          s,steps=control$maxIter,family=family,method=control$method.final,overdispersion=control$overdispersion,
                                          phi=phi,print.iter.final=control$print.iter.final,flushit=control$flushit,eps.final=control$eps.final,
                                          Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
          if(inherits(glmm_fin, "try-error") || glmm_fin$opt>control$maxIter-10)
          {  
            glmm_fin2<-try(glmm_final_smooth(y,Z_fastalles[,aaa],Phi,W,k,n,penal.vec,q_start=q_start,K=K,
                                             Delta_start=Delta_start[c(aaa,rep(T,dim.smooth+n%*%s))],
                                             s,steps=control$maxIter,family=family,method=control$method.final,overdispersion=control$overdispersion,
                                             phi=control$phi,print.iter.final=control$print.iter.final,flushit=control$flushit,eps.final=control$eps.final,
                                             Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
            if(!inherits(glmm_fin2, "try-error"))
            {    
              if(inherits(glmm_fin, "try-error") || (glmm_fin$opt>control$maxIter-10 && glmm_fin2$opt<control$maxIter-10))
                glmm_fin<-glmm_fin2 
            }
          }
        }else{
          glmm_fin<-try(glmm_final_multi_random_smooth(y,Z_fastalles[,aaa],Phi,W,k,penal.vec,q_start=Qfinal,K=K,
                                                       Delta_start=Delta_neu[c(aaa,rep(T,dim.smooth+n%*%s))],
                                                       s,n,steps=control$maxIter,family=family,method=control$method.final,overdispersion=control$overdispersion,
                                                       phi=phi,rnd.len=rnd.len,print.iter.final=control$print.iter.final,flushit=control$flushit,eps.final=control$eps.final,
                                                       Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
          if(inherits(glmm_fin, "try-error") || glmm_fin$opt>control$maxIter-10)
          {  
            glmm_fin2<-try(glmm_final_multi_random_smooth(y,Z_fastalles[,aaa],Phi,W,k,penal.vec,q_start=q_start,K=K,
                                                          Delta_start=Delta_start[c(aaa,rep(T,dim.smooth+n%*%s))],
                                                          s,n,steps=control$maxIter,family=family,method=control$method.final,overdispersion=control$overdispersion,
                                                          phi=control$phi,rnd.len=rnd.len,print.iter.final=control$print.iter.final,flushit=control$flushit,
                                                          eps.final=control$eps.final,Q.max=Q.max,Q.min=Q.min,Q.fac=control$Q.fac),silent = TRUE)
            if(!inherits(glmm_fin2, "try-error"))
            {    
              if(inherits(glmm_fin, "try-error") || (glmm_fin$opt>control$maxIter-10 && glmm_fin2$opt<control$maxIter-10))
                glmm_fin<-glmm_fin2 
            }
          }
        }
        
        if(inherits(glmm_fin, "try-error") || glmm_fin$opt>control$maxIter-10)
        {
          cat("\nWarning:\n")
          cat("Final Fisher scoring reestimation did not converge!\n")
        }
        
        #######
      }else{
        glmm_fin<-NA  
        class(glmm_fin)<-"try-error"
      }  
      Standard_errors<-matrix(NA,length(Delta_neu),length(Delta_neu))
      
      if(!inherits(glmm_fin, "try-error"))
      {
        EDF.matrix<-glmm_fin$EDF.matrix
        complexity.smooth<-sum(diag(EDF.matrix)[c(T,rep(F,sum(aaa)-1),rep(T,dim.smooth),rep(F,n%*%s))])  
        Delta_neu2<-Delta_neu
        Delta_neu2[c(aaa,rep(T,dim.smooth+n%*%s))]<-glmm_fin$Delta
        Standard_errors[c(aaa,rep(T,dim.smooth+n%*%s)),c(aaa,rep(T,dim.smooth+n%*%s))]<-glmm_fin$Standard_errors
        Qfinal<-glmm_fin$Q
        phi<-glmm_fin$phi
        complexity<-glmm_fin$complexity
        Delta_neu<-Delta_neu2
        Eta_opt<-Z_alles%*%Delta_neu
        if(is.null(family$multivariate)){
          Mu_opt<-family$linkinv(Eta_opt)
        }else{
          Mu_opt <- family$linkinv(Eta_opt, K)
        }
      }else{
        glmm_fin<-list()
        glmm_fin$ranef.logLik<-ranef.logLik
        complexity<-df
        
        lin_akt<-q+sum(!is.element(Delta_neu[(q+1):lin],0))
        if(rnd.len==1)
        {
          if(s==1)
          {
            P1a<-diag(c(rep(0,lin_akt+dim.smooth),rep(1/Q1,n*s)))
          }else{
            P1a<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
            for(jf in 1:n)
            {
              P1a[(lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s),
                  (lin_akt+dim.smooth+(jf-1)*s+1):(lin_akt+dim.smooth+jf*s)]<-Q_inv
            }  
          }
        }else{
          if(all(s==1))
          {
            P1a<-diag(c(rep(0,lin_akt+dim.smooth),rep(diag(Q1)^(-1),n)))
          }else{
            P1a<-matrix(0,lin_akt+dim.smooth+n%*%s,lin_akt+dim.smooth+n%*%s)
            for(jf in 1:n[1])
            {
              P1a[(lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1]),
                  (lin_akt+dim.smooth+(jf-1)*s[1]+1):(lin_akt+dim.smooth+jf*s[1])]<-Q_inv[[1]]
            }
            for (zu in 2:rnd.len)
            {
              for(jf in 1:n[zu])
              {
                P1a[(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu]),
                    (lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+(jf-1)*s[zu]+1):(lin_akt+dim.smooth+n[1:(zu-1)]%*%s[1:(zu-1)]+jf*s[zu])]<-Q_inv[[zu]]
              }}
          }
        }
        if(inherits(InvFisher2, "try-error"))
        {
          warning("No EDF's for smooth functions available, as Fisher matrix not invertible!")
          complexity.smooth<-dim.smooth
        }else{  
          ###### EDF of spline; compare Wood's Book on page 167
          if(is.null(family$multivariate)){
            D <- as.vector(D);SigmaInv <- as.vector(SigmaInv)
            EDF.matrix<-InvFisher2%*%(t(Z_aktuell)%*%(Z_aktuell*D*SigmaInv*D)+P1a)
          }else{
            W_opt <- RcppEigenProd2(D, SigmaInv)
            F_gross <- t(Z_aktuell)%*%(W_opt%*%Z_aktuell)+P1a
            EDF.matrix<-InvFisher2%*%F_gross
          }
          complexity.smooth<-sum(diag(EDF.matrix)[c(T,rep(F,sum(aaa)-1),rep(T,dim.smooth),rep(F,n%*%s))])
        }
      }  
      
      if(!(complexity.smooth>=1 && complexity.smooth<=dim.smooth))
        complexity.smooth<-dim.smooth
 
      if(rnd.len==1)
      {
        if(s==1)
          Qfinal<-sqrt(Qfinal)
        
        if(!is.matrix(Qfinal))
          Qfinal<-as.matrix(Qfinal)
        colnames(Qfinal)<-random.labels
        rownames(Qfinal)<-random.labels
      }else{
        Qfinal_old<-Qfinal
        Qfinal<-list()
        Qfinal[[1]]<-as.matrix(Qfinal_old[1:s[1],1:s[1]])
        colnames(Qfinal[[1]])<-random.labels[[1]]
        rownames(Qfinal[[1]])<-random.labels[[1]]
        
        for (zu in 2:rnd.len)
        {
          Qfinal[[zu]]<-as.matrix(Qfinal_old[(sum(s[1:(zu-1)])+1):sum(s[1:zu]),(sum(s[1:(zu-1)])+1):sum(s[1:zu])])
          colnames(Qfinal[[zu]])<-random.labels[[zu]]
          rownames(Qfinal[[zu]])<-random.labels[[zu]]
        }
      }
      
      colnames(Standard_errors) <- rownames(Standard_errors) <- paste0("help.",1:nrow(Standard_errors))
      
      if(dim(X)[2]>0)
      {  
        names(Delta_neu)[1:dim(X)[2]]<-colnames(X)
        colnames(Standard_errors)[1:dim(X)[2]]<-rownames(Standard_errors)[1:dim(X)[2]]<-colnames(X)
      }
      
      if(lin>1)
      {
        names(Delta_neu)[(dim(X)[2]+1):lin]<-colnames(U)
        colnames(Standard_errors)[(dim(X)[2]+1):lin]<-rownames(Standard_errors)[(dim(X)[2]+1):lin]<-colnames(U)
      }
      
      names(Delta_neu)[(lin+1):(lin+dim.smooth)]<-colnames(Phi)
      colnames(Standard_errors)[(lin+1):(lin+dim.smooth)]<-rownames(Standard_errors)[(lin+1):(lin+dim.smooth)]<-colnames(Phi)
      names(Delta_neu)[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]<-colnames(W)
      colnames(Standard_errors)[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]<-colnames(Standard_errors)[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]<-colnames(W)
      colnames(Delta)<-c(final.names,colnames(Phi),colnames(W))
      
      Delta_neu[1:lin] <- Delta_neu[1:lin][transf.names]
      Standard_errors[1:lin,1:lin] <-Standard_errors[1:lin,1:lin][transf.names,transf.names]
      names(Delta_neu)[1:lin] <- transf.names
      colnames(Standard_errors)[1:lin] <- rownames(Standard_errors)[1:lin] <- transf.names
      ## Transform the coefficients back to the original scale if the design
      ## matrix was standardized
      if(standardize){
        if(any.notpen)
        {
          Delta_neu[inotpen.which] <- (1 / scale.notpen) * Delta_neu[inotpen.which]
          if(!inherits(glmm_fin, "try-error"))
            Standard_errors[inotpen.which] <- (1 / scale.notpen) * Standard_errors[inotpen.which]
        }
        ## For df > 1 we have to use a matrix inversion to go back to the
        ## original scale
        for(j in 1:length(ipen.which)){
          ind <- ipen.which[[j]]
          Delta_neu[ind] <- solve(scale.pen[[j]], Delta_neu[ind,drop = FALSE])
          if(!inherits(glmm_fin, "try-error"))
          {
            Sc.help <- solve(scale.pen[[j]])
            Standard_errors[ind,ind] <- t(Sc.help)%*%(Standard_errors[ind,ind]%*%Sc.help)
          }
        }
        if(any(s>1))
          warning("Random slopes are not standardized back!")
      }
      
      Standard_errors <- sqrt(diag(Standard_errors))
      ## Need to adjust intercept if we have performed centering
      if(center){
        Delta_neu[intercept.which] <- Delta_neu[intercept.which] -
          sum(Delta_neu[1:lin][-intercept.which,drop = FALSE] * mu.x)   
      }

      aic<-NaN
      bic<-NaN
      
      if(is.element(family$family,c("gaussian", "binomial", "poisson","acat","cumulative"))) 
      {
        loglik<-logLikeli.glmmLasso(y=y,yhelp=yhelp,mu=Mu_opt,ranef.logLik=glmm_fin$ranef.logLik,family=family,penal=T,K=K, phi = phi)
        
        if(control$complexity!="hat.matrix")  
        {  
          if(rnd.len==1)
          {
            complexity<-0.5*(s*(s+1))
          }else{
            complexity<-0.5*(s[1]*(s[1]+1))
            for(zu in 2:rnd.len)
              complexity<-complexity+0.5*(s[zu]*(s[zu]+1))
          }
          complexity<-complexity+sum(Delta_neu[1:(lin)]!=0)+complexity.smooth
        }      
        aic<--2*loglik+2*complexity
        bic<--2*loglik+log(N)*complexity
      }else{
        warning("For the specified family (so far) no AIC and BIC are available!")  
      }
      
      ret.obj=list()
      ret.obj$aic<-aic
      ret.obj$bic<-bic
      ret.obj$Deltamatrix<-Delta
      ret.obj$smooth<-Delta_neu[(lin+1):(lin+dim.smooth)]
      ret.obj$ranef<-Delta_neu[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]
      ret.obj$coefficients<-Delta_neu[1:(lin)]
      ret.obj$fixerror<-Standard_errors[1:(lin)]
      ret.obj$ranerror<-Standard_errors[(lin+dim.smooth+1):(lin+dim.smooth+n%*%s)]
      ret.obj$Q_long<-Q
      ret.obj$Q<-Qfinal
      ret.obj$y_hat<-Mu_opt
      ret.obj$phi<-phi
      ret.obj$family<-family
      ret.obj$fix<-fix.old
      ret.obj$newrndfrml<-newrndfrml
      ret.obj$subject<-names(rnd)
      ret.obj$data<-data
      ret.obj$rnd.len<-rnd.len
      ret.obj$B<-B
      ret.obj$nbasis<-nbasis
      ret.obj$spline.degree<-spline.degree
      ret.obj$diff.ord<-diff.ord
      ret.obj$knots.no<-knots.no
      ret.obj$conv.step<-conv.step
      ret.obj$phi.med<-phi.med
      ret.obj$complexity.smooth<-complexity.smooth
      ret.obj$y <- y
      ret.obj$X <- cbind(X,U)
      ret.obj$W <- W
      ret.obj$K <- K
      ret.obj$df<-df
      ret.obj$loglik<-loglik
      ret.obj$lambda.max<-lambda.max
      ret.obj$logLik.vec<-logLik.vec
      ret.obj$rnd <- rnd
      return(ret.obj)
    }  
  }
  ##################################################################
  ##################################################################
}

Try the glmmLasso package in your browser

Any scripts or data that you put into this service are public.

glmmLasso documentation built on Aug. 23, 2023, 5:06 p.m.