R/resf_vc.R

Defines functions print.resf_vc resf_vc

Documented in resf_vc

resf_vc<- function( y, x, xconst = NULL, xgroup = NULL, weight = NULL, offset = NULL,
                    x_nvc = FALSE, xconst_nvc = FALSE, x_sel = TRUE,
                    x_nvc_sel = TRUE, xconst_nvc_sel = TRUE, nvc_num = 5, meig,
                    method = "reml", penalty = "bic", nongauss = NULL,
                    miniter=NULL, maxiter = 30, tol = 1e-30){

  erat_min <- 0.001
  {
    af <-function(par,y) par[1]+par[2]*y
    sc <-function(par,y) (y-par[1])/par[2]
    sa <-function(par,y) sinh(par[1]*asinh(y)-par[2])
    bc <-function(par,y, jackup){
      yy     <- y + ifelse(jackup, abs( par[2] ), 0 )
      if(par[1]==0){
        res <- log(yy)
      } else {
        res <- ( yy^par[1] - 1 ) / par[1]
      }
    }
    sal<-function(par,y, noconst=FALSE){
      if(noconst==FALSE){
        par[3]+ par[4]*sinh(par[1]*asinh(y)-par[2])
      } else {
        sinh(par[1]*asinh(y)-par[2])
      }
    }

    d_af <-function(par,y) par[2]
    d_sc <-function(par,y) rep(1/par[2], length(y))###
    d_sa <-function(par,y) par[1]*cosh(par[1]*asinh(y)-par[2])/sqrt(1+y^2)
    d_bc <-function(par,y,jackup) (y + ifelse(jackup, abs( par[2] ), 0 ) )^(par[1]-1)#abs(y)^(par-1)
    d_sal<-function(par,y,noconst=FALSE){
      if(noconst==FALSE){
        par[4]*d_sa(par=par[1:2],y)#d_af(par=par[3:4],sa(par=par[1:2],y))
      } else {
        d_sa(par=par[1:2],y)# d_af(par=c(0, 1),sa(par=par[1:2],y))*
      }
    }

    i_af<-function(par,y) (y-par[1])/par[2]
    i_sc <-function(par,y) par[1]+par[2]*y
    i_sa<-function(par,y) sinh(1/par[1]*(asinh(y)+par[2]))
    i_bc<-function(par,y,jackup){#, ulim=NULL
      y_pre <- y
      if(par[1]==0){
        y   <-exp(y) - ifelse(jackup, abs( par[2] ), 0 )
      } else {
        y   <- (par[1]*y + 1)^(1/par[1]) - ifelse(jackup, abs( par[2] ), 0 )
      }

      #if( !is.null(ulim) ){
      #  ulim_id <- which(y_pre > ulim)
      #  if( length( ulim_id > 1) ){
      #    y_lower    <-sort(unique(y[-ulim_id]),decreasing=TRUE)[2:1]
      #    y_pre_lower<-sort(unique(y_pre[-ulim_id]),decreasing=TRUE)[2:1]
      #    grad_lower <-diff(y_lower)/diff(y_pre_lower)
      #    y[ulim_id] <-y_lower[2]+(y_pre[ulim_id] - y_pre_lower[2])*grad_lower
      #  }
      #}

      if( sum( is.na(y) )>0 ) y[is.na(y)]<-min(y[!is.na(y)]/2)
      return(y)
    }

    i_sal<-function(par,y,noconst=FALSE){
      if(noconst==FALSE){
        i_sa(par=par[1:2],i_af(par=par[3:4],y))
      } else {
        i_sa(par=par[1:2],i_af(par=c(0, 1),y))
      }
    }

    d_sal_k   <-function(par,y,k=2,noconst_last=TRUE,bc_par=NULL,y_ms=NULL,z_ms,jackup){
      d       <-1
      if( !is.null( bc_par[1] ) ){
        d     <- d_bc(bc_par,y,jackup=jackup)*d
        y     <-   bc(bc_par,y,jackup=jackup)
      }

      if(k>0){
        d       <- d_sc(par=y_ms,y)*d
        y       <-   sc(par=y_ms,y)

        for(kk in 1:k){
          nc_dum<-ifelse((noconst_last==TRUE)&(kk==k),TRUE,FALSE)
          d     <-d_sal(par[[kk]],y,noconst=nc_dum)*d
          y     <-  sal(par[[kk]],y,noconst=nc_dum)
        }
      }
      d       <-d_sc(par=z_ms,y)*d
      return(d)
    }
    i_sal_k   <-function(par,y,k=2,noconst_last=TRUE,bc_par=NULL,y_ms=NULL,z_ms, jackup,
                         y_added2=0){#,ulim=NULL
      y_pre   <-y
      y       <-i_sc(par=z_ms,y)
      y       <-y - y_added2
      if(k>0){
        for(kk in k:1) {
          nc_dum<-ifelse((noconst_last==TRUE)&(kk==k),TRUE,FALSE)
          y <- i_sal(par[[kk]],y,noconst=nc_dum)
        }
      }

      if( !is.null( y_ms ) ){
        y     <-i_sc(par=y_ms,y)
      }

      if( !is.null(bc_par[1]) ){
        y <- i_bc(par=bc_par, y, jackup=jackup)
      }

      #if( !is.null(ulim) ){
      #  ulim_id <- which(y_pre > ulim)
      #  if( length( ulim_id > 1) ){
      #    y_ulim     <- y[ulim_id]
      #    y_pre_ulim <-y_pre[ulim_id]
      #    y_lower    <-sort(y[-ulim_id],decreasing=TRUE)[2:1]
      #    y_pre_lower<-sort(y_pre[-ulim_id],decreasing=TRUE)[2:1]
      #    grad_lower <-diff(y_lower)/diff(y_pre_lower)
      #    y[ulim_id] <-y_lower[2]+(y_pre_ulim - y_pre_lower[2])*grad_lower
      #  }
      #}

      if( sum( is.na(y) )>0 ) y[is.na(y)]<-min(y[!is.na(y)]/2)
      return(y)
    }
    sal_k     <-function(par,y,k=2,noconst_last=TRUE,bc_par=NULL, jackup,y_added2=0){

      if( !is.null(bc_par[1]) ){
        y   <- bc(par=bc_par,y,jackup=jackup)
      }

      if( k > 0){
        y_ms<- c( mean(y), sd(y) )
        y   <- sc(par=y_ms,y)

        for(kk in 1:k) {
          nc_dum<-ifelse((noconst_last==TRUE)&(kk==k),TRUE,FALSE)
          y <- sal(par[[kk]],y,noconst=nc_dum)
        }
      } else {
        y_ms<- NULL
      }
      y     <- y + y_added2

      z_ms  <- c( mean(y), sd(y) )
      y     <- sc(par=z_ms,y)

      return(list(y=y,y_ms=y_ms, z_ms=z_ms))
    }

    ######## Negative log-likelihood (SAL distribution)
    NLL_sal<-function(par,y,M,Minv,m0,k=2,noconst_last=TRUE,y_nonneg=FALSE,jackup,weight=NULL,plim,
                      bc_value=NULL,y_added2=0, y_type){
      n   <-length(y)
      par[par > plim]  <-  plim
      par[par < -plim] <- -plim

      par2<-list(NULL)
      for(kk in 1:k){
        if(noconst_last & (kk==k)){
          par[(4*(kk-1)+1)]<- abs( par[(4*(kk-1)+1)] )
          par2[[kk]]        <- par[(4*(kk-1)+1):(4*(kk-1)+2)]

        } else {
          par[(4*(kk-1)+1)]<- abs( par[(4*(kk-1)+1)] )
          par[(4*(kk-1)+4)]<- abs( par[(4*(kk-1)+4)] )
          par2[[kk]]        <- par[(4*(kk-1)+1):(4*(kk-1)+4)]
        }
      }

      if(y_nonneg){
        np_b  <-length(par)
        bc_par<-par[(np_b-1):np_b]
        bc_par[2]<-abs(bc_par[2])
        if( !is.null( bc_value ) ){
          bc_par[1]<- ifelse( bc_value < -5, -5, ifelse( bc_value > 5, 5, bc_value) )
        } else {
          bc_par[1]<- ifelse( bc_par[1] < -5, -5, ifelse( bc_par[1] > 5, 5, bc_par[1]) )
        }

      } else {
        bc_par<-NULL
      }

      z0  <- sal_k(par=par2,y=y,k=k,noconst_last=noconst_last,bc_par=bc_par,jackup=jackup,
                   y_added2=y_added2)
      z   <- z0$y
      z_ms<- z0$z_ms
      y_ms<- z0$y_ms

      #if(1<0){
      steepen <-0
      if( y_type=="count" ){
        yz0    <-cbind(y,z)

        yz_ord <-yz0[order(yz0[,1]),]
        yz     <-yz_ord[!duplicated(yz_ord[,1]),]
        uni_y  <-yz[,1];duni_y<-diff(uni_y)
        uni_z  <-yz[,2];duni_z<-diff(uni_z)
        duni_n <-length(duni_y)
        steep1a<-max( (duni_z/duni_y)[-duni_n]/(duni_z/duni_y)[-1] )
        steep1b<-max( (duni_z/duni_y)[-1]/(duni_z/duni_y)[-duni_n] )

        yz_ord <-yz0[order(yz0[,2]),]
        yz     <-yz_ord[!duplicated(yz_ord[,2]),]
        uni_y  <-yz[,1];duni_y<-diff(uni_y)
        uni_z  <-yz[,2];duni_z<-diff(uni_z)
        duni_n <-length(duni_y)
        steep1c<-max( (duni_y/duni_z)[-duni_n]/(duni_y/duni_z)[-1] )
        steep1d<-max( (duni_y/duni_z)[-1]/(duni_y/duni_z)[-duni_n] )

        steep1 <-max( c(steep1a, steep1b) )#steep1c, steep1d

        if( length(duni_y) >=4 ){
          steep2   <- (duni_z[2]/duni_y[2]) / (duni_z[3]/duni_y[3])
          steep3   <- (duni_z[3]/duni_y[3]) / (duni_z[4]/duni_y[4])
          dsteep123<- c(steep2/steep1, steep3/steep2)
          dsteepen <- max(dsteep123)/min(dsteep123)
          dsteep   <- (steep2 > steep1 & steep3 < steep2)|(steep2 < steep1 & steep3 > steep2)
        } else {
          dsteep   <- FALSE
        }

        duni_y_ex_last<-diff( uni_z[(duni_n):(duni_n+1)] )
        duni_y_last   <-duni_y[duni_n]
        dsteep_last   <-duni_y_last/duni_y_ex_last

        steepen<-ifelse( steep1 < 10,0,10^10 + steep1^3) +
          ifelse( dsteep==FALSE,0,10^10 + dsteepen^3) +
          ifelse( dsteep_last < 1.5, 0, 10^10 + dsteep_last^2)
      }
      #}

      m   <- t(m0)%*%z
      b   <- Minv%*%m
      if( !is.null(weight)){
        ee	<- sum( weight*z ^ 2 ) - 2 * t( b ) %*% m + t( b ) %*% M %*% b
      } else {
        ee	<- sum( z ^ 2 ) - 2 * t( b ) %*% m + t( b ) %*% M %*% b
      }
      d_abs<-abs( d_sal_k(par=par2,y,k=k,noconst_last=noconst_last,bc_par=bc_par,y_ms=y_ms,z_ms=z_ms,jackup=jackup) )
      comp<- sum( log( d_abs ) )
      nll <- n * ( 1 + log( 2 * pi * ee / n ) ) -2*comp + steepen
    }

    NLL_sal2<-function(par,y,M,Minv,m0,k=2,noconst_last=TRUE,y_nonneg=FALSE,jackup,weight=NULL,
                       plim, bc_value=NULL,y_added2=0){
      par[par > plim]  <-  plim
      par[par < -plim] <- -plim

      n   <-length(y)
      par2<-list(NULL)
      for(kk in 1:k){
        if(noconst_last & (kk==k)){
          par[(4*(kk-1)+1)]<- abs( par[(4*(kk-1)+1)] )
          par2[[kk]]        <- par[(4*(kk-1)+1):(4*(kk-1)+2)]
        } else {
          par[(4*(kk-1)+1)]<- abs( par[(4*(kk-1)+1)] )
          par[(4*(kk-1)+4)]<- abs( par[(4*(kk-1)+4)] )
          par2[[kk]]        <- par[(4*(kk-1)+1):(4*(kk-1)+4)]
        }
      }

      if(y_nonneg){
        np_b    <-length(par)
        bc_par<-par[(np_b-1):np_b]
        bc_par[2]<-abs(bc_par[2])
        if( !is.null( bc_value ) ){
          bc_par[1]<- ifelse( bc_value < -5, -5, ifelse( bc_value > 5, 5, bc_value) )
        } else {
          bc_par[1]<- ifelse( bc_par[1] < -5, -5, ifelse( bc_par[1] > 5, 5, bc_par[1]) )
        }

      } else {
        bc_par<-NULL
      }

      z0  <- sal_k(par=par2,y=y,k=k,noconst_last=noconst_last,bc_par=bc_par,jackup=jackup,
                   y_added2=y_added2)
      z   <- z0$y
      z_ms<- z0$z_ms
      y_ms<- z0$y_ms

      m   <- t(m0)%*%z
      b   <- Minv%*%m
      if( !is.null(weight)){
        ee	<- sum( weight*z ^ 2 ) - 2 * t( b ) %*% m + t( b ) %*% M %*% b
      } else {
        ee	<- sum( z ^ 2 ) - 2 * t( b ) %*% m + t( b ) %*% M %*% b
      }
      d_abs <- abs( d_sal_k(par=par2,y,k=k,noconst_last=noconst_last,bc_par=bc_par,y_ms=y_ms,z_ms=z_ms,jackup=jackup) )
      comp  <- sum( log( d_abs ) )
      #loglik<- ee/2 - comp
      nll   <- n * ( 1 + log( 2 * pi * ee / n ) ) + (-2)*comp
      loglik<- (-1/2)*nll
      return(list(z=z, b=b, loglik=loglik,comp=comp,y_ms=y_ms,z_ms=z_ms))
    }

    ######## Negative log-likelihood (BC distribution)
    NLL_bc <-function(par,y,M,Minv,m0,jackup,weight=NULL,bc_value=NULL,y_added2=0){
      n     <-length(y)
      par[2]<-abs(par[2])
      if( !is.null( bc_value ) ){
        par[1]<- ifelse( bc_value < -5, -5, ifelse( bc_value > 5, 5, bc_value) )
      } else {
        par[1]<- ifelse( par[1] < -5, -5, ifelse( par[1] > 5, 5, par[1]) )
      }
      z0  <- bc(par=par,y=y,jackup=jackup) + y_added2
      z_ms<- c(mean(z0),sd(z0))
      z   <- (z0-z_ms[1])/z_ms[2]

      #if(1<0){
      steepen <-0
      if( y_type=="count" ){
        yz0   <-cbind(y,z)

        yz_ord<-yz0[order(yz0[,1]),]
        yz    <-yz_ord[!duplicated(yz_ord[,1]),]
        uni_y <-yz[,1];duni_y<-diff(uni_y)
        uni_z <-yz[,2];duni_z<-diff(uni_z)
        duni_n<-length(duni_y)
        steep1a<-max( (duni_z/duni_y)[-duni_n]/(duni_z/duni_y)[-1] )
        steep1b<-max( (duni_z/duni_y)[-1]/(duni_z/duni_y)[-duni_n] )

        yz_ord<-yz0[order(yz0[,2]),]
        yz    <-yz_ord[!duplicated(yz_ord[,2]),]
        uni_y <-yz[,1];duni_y<-diff(uni_y)
        uni_z <-yz[,2];duni_z<-diff(uni_z)
        duni_n<-length(duni_y)
        steep1c<-max( (duni_y/duni_z)[-duni_n]/(duni_y/duni_z)[-1] )
        steep1d<-max( (duni_y/duni_z)[-1]/(duni_y/duni_z)[-duni_n] )

        steep1 <-max( c(steep1a, steep1b) )#,steep1c, steep1d

        if( length(duni_y) >=4 ){
          steep2   <- (duni_z[2]/duni_y[2]) / (duni_z[3]/duni_y[3])
          steep3   <- (duni_z[3]/duni_y[3]) / (duni_z[4]/duni_y[4])
          dsteep123<- c(steep2/steep1, steep3/steep2)
          dsteepen <- max(dsteep123)/min(dsteep123)
          dsteep   <- (steep2 > steep1 & steep3 < steep2)|(steep2 < steep1 & steep3 > steep2)
        } else {
          dsteep   <- FALSE
        }

        duni_y_ex_last<-diff( uni_z[(duni_n):(duni_n+1)] )
        duni_y_last   <-duni_y[duni_n]
        dsteep_last   <-duni_y_last/duni_y_ex_last

        steepen<-ifelse( steep1 < 10,0,10^10 + steep1^3) +
          ifelse( dsteep==FALSE,0,10^10 + dsteepen^3) +
          ifelse( dsteep_last < 1.5, 0, 10^10 + dsteep_last^2)
      }
      #}

      m   <- t(m0)%*%z
      b   <- Minv%*%m
      if( !is.null(weight)){
        ee	<- sum( weight*z ^ 2 ) - 2 * t( b ) %*% m + t( b ) %*% M %*% b
      } else {
        ee	<- sum( z ^ 2 ) - 2 * t( b ) %*% m + t( b ) %*% M %*% b
      }
      d_abs<-abs( d_bc( par=par,y,jackup=jackup )*d_sc( par=z_ms,z0 ) )# y -> z0 ####### Fix this ############################
      comp<- sum( log( d_abs ) )
      nll   <- n * ( 1 + log( 2 * pi * ee / n ) )  -2 * comp + steepen
    }

    NLL_bc2<-function(par,y, M, Minv,m0,jackup,weight=NULL,bc_value=NULL,y_added2=0){
      n     <-length(y)
      par[2]<-abs(par[2])
      if( !is.null( bc_value ) ){
        par[1]<- ifelse( bc_value < -5, -5, ifelse( bc_value > 5, 5, bc_value) )
      } else {
        par[1]<- ifelse( par[1] < -5, -5, ifelse( par[1] > 5, 5, par[1]) )
      }

      z0  <- bc(par=par,y=y,jackup=jackup) + y_added2
      z_ms<- c(mean(z0),sd(z0))
      z   <- (z0-z_ms[1])/z_ms[2]#z0#
      m   <- t(m0)%*%z
      b   <- Minv%*%m
      if( !is.null(weight)){
        ee	<- sum( weight*z ^ 2 ) - 2 * t( b ) %*% m + t( b ) %*% M %*% b
      } else {
        ee	<- sum( z ^ 2 ) - 2 * t( b ) %*% m + t( b ) %*% M %*% b
      }
      d_abs <-abs( d_bc( par=par,y,jackup=jackup )*d_sc( par=z_ms,z0 ) )# y->z0 ###### Fix this ############################
      comp  <- sum( log( d_abs ) )
      nll   <- n * ( 1 + log( 2 * pi * ee / n ) ) -2*comp
      loglik<- (-1/2)*nll
      return(list(z=z, b=b, loglik=loglik,comp=comp,z_ms=z_ms))#
    }

    lik_resf_vc_tr <- function( par, k, y_nonneg=FALSE,noconst_last=TRUE, X, weight=NULL, bc_value=NULL,
                                evSqrt, y0, n, nx, emet,M0,Minv,term1,null_dum3=NULL,jackup, plim,
                                y_added2,y_type){#,omit_id
      tr_par                 <-  par
      tr_par[tr_par > plim]  <-  plim
      tr_par[tr_par < -plim] <- -plim

      steepen <-0
      if(k > 0){
        tr_par2<-list(NULL)
        for(kk in 1:k){
          if(noconst_last & (kk==k)){
            tr_par2[[kk]]        <- tr_par[(4*(kk-1)+1):(4*(kk-1)+2)]
            tr_par2[[kk]][1]     <- abs(tr_par2[[kk]][1])
          } else {
            tr_par2[[kk]]        <- tr_par[(4*(kk-1)+1):(4*(kk-1)+4)]
            tr_par2[[kk]][c(1,4)]<- abs(tr_par2[[kk]][c(1,4)])
          }
        }

        if(y_nonneg){
          np_b  <-length(tr_par)
          bc_par<-tr_par[(np_b-1):np_b]
          #if(bc_par[1] < -0.5) bc_par[1] <- -0.5###################### added 2020/12/14
          bc_par[2]<-abs(bc_par[2])
          if( !is.null( bc_value ) ){
            bc_par[1]<- ifelse( bc_value < -5, -5, ifelse( bc_value > 5, 5, bc_value) )
          } else {
            bc_par[1]<- ifelse( bc_par[1] < -5, -5, ifelse( bc_par[1] > 5, 5, bc_par[1]) )
          }

        } else {
          bc_par<-NULL
        }
        z0      <- sal_k(par=tr_par2,y=y0,k=k,noconst_last=noconst_last,bc_par=bc_par,
                         jackup=jackup,y_added2 = y_added2)
        z       <- z0$y
        z_ms    <- z0$z_ms
        y_ms    <- z0$y_ms

        #if(1<0){
        steepen <-0
        if( y_type=="count" ){
          yz0   <-cbind(y0,z)

          yz_ord<-yz0[order(yz0[,1]),]
          yz    <-yz_ord[!duplicated(yz_ord[,1]),]
          uni_y <-yz[,1];duni_y<-diff(uni_y)
          uni_z <-yz[,2];duni_z<-diff(uni_z)
          duni_n<-length(duni_y)
          steep1a<-max( (duni_z/duni_y)[-duni_n]/(duni_z/duni_y)[-1] )
          steep1b<-max( (duni_z/duni_y)[-1]/(duni_z/duni_y)[-duni_n] )

          yz_ord<-yz0[order(yz0[,2]),]
          yz    <-yz_ord[!duplicated(yz_ord[,2]),]
          uni_y <-yz[,1];duni_y<-diff(uni_y)
          uni_z <-yz[,2];duni_z<-diff(uni_z)
          duni_n<-length(duni_y)
          steep1c<-max( (duni_y/duni_z)[-duni_n]/(duni_y/duni_z)[-1] )
          steep1d<-max( (duni_y/duni_z)[-1]/(duni_y/duni_z)[-duni_n] )
          steep1 <-max( c(steep1a, steep1b) )#,steep1c, steep1d

          if( length(duni_y) >=4 ){
            steep2   <- (duni_z[2]/duni_y[2]) / (duni_z[3]/duni_y[3])
            steep3   <- (duni_z[3]/duni_y[3]) / (duni_z[4]/duni_y[4])
            dsteep123<- c(steep2/steep1, steep3/steep2)
            dsteepen <- max(dsteep123)/min(dsteep123)
            dsteep   <- (steep2 > steep1 & steep3 < steep2)|(steep2 < steep1 & steep3 > steep2)
          } else {
            dsteep   <- FALSE
          }

          duni_y_ex_last<-diff( uni_z[(duni_n):(duni_n+1)] )
          duni_y_last   <-duni_y[duni_n]
          dsteep_last   <-duni_y_last/duni_y_ex_last

          steepen<-ifelse( steep1 < 10,0,10^10 + steep1^3) +
            ifelse( dsteep==FALSE,0,10^10 + dsteepen^3) +
            ifelse( dsteep_last < 1.5, 0, 10^10 + dsteep_last^2)
        }
        #}

        d_abs   <- abs( d_sal_k( par=tr_par2 ,y=y0, k=k,noconst_last=noconst_last,
                                 bc_par=bc_par,y_ms=y_ms,z_ms=z_ms,jackup=jackup ) )
        comp    <- (-2)*sum( log( d_abs ))

        if( is.na(max(abs(z_ms))) ){
          comp  <- 10^100
        } else if(max(abs(z_ms)) > 10^10){######### deleted
          comp  <- 10^100
        }

      } else if(y_nonneg){
        if( !is.null( bc_value ) ){
          tr_par[1]<- ifelse( bc_value < -5, -5, ifelse( bc_value > 5, 5, bc_value) )
        } else {
          tr_par[1]<- ifelse( tr_par[1] < -5, -5, ifelse( tr_par[1] > 5, 5, tr_par[1]) )
        }

        z0      <- bc(par=tr_par,y=y0,jackup=jackup) + y_added2
        z_ms    <- c(mean(z0),sd(z0))
        z       <- (z0-z_ms[1])/z_ms[2]
        d_abs   <- abs( d_bc(par=tr_par,y0,jackup=jackup)*d_sc(par=z_ms,z0) )#y0 ->z0###### Fix this ############################
        comp    <- (-2)*sum( log( d_abs ) )
      } else {
        z       <- y0
        comp    <- 0
      }

      if( !is.null(weight)){
        m	        <- crossprod( X[,null_dum3], sqrt( weight )*z )#X is sqrt(weight)*X
        zz        <- sum( weight*z^2 )
      } else {
        m	        <- crossprod( X[,null_dum3], z )
        zz        <- sum( z^2 )
      }

      ########################## reml
      m[-(1:nx)]		<- m[ -( 1:nx ) ] * evSqrt[null_dum3[-( 1:nx )]]#[omit_id[-( 1:nx )]]
      b		<- Minv %*% m
      sse		<- zz - 2 * t( b ) %*% m + t( b ) %*% M0 %*% b
      dd		<- abs(sse) + sum( b[ -( 1:nx ) ] ^ 2 )
      if( emet == "reml" ){
        term2	<- ( n - nx ) * ( 1 + log( 2 * pi * dd / ( n - nx ) ) )
      } else if( emet == "ml" ){
        term2	<- n * ( 1 + log( 2 * pi * dd / n ) )
      }
      loglik		<- term1 + term2 + comp + steepen

      if(is.nan(loglik)) loglik <-10^100
      if(is.na(loglik))  loglik <-10^100

      return( loglik[ 1 ] )
    }

    lm_cw<-function(y, M, Minv, m0, k=2,noconst_last=TRUE,y_nonneg=FALSE,jackup,weight=NULL,
                    plim, bc_value=NULL,y_added2=0){

      if(k != 0){
        par00<-rep(c(1,0,0,1),k)
        #lower<-rep(c(1e-10,-3,-Inf,1e-10),k)###### check!!!
        #upper<-rep(c(Inf, 3, Inf, Inf),k)

        if(noconst_last){
          par00<-par00[1:(4*k-2)]
          #lower<-lower[1:(4*k-2)]
          #upper<-upper[1:(4*k-2)]
        }

        if(y_nonneg){
          if(!is.null(bc_value)) {
            par00<-c(par00,bc_value, 0.001)
          } else {
            par00<-c(par00,1, 0.001)
          }
          #lower<-c(lower, -1, 0)
          #upper<-c(upper,  5, 10)
        }
        res    <-optim(par=par00,fn=NLL_sal, y=y,M=M,Minv=Minv,m0=m0,k=k,weight=weight,
                       noconst_last=noconst_last,y_nonneg=y_nonneg,jackup=jackup,
                       plim=plim,bc_value=bc_value,y_added2=y_added2,y_type=y_type)
        #method="L-BFGS-B",lower=lower,upper=upper)
        est0   <-res$par
        est0[est0 >  plim] <-  plim
        est0[est0 < -plim] <- -plim

        est    <-list(NULL)
        for(kk in 1:k){
          if(noconst_last & (kk==k)){
            est0[(4*(kk-1)+1)]<- abs( est0[(4*(kk-1)+1)] )
            est[[kk]]    <- est0[(4*(kk-1)+1):(4*(kk-1)+2)]
          } else {
            est0[(4*(kk-1)+1)]<- abs( est0[(4*(kk-1)+1)] )
            est0[(4*(kk-1)+4)]<- abs( est0[(4*(kk-1)+4)] )
            est[[kk]]    <- est0[(4*(kk-1)+1):(4*(kk-1)+4)]
          }
        }

        if(y_nonneg){
          np_be      <- length(est0)
          bc_par     <- est0[(np_be-1):np_be]
          if(jackup) bc_par[2]  <- abs(bc_par[2])
          if( !is.null( bc_value ) ){
            bc_par[1]<- ifelse( bc_value < -5, -5, ifelse( bc_value > 5, 5, bc_value) )
          } else {
            bc_par[1]<- ifelse( bc_par[1] < -5, -5, ifelse( bc_par[1] > 5, 5, bc_par[1]) )
          }
          est[[kk+1]]<- bc_par
        } else {
          bc_par<-NULL
        }
        res2  <-NLL_sal2(par=unlist(est), y=y, M=M, Minv=Minv, m0=m0,k=k,weight=weight,
                         noconst_last=noconst_last, y_nonneg=y_nonneg,jackup=jackup,
                         plim=plim, bc_value = bc_value,y_added2=y_added2)
        b     <-res2$b
        z     <-res2$z
        loglik<-res2$loglik
        comp  <-res2$comp
        y_ms  <-res2$y_ms
        z_ms  <-res2$z_ms
        #if(y_nonneg){
        #  est <- est[1:(length(est)-2)]
        #}

      } else if(y_nonneg){
        #res   <-optimize(f=NLL_bc, interval=c(-5,5),y=y,M=M,Minv=Minv,m0=m0 )#c(-5,5)
        #bc_par<-res$minimum
        #est   <-NULL

        if(!is.null(bc_value)){
          bc_par0<-c(bc_value, 0.001)
        } else {
          bc_par0<-c(1, 0.001)
        }
        res   <-optim(par=bc_par0,fn=NLL_bc,y=y,M=M,Minv=Minv,m0=m0,jackup=jackup,
                      weight=weight, bc_value=bc_value,y_added2=y_added2)# interval=c(-5,5)
        bc_par<-res$par
        bc_par[2]<-abs(bc_par[2])
        if( !is.null( bc_value ) ){
          bc_par[1]<- ifelse( bc_value < -1, -1, ifelse( bc_value > 5, 5, bc_value) )
        } else {
          bc_par[1]<- ifelse( bc_par[1] < -5, -5, ifelse( bc_par[1] > 5, 5, bc_par[1]) )
        }

        est   <-NULL

        res2  <-NLL_bc2(par=bc_par, y=y, M=M, Minv=Minv, m0=m0,jackup=jackup,weight=weight,
                        bc_value=bc_value,y_added2=y_added2 )
        b     <-res2$b
        z     <-res2$z
        loglik<-res2$loglik
        comp  <-res2$comp
        y_ms  <-NULL
        z_ms  <-res2$z_ms# this part is different from Kaihatsuchushi

      } else {
        est    <- NULL
        res2   <- NULL
        b      <- NULL
        z      <- y
        loglik <- NULL
        comp   <- 0
        bc_par <- NULL
        y_ms  <-NULL
        z_ms  <-NULL
      }

      return(list(vpar=est,bc_par=bc_par,b=b,z=z,loglik=loglik,comp=comp,
                  y_ms=y_ms,z_ms=z_ms))
    }

    Mdet_f0	  	<- function( M, M0, id, par0_sel, emet ){
      if( emet == "ml" ){
        M	  <- M [ id != 0, id != 0 ]
        M0	<- M0[ id != 0, id != 0 ]
        id	<- id[ id != 0 ]
      }
      if(sum( id != par0_sel ) == 0 ){
        term2	  <- NULL
        term3_0	<- M[ id == par0_sel, id == par0_sel ]
      } else {
        M0_sub	<- as.matrix(M0[ id != par0_sel, id != par0_sel ])
        term2		<- determinant( M0_sub )$modulus
        Msub_00	<- M[ id == par0_sel, id == par0_sel ]
        Msub_01	<- M[ id == par0_sel, id != par0_sel ]
        if(sum(id == par0_sel)==1){
          term3_0	<- Msub_00 - c( t( Msub_01 ) %*% solve2( M0_sub, tol = tol ) %*% Msub_01 )
        } else {
          term3_0	<- Msub_00 - Msub_01 %*% solve2( M0_sub, tol = tol ) %*% t( Msub_01 )
        }
      }
      return(list(term2 = term2, term3_0 = term3_0))
    }

    Mdet_f	  	<- function( evSqrt, id, term2, term3_0, par0_sel ){
      term1		  <- sum( log( evSqrt ) ) * 2
      if( dim( as.matrix( term3_0 ) )[1] == 1 ){
        term3_0 <- term3_0 + 1/evSqrt[ id[ id != 0 ] == par0_sel ] ^ 2
        term3_0 <- as.matrix(term3_0)
      } else {
        diag( term3_0 ) <- diag( term3_0 ) + 1/evSqrt[ id[ id != 0 ] == par0_sel ] ^ 2
      }

      term3	  <- determinant( term3_0 )$modulus
      if( is.null( term2 )  ){
        Mdet    <- term1 + term3
      } else {
        Mdet	  <- term1 + term2 + term3
      }
      return(Mdet)
    }

  }

  {
    solve2    <-function(x, tol=tol){#.Machine$double.eps
      test    <- try(xinv<-chol2inv(chol(x)), silent=TRUE)
      if( inherits(test, "try-error") ){
        xinv  <-solve(x, tol=tol)
      }
      return(xinv)
    }

    reg_schur    <- function(x11inv,x11det,x12,x22,xy1, xy2, tol = tol){
      x11inv_x12    <- x11inv %*% x12
      Q_shur        <- x22 - t(x12)%*% x11inv_x12
      Qinv          <- solve(Q_shur, tol = tol)
      x11inv_xy1    <- x11inv%*%xy1
      x12_x11inv_xy1<-t(x12)%*%x11inv_xy1
      Qinv_xy2  <- Qinv %*% xy2
      Qinv_xxxy <- Qinv%*%x12_x11inv_xy1
      b1        <- x11inv_xy1 + x11inv_x12%*%Qinv_xxxy - x11inv_x12 %*% Qinv_xy2
      b2        <- Qinv_xy2 - Qinv_xxxy
      b         <- c(b1,b2)
      Mdet      <- x11det + determinant(Q_shur)$modulus[1]
      return(list(b=b,Q_shur=Q_shur,Mdet=Mdet))
    }

    reg_schur2   <- function(b0, b0_x,x11det,x12,x22,xy2,tol = tol){
      Q_shur    <- x22 - t(x12)%*% b0_x
      Qinv      <- solve(Q_shur, tol = tol)
      Qinv_xy2  <- Qinv %*% xy2
      Qinv_xxxy <- Qinv%*%t(x12)%*%b0
      b1        <- b0 + b0_x%*%Qinv_xxxy - b0_x %*% Qinv_xy2
      b2        <- Qinv_xy2 - Qinv_xxxy
      b         <- c(b1,b2)

      Mdet      <- x11det + determinant(Q_shur)$modulus[1]
      trans     <- list(x12=x12,Qinv=Qinv)
      return(list(b=b,Q_shur=Q_shur, Qinv=Qinv,Mdet=Mdet,trans=trans))
    }

    solve_shur<-function(x11inv,x12,x22, tol = tol){#xy1, xy2,
      Q_shur    <- x22 - t(x12)%*% x11inv %*% x12
      Qinv      <- solve(Q_shur, tol = tol)
      x11inv_x2 <- x11inv %*% x12
      x11inv_x2_Q<- x11inv_x2 %*% Qinv
      AA        <-  x11inv + x11inv_x2_Q %*% t(x11inv_x2)
      BB        <- -x11inv_x2_Q
      CC        <- t(BB)
      DD        <- Qinv
      M_inv   <- rbind(cbind(AA,BB),cbind(CC,DD))
      return(list(M_inv=M_inv,Q_shur=Q_shur))
    }

    reg_schur2_x<- function(x12,x22, b0_x, xy2_x){# Qinv,
      Q_shur    <- x22 - t(x12)%*% b0_x
      Qinv      <- solve(Q_shur, tol = tol)
      Qinv_xx   <- Qinv %*% xy2_x
      Qinv_xxxx <- Qinv%*%t(x12)%*%b0_x
      b1_x      <- b0_x + b0_x%*%Qinv_xxxx - b0_x %*% Qinv_xx
      b2_x      <- Qinv_xx - Qinv_xxxx
      b_x       <- rbind(b1_x,b2_x)
      return(list(b_x=b_x))
    }

    bse_vc_adj    <- function(meig,bse_vc){
      cminfun  <-function(par, sfsf, bse_vc){
        val    <-bse_vc^2 + par*( 1 - sfsf )
        weighted_corr <- cov.wt(data.frame(val,sfsf), wt = sfsf, cor = TRUE)
        weighted_corr$cor[1,2]^2
      }
      opt_v    <- optimize(interval=c(0,1.5),cminfun,sfsf=meig$other$sfsf,bse_vc=bse_vc)
      add      <- sqrt(bse_vc^2 + opt_v$minimum*(1 - meig$other$sfsf )) - bse_vc
      add[add > 0.5*bse_vc]<- 0.5*bse_vc[add > 0.5*bse_vc]
      return(list(add=add,opt_v=opt_v))
    }

    proc_evSqrt   <- function(nsv,ng,nnxf,nnsv,ntv,ntcv,M,par0,id,meig,
                              par0_sel=Inf,null_dum=NULL){
      evSqrt	 <- NULL
      par0_sq	 <- par0 ^ 2
      omit_list<-list(NULL)

      seq_id  <- 1:nsv
      if(!is.null(null_dum)) seq_id  <- seq_id[ null_dum[ seq_id ] == 0 ]
      for( i in seq_id ){
        evv  	  <- meig$ev ^ par0_sq[ nsv + i ] * sum( meig$ev ) / sum( meig$ev ^ par0_sq[ nsv + i ] )
        evSqrt	<- c( evSqrt, par0_sq[ i ] * sqrt( evv ) )
        omit_list[[i]]<- evv/max(evv) >= erat_min
      }
      M0	<- M
      for( j in seq_id ){
        if( j != par0_sel ){
          id_sub	<- ( id == j )
          diag( M0[ id_sub, id_sub ] )<-
            diag( M0[ id_sub, id_sub ] ) + 1/evSqrt[ id[ - ( 1:nx ) ] == j ] ^ 2
        }
      }

      if( ng != 0 ){
        seq_id  <- 1:ng
        if(!is.null(null_dum)) seq_id  <- seq_id[ null_dum[ nsv + seq_id ] == 0 ]

        for( i2 in seq_id ){
          xgg	<- rep( 1, sum( id == ( nsv + i2 ) ) )
          evSqrt	<- c( evSqrt, par0_sq[ 2 * nsv + i2 ] * xgg )
          omit_list[[nsv + i2]]<- rep(TRUE, length(xgg))
        }
        for( j2 in seq_id ){
          if( j2 != ( par0_sel - nsv ) ){
            id_sub	<- ( id == ( j2 + nsv ) )
            if( sum(id_sub)==1 ){
              M0[ id_sub, id_sub ] <- M0[ id_sub, id_sub ] + 1/evSqrt[ id[ - ( 1:nx ) ] == (j2 + nsv) ] ^ 2
            } else {
              diag( M0[ id_sub, id_sub ] )<-
                diag( M0[ id_sub, id_sub ] ) + 1/evSqrt[ id[ - ( 1:nx ) ] == (j2 + nsv) ] ^ 2
            }
          } else {
            id_sub	<- ( id == ( j2 + nsv ) )
            g_add   <- (1/evSqrt[ id[ - ( 1:nx ) ] == (j2 + nsv) ] ^ 2)/100

            if( sum(id_sub)==1 ){
              M0[ id_sub, id_sub ]        <- M0[ id_sub, id_sub ] + g_add
            } else {
              diag( M0[ id_sub, id_sub ] )<- diag( M0[ id_sub, id_sub ] ) + g_add
            }
          }
        }
      }

      if( nnxf != 0 ){
        seq_id  <- 1:nnxf
        if(!is.null(null_dum)) seq_id  <- seq_id[ null_dum[ nsv+ng+seq_id ] == 0 ]

        for( i2 in seq_id ){
          xgg	<- rep( 1, sum( id == ( nsv+ ng + i2 ) ) )
          evSqrt	<- c( evSqrt, par0_sq[ (2 * nsv + ng ) + i2 ] * xgg )
          omit_list[[nsv+ ng + i2]]<- rep(TRUE, length(xgg))
        }

        for( j2 in seq_id ){
          if( j2 != ( par0_sel - nsv - ng ) ){
            id_sub	<- ( id == ( j2 + nsv + ng ) )
            diag( M0[ id_sub, id_sub ] )<-
              diag( M0[ id_sub, id_sub ] ) + 1/evSqrt[ id[ - ( 1:nx ) ] == (j2 + nsv + ng) ] ^ 2
          }
        }
      }

      if( nnsv != 0 ){
        seq_id  <- 1:nnsv
        if(!is.null(null_dum)) seq_id  <- seq_id[ null_dum[ nsv+ng+nnxf+seq_id ] == 0 ]

        for( i2 in seq_id ){
          xgg	<- rep( 1, sum( id == ( nsv + ng + nnxf + i2 ) ) )
          evSqrt	<- c( evSqrt, par0_sq[ (2 * nsv + ng + nnxf ) + i2 ] * xgg )
          omit_list[[nsv + ng + nnxf + i2]]<- rep(TRUE, length(xgg))
        }

        for( j2 in seq_id ){
          if( j2 != ( par0_sel - nsv - ng - nnxf ) ){
            id_sub	<- ( id == ( j2 + nsv + ng + nnxf) )
            diag( M0[ id_sub, id_sub ] )<-
              diag( M0[ id_sub, id_sub ] ) + 1/evSqrt[ id[ - ( 1:nx ) ] == (j2 + nsv + ng + nnxf) ] ^ 2
          }
        }
      }

      if( ntv != 0 ){
        seq_id  <- 1:ntv
        if(!is.null(null_dum)) seq_id  <- seq_id[ null_dum[ nsv+ng+nnxf+nnsv+seq_id ] == 0 ]

        for( i2 in seq_id ){
          n_base  <- 2 * nsv + ng + nnxf + nnsv
          evv  	  <- meig$ev_z[[1]] ^ par0_sq[ n_base + ntv + i2 ] * sum( meig$ev_z[[1]] ) / sum( meig$ev_z[[1]] ^ par0_sq[ n_base+ntv+i2 ] )
          evSqrt	<- c( evSqrt, par0_sq[ n_base + i2 ] * sqrt( evv ) )
          omit_list[[nsv + ng + nnxf + nnsv + i2]]<- evv/max(evv) >= erat_min
        }
        for( j2 in seq_id ){
          if( j2 != ( par0_sel - nsv - ng - nnxf - nnsv ) ){
            id_sub	<- ( id == ( j2 + nsv + ng + nnxf + nnsv) )
            diag( M0[ id_sub, id_sub ] )<-
              diag( M0[ id_sub, id_sub ] ) + 1/evSqrt[ id[ - ( 1:nx ) ] == (j2 + nsv + ng + nnxf + nnsv) ] ^ 2
          }
        }
      }

      if( ntcv != 0 ){
        seq_id  <- 1:ntcv
        if(!is.null(null_dum)) seq_id  <- seq_id[ null_dum[ nsv+ng+nnxf+nnsv+ntv+seq_id ] == 0 ]

        for( i2 in seq_id ){
          n_base  <- 2 * nsv + ng + nnxf + nnsv + 2*ntv
          evv  	  <- meig$ev_z[[2]] ^ par0_sq[ n_base + ntcv + i2 ] * sum( meig$ev_z[[2]] ) / sum( meig$ev_z[[2]] ^ par0_sq[ n_base+ntcv+i2 ] )
          evSqrt	<- c( evSqrt, par0_sq[ n_base + i2 ] * sqrt( evv ) )
          omit_list[[nsv + ng + nnxf + nnsv + ntv + i2]]<- evv/max(evv) >= erat_min
        }
        for( j2 in seq_id ){
          if( j2 != ( par0_sel - nsv - ng - nnxf - nnsv - ntv) ){
            id_sub	<- ( id == ( j2 + nsv + ng + nnxf + nnsv + ntv) )
            diag( M0[ id_sub, id_sub ] )<-
              diag( M0[ id_sub, id_sub ] ) + 1/evSqrt[ id[ - ( 1:nx ) ] == (j2 + nsv + ng + nnxf + nnsv + ntv) ] ^ 2
          }
        }
      }

      return(list(evSqrt=evSqrt, omit_list=omit_list, M0=M0))
    }

    lik_resf_vc		<- function( par0, par0_est, par0_id, par0_sel, meig, omit_list,#ev, ev_z=NULL, ev_zc=NULL,
                              M, M0inv, M0inv_01, M0inv_00,
                              m, yy, b_01, b_02, n, nx, nsv, ntv, ntcv,nnxf,
                              nnsv, ng, emet, term2, term3_0, null_dum2,
                              id, tr_comp=0, max_sigma=NULL, max_alpha=NULL ){ ##### erat>0.001

      if(!is.null(max_sigma)) par0_est[1]<-min(par0_est[2],max_sigma)
      if(!is.null(max_alpha)) par0_est[2]<-min(par0_est[2],max_alpha)

      par		    <- par0 ^ 2
      par_est		<- par0_est ^ 2
      par[ par0_id == par0_sel ]  <- par_est
      evSqrt	  <- NULL
      for( i in ( 1:nsv )[ null_dum2[ 1:nsv ] == 0 ] ){
        ev_sel  <- meig$ev[ omit_list[[i]] ]
        evv  	  <- ev_sel ^ par[ nsv + i ] * sum( ev_sel ) / sum( ev_sel ^ par[ nsv + i ] )

        evSqrt	<- c( evSqrt, par[ i ] * sqrt( evv ) )

        if(par0_sel==i) omit_sub<- evv/max(evv) >= erat_min
      }

      if( ng != 0 ){
        for( j in (1:ng)[ null_dum2[ (( nsv + 1 ):( nsv + ng )) ] == 0 ] ){
          xgg	  <- rep( 1, sum( id == nsv + j ) )
          evSqrt<- c( evSqrt, par[ 2 * nsv + j ] * xgg )

          if(par0_sel==nsv + j) omit_sub<- rep(TRUE,length(xgg))
        }
      }
      if( nnxf != 0 ){
        for( i2 in (1:nnxf)[ null_dum2[ (nsv+ng+1):(nsv+ng+nnxf) ] == 0 ] ){
          xgg	  <- rep( 1, sum( id == ( nsv+ ng + i2 ) ) )
          evSqrt<- c( evSqrt, par[ 2 * nsv + ng + i2 ] * xgg )

          if(par0_sel==nsv + ng + i2) omit_sub<- rep(TRUE,length(xgg))
        }
      }
      if( nnsv != 0 ){
        for( i2 in (1:nnsv)[ null_dum2[ (nsv+ng+nnxf+1):(nsv+ng+nnxf+nnsv) ] == 0 ] ){
          xgg	  <- rep( 1, sum( id == ( nsv + ng + nnxf + i2 ) ) )
          evSqrt<- c( evSqrt, par[ 2 * nsv + ng + nnxf + i2 ] * xgg )

          if(par0_sel==nsv + ng + nnxf + i2) omit_sub<- rep(TRUE,length(xgg))
        }
      }
      if( ntv != 0 ){
        for( i2 in (1:ntv)[ null_dum2[ (nsv+ng+nnxf+nnsv+1):(nsv+ng+nnxf+nnsv+ntv) ] == 0 ] ){
          n_base  <- 2 * nsv + ng + nnxf + nnsv
          ev_sel  <- meig$ev_z[[1]][ omit_list[[nsv+ng+nnxf+nnsv+i2]] ]
          evv  	  <- ev_sel ^ par[ n_base + ntv + i2 ] * sum( ev_sel ) / sum( ev_sel ^ par[ n_base+ntv+i2 ] )
          evSqrt	<- c( evSqrt, par[ n_base + i2 ] * sqrt( evv ) )
          if(par0_sel==nsv+ng+nnxf+nnsv+i2) omit_sub<- evv/max(evv) >= erat_min
        }
      }
      if( ntcv != 0 ){
        for( i2 in (1:ntcv)[ null_dum2[ (nsv+ng+nnxf+nnsv+ntv+1):(nsv+ng+nnxf+nnsv+ntv+ntcv) ] == 0 ] ){
          n_base  <- 2 * nsv + ng + nnxf + nnsv + 2*ntv
          ev_sel  <- meig$ev_z[[2]][ omit_list[[nsv+ng+nnxf+nnsv+ntv+i2]] ]
          evv  	  <- ev_sel ^ par[ n_base + ntcv + i2 ] * sum( ev_sel ) / sum( ev_sel ^ par[ n_base+ntcv+i2 ] )
          evSqrt	<- c( evSqrt, par[ n_base + i2 ] * sqrt( evv ) )
          if(par0_sel==nsv+ng+nnxf+nnsv+ntv+i2) omit_sub<- evv/max(evv) >= erat_min
        }
      }

      sel     <- rep(TRUE,length(evSqrt))
      sel[which(id[id!=0]==par0_sel)]<-omit_sub
      #sel     <- which(unlist(omit_list)==TRUE)#!(ev_rat < erat_min & id[id!=0] == par0_sel)
      sel_sub <- omit_sub#sel[id[id!=0] == par0_sel]
      sel_all <- c(rep(TRUE,nx),sel)
      evv     <- evv[sel]
      evSqrt  <- evSqrt[sel]
      id      <- id[sel_all]
      m       <- m[sel_all]
      M       <- M[sel_all,sel_all]
      b_01    <- b_01[sel_all]
      b_02    <- b_02[sel_sub]

      if(length(sel_sub)==1){
        M0inv_00<-as.matrix(M0inv_00)[sel_sub,sel_sub]
        M0inv_01<-as.matrix(M0inv_01)[sel_all,sel_sub]
        term3_0 <-as.matrix(term3_0)[sel_sub,sel_sub]
      } else {
        M0inv_00<-M0inv_00[sel_sub,sel_sub]
        M0inv_01<-M0inv_01[sel_all,sel_sub]
        term3_0 <- term3_0[sel_sub,sel_sub]
      }
      M0		  <- M

      Mdet		<- Mdet_f( id = id, par0_sel=par0_sel, term2 = term2, term3_0 = term3_0, evSqrt = evSqrt )
      for( j in 1:max( id ) ){
        if( sum( id == j )==1){
          M[ id == j, id == j ]       <- M[ id == j, id == j ] + 1/evSqrt[ id[ -c( 1:nx ) ] == j ] ^ 2
        } else {
          diag( M[ id == j, id == j ])<- diag( M[ id == j, id == j ] ) + 1/evSqrt[ id[ -c( 1:nx ) ] == j ] ^ 2
        }
      }

      if( dim(as.matrix(M0inv_00))[1]==1){
        M0inv_00	      <- M0inv_00 + evSqrt[ id[ id != 0 ] == par0_sel ] ^ 2
      } else {
        diag(M0inv_00)	<- diag( M0inv_00 ) + evSqrt[ id[ id != 0 ] == par0_sel ] ^ 2
      }

      b_02_b<- solve2( M0inv_00, tol = tol ) %*% b_02
      b_02_c<- M0inv_01 %*% b_02_b
      b		  <- b_01 - b_02_c
      sse		<- yy - 2 * t( b ) %*% m + t( b ) %*% M0 %*% b
      dd		<- abs(sse) + sum( ( b[ -( 1:nx ) ] / evSqrt ) ^ 2 )
      if( emet == "reml" ){
        loglik	<- Mdet + ( n - nx ) * ( 1 + log( 2 * pi * dd / ( n - nx ) ) ) - 2*tr_comp
      } else if( emet == "ml" ){
        loglik	<- Mdet + n * ( 1 + log( 2 * pi * dd / n ) ) - 2*tr_comp
      }
      return( loglik )
    }

    lik_resf_vc0	<- function( par0, M, m, yy, n, nx, nsv, ntv, ntcv, meig, omit_list,
                              ng, nnxf, nnsv, emet, null_dum4, tr_comp=0 ){
      par		<- par0 ^ 2
      evSqrt		<- NULL
      for( i in ( 1:nsv )[ null_dum4[ 1:nsv ] == 0 ] ){
        ev_sel  <- meig$ev[ omit_list[[i]] ]
        evv  	  <- ev_sel ^ par[ nsv + i ] * sum( ev_sel ) / sum( ev_sel ^ par[ nsv + i ] )
        evSqrt	<- c( evSqrt, par[ i ] * sqrt( evv ) )
      }

      if( ng != 0 ){
        for( j in (1:ng)[ null_dum4[ (( nsv + 1 ):( nsv + ng )) ] == 0 ] ){#null_dum4[ -(1:nsv) ] == 0
          xgg	<- rep( 1, sum( id == nsv + j ) )
          evSqrt	<- c( evSqrt, par[ 2 * nsv + j ] * xgg )
        }
      }
      if( nnxf != 0 ){
        for( i2 in (1:nnxf)[ null_dum4[ (nsv+ng+1):(nsv+ng+nnxf) ] == 0 ] ){
          xgg	<- rep( 1, sum( id == ( nsv+ ng + i2 ) ) )
          evSqrt	<- c( evSqrt, par[ 2 * nsv + ng + i2 ] * xgg )
        }
      }
      if( nnsv != 0 ){
        for( i2 in (1:nnsv)[ null_dum4[ (nsv+ng+nnxf+1):(nsv+ng+nnxf+nnsv) ] == 0 ] ){
          xgg	<- rep( 1, sum( id == ( nsv + ng + nnxf + i2 ) ) )
          evSqrt	<- c( evSqrt, par[ 2 * nsv + ng + nnxf + i2 ] * xgg )
        }
      }
      if( ntv != 0 ){
        for( i2 in (1:ntv)[ null_dum4[ (nsv+ng+nnxf+nnsv+1):(nsv+ng+nnxf+nnsv+ntv) ] == 0 ] ){
          n_base  <- 2 * nsv + ng + nnxf + nnsv
          ev_sel  <- meig$ev_z[[1]][ omit_list[[nsv+ng+nnxf+nnsv+i2]] ]
          evv  	  <- ev_sel ^ par[ n_base + ntv + i2 ] * sum( ev_sel ) / sum( ev_sel ^ par[ n_base+ntv+i2 ] )
          evSqrt	<- c( evSqrt, par[ n_base + i2 ] * sqrt( evv ) )
        }
      }
      if( ntcv != 0 ){
        for( i2 in (1:ntcv)[ null_dum4[ (nsv+ng+nnxf+nnsv+ntv+1):(nsv+ng+nnxf+nnsv+ntv+ntcv) ] == 0 ] ){
          n_base  <- 2 * nsv + ng + nnxf + nnsv + 2*ntv
          ev_sel  <- meig$ev_z[[2]][ omit_list[[nsv+ng+nnxf+nnsv+ntv+i2]] ]
          evv  	  <- ev_sel ^ par[ n_base + ntcv + i2 ] * sum( ev_sel ) / sum( ev_sel ^ par[ n_base+ntcv+i2 ] )
          evSqrt	<- c( evSqrt, par[ n_base + i2 ] * sqrt( evv ) )
        }
      }

      M[ -( 1:nx ), -( 1:nx ) ]	<- t(M[ -( 1:nx ), -( 1:nx ) ] * evSqrt ) * evSqrt
      M[ -( 1:nx ),    1:nx   ]	<-   M[ -( 1:nx ),    1:nx   ] * evSqrt
      M[    1:nx  , -( 1:nx ) ]	<- t(M[ -( 1:nx ),    1:nx   ] )
      M0		<- M
      #diag( M [ -( 1:nx ), -( 1:nx ) ] ) <- diag( M[ -( 1:nx ), -( 1:nx ) ] ) + 1
      diag( M )[ -( 1:nx ) ]  <- diag( M )[ -( 1:nx ) ] + 1

      m[-(1:nx)]		<- m[ -( 1:nx ) ] * evSqrt
      test			<-try(Minv	<- solve2( M, tol = tol ))
      if( inherits(test, "try-error") ){#class(test)[1]=="try-error"
        loglik  	<- Inf
      } else {
        b		<- Minv %*% m
        sse		<- yy - 2 * t( b ) %*% m + t( b ) %*% M0 %*% b
        dd		<- abs(sse) + sum( b[ -( 1:nx ) ] ^ 2 )
        if( emet == "reml" ){
          term1	<- determinant( M )$modulus
          term2	<- ( n - nx ) * ( 1 + log( 2 * pi * dd / ( n - nx ) ) )
        } else if( emet == "ml" ){
          term1	<- determinant( as.matrix( M[ -( 1:nx ), -( 1:nx ) ] ) )$modulus
          term2	<- n * ( 1 + log( 2 * pi * dd / n ) )
        }
        loglik		<- term1 + term2 -2*tr_comp
      }
      return( loglik[ 1 ] )
    }


    lik_resf_vc_int_Schur<- function( par0_est, ev_z,#par0_sel,
                                      M_t,MM_t,M_old,M_x_t,M_old2_inv,M_old2_det,
                                      evSqrt_sub,m_t,m_old,yy, n, nx, emet, tr_comp=0,max_alpha=NULL ){ ##### erat>0.001

      if(!is.null(max_alpha)) par0_est[2]<-min(par0_est[2],max_alpha)

      par0_sq	    <- par0_est ^ 2
      evv_z0      <- ev_z ^ par0_sq[ 2 ]
      evv_z0_max  <- max(evv_z0[is.finite(evv_z0)])
      evv_z1      <- exp( par0_sq[ 2 ]*log(ev_z)-log(evv_z0_max) )
      evv_z  	    <- evv_z1 * sum( ev_z ) / sum( evv_z1 )
      #evv_z  	    <- ev_z ^ par0_sq[ 2 ] * sum( ev_z ) / sum( ev_z ^ par0_sq[ 2 ] )
      omit_id_int <- evv_z/max(evv_z) >= erat_min

      evSqrt_t    <- par0_sq[ 1 ] * sqrt( evv_z[omit_id_int] )
      evSqrt      <- c(evSqrt_sub, evSqrt_t)
      if(sum(omit_id_int)==1){
        M_t_evSqrt  <- M_t[omit_id_int,omit_id_int]+ 1/evSqrt_t^2
      } else {
        M_t_evSqrt  <- M_t[omit_id_int,omit_id_int]+ diag( 1/evSqrt_t^2 )
      }
      m_t2        <- m_t[omit_id_int]
      n_old       <- nrow(M_old)
      M_x_t2      <- M_x_t[,omit_id_int]
      MM_t2       <- MM_t[c(1:n_old, n_old+which(omit_id_int)),c(1:n_old, n_old+which(omit_id_int))]

      #x11inv=M_old2_inv;x11det=M_old2_det;x12=M_x_t2;x22=M_t_evSqrt;xy1=m_old;xy2=m_t2
      b_hat0      <- reg_schur(x11inv=M_old2_inv,x11det=M_old2_det,x12=M_x_t2,
                              x22=M_t_evSqrt,xy1=m_old, xy2=m_t2, tol = 1e-30)
      b_hat       <- b_hat0$b
      Q_shur      <- b_hat0$Q_shur
      Mdet        <- b_hat0$Mdet + 2*sum(log(evSqrt))
      mm_t        <- c(m_old, m_t2)

      sse		      <- yy - 2 * t( b_hat ) %*% mm_t + t( b_hat ) %*% MM_t2 %*% b_hat
      dd		      <- abs(sse) + sum( ( b_hat[ -( 1:nx ) ]/ evSqrt ) ^ 2 )#
      if( emet == "reml" ){
        loglik	  <- Mdet + ( n - nx ) * ( 1 + log( 2 * pi * dd / ( n - nx ) ) ) - 2*tr_comp
      } else if( emet == "ml" ){
        loglik	  <- Mdet + n * ( 1 + log( 2 * pi * dd / n ) ) - 2*tr_comp
      }
      return(loglik)
    }


    lik_resf_vc_int_Schur2<- function( par0_est, ev_z,#par0_sel,
                                       M_t,MM_t,M_x_t,b0, b0_x, M_old2_det,
                                       evSqrt_sub,m_t,m_old,yy, n, nx, emet, tr_comp=0 ){
      par0_sq	    <- par0_est ^ 2
      evv_z  	    <- ev_z ^ par0_sq[ 2 ] * sum( ev_z ) / sum( ev_z ^ par0_sq[ 2 ] )
      evSqrt_t    <- par0_sq[ 1 ] * sqrt( evv_z )

      if(par0_sel==i) omit_sub<- evv_z/max(evv_z) >= erat_min

      if(1<0){
        sel     <- rep(TRUE,length(evSqrt))
        sel[which(id[id!=0]==par0_sel)]<-omit_sub
        #sel     <- which(unlist(omit_list)==TRUE)#!(ev_rat < erat_min & id[id!=0] == par0_sel)
        sel_sub <- omit_sub#sel[id[id!=0] == par0_sel]
        sel_all <- c(rep(TRUE,nx),sel)
        evv     <- evv[sel]
        evSqrt  <- evSqrt[sel]
        id      <- id[sel_all]
        m       <- m[sel_all]
        M       <- M[sel_all,sel_all]
        M0inv_00<-M0inv_00[sel_sub,sel_sub]
        M0inv_01<-M0inv_01[sel_all,sel_sub]
        b_01    <- b_01[sel_all]
        b_02    <- b_02[sel_sub]
        term3_0 <- term3_0[sel_sub,sel_sub]
        M0		  <- M
      }

      evSqrt      <- c(evSqrt_sub, evSqrt_t)
      M_t_evSqrt  <- M_t + diag(1/evSqrt_t^2)
      b_hat0      <- reg_schur2(b0=b0,b0_x=b0_x,x11det=M_old2_det,
                               x12=M_x_t,x22=M_t_evSqrt,xy2=m_t, tol = 1e-30)
      b_hat       <- b_hat0$b
      Q_shur      <- b_hat0$Q_shur
      Mdet        <- b_hat0$Mdet + 2*sum(log(evSqrt))
      mm_t        <- c(m_old, m_t)

      sse		      <- yy - 2 * t( b_hat ) %*% mm_t + t( b_hat ) %*% MM_t %*% b_hat
      dd		      <- abs(sse) + sum( ( b_hat[ -( 1:nx ) ]/ evSqrt ) ^ 2 )#
      if( emet == "reml" ){
        loglik	  <- Mdet + ( n - nx ) * ( 1 + log( 2 * pi * dd / ( n - nx ) ) ) - 2*tr_comp
      } else if( emet == "ml" ){
        loglik	  <- Mdet + n * ( 1 + log( 2 * pi * dd / n ) ) - 2*tr_comp
      }
      return(loglik)
    }


    sel_basis_vif <-function(test0, vif_max=15){
      tres3<-TRUE
      sel  <-1:dim(test0)[2]
      while(tres3){
        suppressWarnings(cor_test0<-cor(test0))
        testt2<-try(vif <- diag(solve2(cor_test0, tol=tol)), silent=TRUE)
        if( inherits(testt2, "try-error" ) ){#class(testt2)=="try-error"
          suppressWarnings(corr_0<-cor(test0))
          diag(corr_0)<-0
          cor_max <-apply(corr_0,2,max)
          test0   <-test0[, -(which(cor_max == max(cor_max))[1])]
          sel     <-sel[-(which(cor_max == max(cor_max))[1])]

        } else {
          if( max(vif)>vif_max ){
            test0<-test0[,-(which(vif==max(vif))[1])]
            sel     <-sel[-(which(vif==max(vif))[1])]
          } else {
            tres3<-FALSE
          }
        }
      }
      return(list(basis=test0,sel=sel))
    }
  }

  {

  n     	   <- length( y )
  resf_flag  <- ifelse( is.null(x), 1, 0)
  if( !is.null(nongauss)){
    y_type   <- nongauss$y_type
    tr_num   <- nongauss$tr_num
    y_nonneg <- nongauss$y_nonneg
    if( y_nonneg ){
      if( min( y ) < 0 ) stop(" y must be non-negative when y_nonneg = TRUE")
      bc_value<- NULL
    }

  } else {
    y_type   <- "continuous"
    tr_num   <- 0
    y_nonneg <- FALSE
  }

  #if( x_nvc[1]!=FALSE & !is.null(weight) & !is.null(x)){
  #  if( max( na.omit(cor(as.matrix(x),weight)) )>=0.999999 ){
  #    stop("x and weight cannot have the same collumn when x_nvc=TRUE")
  #  }
  #}
  #if( xconst_nvc[1]!=FALSE & !is.null(weight) & !is.null(xconst)){
  #  if( max( na.omit(cor(as.matrix(xconst),weight, na.rm=TRUE)) )>=0.999999 ){
  #    stop("xconst and weight cannot have the same collumn when xconst_nvc=TRUE")
  #  }
  #}

  y_added <- 0
  y_added2<- 0
  if(y_type  == "count"){
    bc_value <- NULL
    y_nonneg <- FALSE
    y_org    <- y
    zrat     <- sum(y_org==0)/length(y_org)
    #y_added  <- 0#0.5
    y_added2 <- -(1+0.5*zrat)/(y_org+0.5)

    if( is.null(offset) ){
      y      <- log( y_org + 0.5 ) + y_added2
    } else {
      if( min( offset ) <= 0 ) stop( "offset must be positive" )
      y      <- log( ( y_org+0.5 )/offset ) + y_added2
    }

    if( is.null(weight) ){
      weight0<- 1
      weight <- y_org + 0.5
    } else {
      weight0<- weight
      weight <- weight0*(y_org+0.5)
    }

  } else if( y_type =="continuous"){
    if( !is.null(offset) ) {
      message("Note: offset is available only for count data (y_type = count). It is ignored")
      offset <- NULL
    }

    if(min(y) <= 0.00000001){
      if(y_nonneg){
        y_added <- 0.00000001
        y       <- y + y_added
        #jackup  <- TRUE
      }
    }

  } else {
    stop( "y_type must be continuous or count" )
  }

  jackup  <- FALSE
  noconst_last <- TRUE

  if( !is.null(weight) ){
    if( min(weight) <= 0) stop( "weight must be positive" )

    if(length(weight)==1){
      w_scale <- 1/weight
    } else {
      w_scale <- n/sum(weight)
    }

    weight  <-c( weight*w_scale )
  } else {
    weight  <- 1
    w_scale <- 1
  }

  nvc_x       <- x_nvc
  nvc_xconst  <- xconst_nvc

  if( is.logical( x_sel ) ){
    allsvc     <- !x_sel
    allsvc_flag<-0
  } else {
    allsvc    <- x_sel
    allsvc_flag<-1
  }

  if( is.logical( x_nvc_sel ) ){
    allnvc_x    <- !x_nvc_sel
    allnvc_x_flag<-0
  } else {
    allnvc_x    <- x_nvc_sel
    allnvc_x_flag<-1
  }

  if( is.logical( xconst_nvc_sel ) ){
    allnvc_xconst     <- !xconst_nvc_sel
    allnvc_xconst_flag<-0
  } else {
    allnvc_xconst<- xconst_nvc_sel
    allnvc_xconst_flag<-1
  }

  if( is.logical( x_sel ) ){
    alltvc     <- !x_sel
    alltvc_flag<- 0
  } else {
    alltvc     <- x_sel
    alltvc_flag<- 1
  }
  if( is.logical( x_sel ) ){
    alltcvc     <- !x_sel
    alltcvc_flag<- 0
  } else {
    alltcvc     <- x_sel
    alltcvc_flag<- 1
  }

  if( x_nvc[1]      == FALSE ) x_nvc_sel = FALSE
  if( xconst_nvc[1] == FALSE ) xconst_nvc_sel = FALSE

  if( method == "reml" ){
    lik_nam	<- "rlogLik"
  } else if( method == "ml" ){
    lik_nam	<- "logLik"
  }

  if(is.null(xconst)) nvc_xconst <- FALSE

  nx0   <- 0
  X1	  <- x
  if( is.null( X1 ) == FALSE ){
    X1	<- as.matrix( X1 )
    if( is.numeric( X1 ) == FALSE ){
      mode( X1 ) <- "numeric"
    }
    x_id	<- apply( X1, 2, sd ) != 0
    nx0	  <- sum( x_id )
    if( nx0 == 0 ){
      X1	  <- NULL
      xname	<- NULL
      x_id	<- NULL
    } else {
      X1	<- as.matrix( X1[ , x_id ] )
      xname	<- names( as.data.frame( X1 ) )
    }

    if( allsvc_flag == 1 ){
      allsvc_id        <- rep(FALSE,length(x_id))
      allsvc_id[allsvc]<- TRUE
      allsvc_id        <- allsvc_id[ x_id ]

      if( !is.null(meig$ev_z)){
        alltvc_id        <- rep(FALSE,length(x_id))
        alltvc_id[allsvc]<- TRUE
        alltvc_id        <- alltvc_id[ x_id ]
      }
      if( !is.null(meig$ev_z) & length(meig$ev_z)>1 ){
        alltcvc_id        <- rep(FALSE,length(x_id))
        alltcvc_id[allsvc]<- TRUE
        alltcvc_id        <- alltcvc_id[ x_id ]
      }

    }

    if( allnvc_x_flag == 1 ){
      allnvc_x_id        <- rep(FALSE,length(x_id))
      allnvc_x_id[allnvc_x]<- TRUE
      allnvc_x_id        <- allnvc_x_id[ x_id ]
    }

  } else {
    xname	<- NULL
    x_id	<- NULL
  }

  Xconst	<- xconst
  if( is.null( xconst ) == FALSE ){
    Xconst	<- as.matrix( Xconst )
    if( is.numeric( Xconst ) == FALSE ){
      mode( Xconst ) <- "numeric"
    }
    xf_id	<- apply( Xconst, 2, sd ) != 0
    nxf	<- sum( xf_id )
    if( nxf == 0 ){
      Xconst	<- NULL
      xfname	<- NULL
      xf_id	  <- NULL
    } else {
      Xconst	<- as.matrix( Xconst[ , xf_id ] )
      xfname	<- names( as.data.frame( Xconst ) )
    }

    if( allnvc_xconst_flag == 1 ){
      allnvc_xconst_id               <- rep(FALSE,length(xf_id))
      allnvc_xconst_id[allnvc_xconst]<- TRUE
      allnvc_xconst_id  <- allnvc_xconst_id[ xf_id ]
    }
  } else {
    xfname	<- NULL
    xf_id	  <- NULL
    nxf	    <- 0
  }

  if(nx0==0) nvc_x <- FALSE
  XX1_0	           <- list(NULL)
  XX1	             <- NULL
  if( is.logical( nvc_x[ 1 ] )&( nvc_x[ 1 ] == TRUE) ) nvc_x     <- 1:nx0

  sel_basis_n <-list( NULL )
  if( is.logical( nvc_x[ 1 ] ) == FALSE ){
    X1_nvc  <- as.matrix(X1[ , nvc_x ])
    xxname	<- names( as.data.frame( X1 ) )[ nvc_x ]
    nnsv    <- length( xxname )
    np_xx     <-apply( X1_nvc,2,function( x ) length( unique( x )))
    np_xx     <-ifelse( np_xx < nvc_num/0.7, round( np_xx * 0.7 ) ,nvc_num )
    np_xx_max <-round( n/nnsv ) - 2
    np_xx[ np_xx > np_xx_max ] <-np_xx_max
    np_xx[ np_xx < 2 ] <- 2

    for( ii in 1:dim( X1_nvc )[ 2 ] ){
      test<-TRUE
      iiii<-0
      while(test){
        kkk  <- np_xx[ ii ]-iiii
        knots<-seq(min( X1_nvc[ , ii ] ),max( X1_nvc[ , ii ] ),len=kkk+2)[2:(kkk+1)]
        testt<- try(XX1_00<- ns( X1_nvc[ ,ii ], knots = knots ), silent=TRUE)#replaced
        test <- inherits(testt, "try-error")#class(testt)[1] == "try-error"
        iiii <- iiii+1
      }

      XX1_00      <- cbind( X1_nvc[,ii], XX1_00)
      sel_basis_n[[ ii ]]    <- 1:dim(XX1_00)[2]
      testt2<-try(XX1_00_mod <- sel_basis_vif( XX1_00, vif_max=15),silent=TRUE)
      if( !inherits(testt2, "try-error") ){#class(testt2)[1] != "try-error"
        XX1_00               <- XX1_00_mod$basis
        sel_basis_n[[ ii ]]  <- XX1_00_mod$sel
        np_xx[ ii ]          <- dim(XX1_00)[2]
      } else {
        np_xx[ ii ]          <- 0
      }
      if( np_xx[ ii ] > 2 ){
        XX1_0[[ii]]<- scale( XX1_00 )
        XX1	      <- cbind( XX1, X1_nvc[, ii ] * XX1_0[[ii]] )
      } else {
        XX1_0[[ii]]<- 0
        np_xx[ ii ]<- 0
      }
    }
  } else {
    nnsv      <- 0
    np_xx     <- NULL
  }

  XXconst_0	  <- list( NULL )
  XXconst	    <- NULL
  if( is.logical( nvc_xconst[ 1 ] )&( nvc_xconst[ 1 ] == TRUE) ) nvc_xconst<- 1:nxf

  sel_basis_c <-list( NULL )
  if( is.logical( nvc_xconst[ 1 ] ) ==FALSE ){
    Xconst_nvc<- as.matrix( Xconst[ , nvc_xconst ] )
    xxfname	  <- names( as.data.frame( Xconst ) )[ nvc_xconst ]
    nnxf      <- length( xxfname )
    np_xxconst       <-apply( Xconst_nvc,2,function( x ) length( unique( x )))
    np_xxconst       <-ifelse( np_xxconst < nvc_num/0.7, round( np_xxconst * 0.7 ) ,nvc_num )
    np_xxconst_max <-round( n/nnxf ) - 2
    np_xxconst[ np_xxconst > np_xxconst_max ] <-np_xxconst_max
    np_xxconst[ np_xxconst < 2 ] <- 2

    for( ii in 1:dim( Xconst_nvc )[ 2 ] ){
      test<-TRUE
      iiii<-0
      while(test){
        kkk      <- np_xxconst[ ii ]-iiii
        knots    <-seq(min( Xconst_nvc[ ,ii ] ),max( Xconst_nvc[ ,ii ] ),len=kkk+2)[2:(kkk+1)]
        #testt<-try(XXconst_00<- scale( poly( Xconst_nvc[ , ii ], np_xxconst[ ii ]-iiii )[ , -1 ] ), silent=TRUE)
        testt<- try(XXconst_00<- ns( Xconst_nvc[ , ii], knots = knots ), silent=TRUE)
        test <- inherits(testt, "try-error")#class(testt)[1] == "try-error"
        iiii <- iiii+1
      }

      XXconst_00<- cbind( Xconst_nvc[,ii], XXconst_00)
      sel_basis_c[[ ii ]]    <- 1:dim(XXconst_00)[2]
      testt2<-try(XXconst_00_mod     <- sel_basis_vif( XXconst_00, vif_max=15),silent=TRUE)
      if( !inherits(testt2, "try-error") ){#class(testt2)[1] != "try-error"
        XXconst_00       <- XXconst_00_mod$basis
        sel_basis_c[[ii]]<- XXconst_00_mod$sel
        np_xxconst[ ii ] <- dim(XXconst_00)[2]
      } else {
        np_xxconst[ ii ]<- 0
      }
      if( np_xxconst[ ii ] > 2 ){
        XXconst_0[[ii]]<- scale( XXconst_00 )
        XXconst	       <- cbind( XXconst, Xconst_nvc[, ii ] * XXconst_0[[ii]] )
        #np_xx[ ii ]<- np_xx[ ii ]  - iiii
      } else {
        XXconst_0[[ii]]<- 0
        np_xxconst[ ii ]<- 0
      }
    }
  } else {
    nnxf         <- 0
    np_xxconst   <- NULL # added
  }

  if( (nx0 > 0) & !is.null(meig$other$coords_z) ){
    dup      <-duplicated(cbind(X1, meig$other$coords_z),MARGIN=2)
    if(sum(dup)>0) stop( "x and coords_z in meig cannot have the same column" )
  }

  if( (nxf > 0) & !is.null(meig$other$coords_z) ){
    dup      <-duplicated(cbind(Xconst, meig$other$coords_z),MARGIN=2)
    if(sum(dup)>0) stop( "xconst and coords_z in meig cannot have the same column" )
  }

  if( ( nx0 > 0 ) & ( nxf > 0 ) ){
    dup      <-duplicated(cbind(X1,Xconst),MARGIN=2)
    if(sum(dup)>0) stop( "x and xconst cannot have the same column" )
  }

  nsv	  <- ifelse( is.null( X1 ),1, dim( X1 )[ 2 ] + 1 )
  nev0	<- min( round( n /nsv ) - 2, length( meig$ev ))
  nev0  <- max(nev0, 3)
  meig$sf	<- meig$sf[ , 1 : nev0 ]
  meig$ev	<- meig$ev[   1 : nev0 ]

  if( !is.null( xgroup ) ){
    xgroup  <- data.frame(xgroup)
    ng	    <- dim( xgroup )[ 2 ]

    if( sum(is.na(xgroup)) > 0) stop("NA is not allowed in xgroup")

    if( ng == 1 ){
      xg_num<- length( unique( xgroup[,1] ) ) - ng
    } else {
      xg_num<- sum( apply(xgroup,2,function(x) length(unique(x)))) - ng
    }

    nev0g   <- round( ( n - xg_num - nxf )/nsv ) - 2
    if( nev0 > nev0g ){
      if( nev0g < 2 ){
        stop("Too many groups. Reduce variables in xgroup")
      } else {
        message( paste( "Note: ", nev0 - nev0g, " eigenvectors are omitted to stablize the estimates",sep=""))
        nev0    <- nev0g
        meig$sf	<- meig$sf[ , 1 : nev0 ]
        meig$ev	<- meig$ev[   1 : nev0 ]
      }
    }
  }

  X2	    <- meig$sf
  #X2_0    <- meig0$sf
  ev	    <- meig$ev
  id	    <- c( rep( 0, nsv ), rep( 0, nxf ), rep( 1, length( ev ) ) )
  if( nsv >= 2 ){
    for( i in 1:( nsv - 1 ) ){
      X2  <- cbind( X2  , X1[, i ] * meig$sf )
      #X2_0<- cbind( X2_0, x0[, i ] * meig0$sf )
      id	<- c( id, rep( i+1, length( ev )))
    }
  }

  if( is.null( xgroup ) == FALSE ){
    xg_id0	<- nsv + 1
    xg_idid2<- NULL
    xg_link_id <- list( NULL )
    xg_levels  <- list(NULL)
    for( ff in 1:ng){
      xg00	<- factor( xgroup[ , ff ] )
      xg_levels[[ff]] <- levels( xg00 )

      xg0	  <- sparse.model.matrix( ~ 0 + xg00 )
      xg0s	<- Matrix::colSums( xg0 )
      xg_idid	<- max( which( xg0s == max( xg0s ) ) )
      Xg0	  <- as.matrix( xg0[ , -xg_idid ] )
      xg_id1	<- rep( xg_id0, dim( Xg0 )[2] )
      if( ff == 1 ){
        Xg	<- Xg0
        xg_id	<- xg_id1
      } else {
        Xg	<- cbind( Xg, Xg0 )
        xg_id	<- c( xg_id, xg_id1 )
      }
      xg_id0	<- xg_id0 + 1
      xg_idid2<-c(xg_idid2, xg_idid)

      xgroup_datt <- data.frame(id=1:length(xgroup[,ff]),xgroup_sub=xgroup[,ff])
      xg00_datt   <- data.frame(id_b_g = 1:length(levels(xg00)), xgroup_sub=levels(xg00))
      xgroup_datt2<- merge(xgroup_datt,xg00_datt,by="xgroup_sub",all.x=TRUE)
      xg_link_id[[ ff ]]  <- xgroup_datt2[order(xgroup_datt2$id),"id_b_g"]
    }
    X2	         <- cbind( X2, Xg )
    id	         <- c(id, xg_id )

  } else {
    ng	<- 0
    xg_levels <- NULL
  }

  pp_id    <-max( id )
  if( is.logical( nvc_xconst[ 1 ] ) == FALSE ){
    xxc_id   <-NULL
    for(pp in 1:length(np_xxconst)){
      xxc_id <- c(xxc_id, rep(pp_id + pp, np_xxconst[pp]))
    }
    id	   <- c( id, xxc_id )
    X2     <- cbind( X2, XXconst )
    pp_id  <- pp_id + length(np_xxconst)
  }

  if( is.logical( nvc_x[ 1 ] ) == FALSE ){
    xx_id   <-NULL
    for(pp in 1:length(np_xx)){
      xx_id <- c(xx_id, rep(pp_id + pp, np_xx[pp]))
    }
    id	   <- c( id, xx_id )
    X2     <- cbind( X2, XX1 )
  }

  pp_id    <-max( id )
  if( !is.null( meig$ev_z ) ){
    id	   <- c( id, rep( pp_id+1, length( meig$ev_z[[1]] )) )
    X2     <- cbind( X2  , meig$sf_z[[1]] )
    #X2_0   <- cbind( X2_0, meig0$sf_t )
    if( nsv >= 2 ){
      for( pp in 2:nsv ){
        X2  <- cbind( X2  , X1[, pp-1 ] * meig$sf_z[[1]] )
        #X2_0<- cbind( X2_0, x0[, pp-1 ] * meig0$sf_t )
        id	<- c( id, rep( pp_id+pp, length( meig$ev_z[[1]] )))
      }
    }
    ntv     <- nsv
  } else {
    ntv     <- 0
  }

  pp_id    <-max( id )
  if( !is.null( meig$ev_z )& length( meig$ev_z )>1  ){
    id	   <- c( id, rep( pp_id+1, length( meig$ev_z[[2]] )) )
    X2     <- cbind( X2  , meig$sf_z[[2]] )
    #X2_0   <- cbind( X2_0, meig0$sf_tc )
    if( nsv >= 2 ){
      for( pp in 2:nsv ){
        X2  <- cbind( X2, X1[, pp-1 ] * meig$sf_z[[2]] )
        #X2_0<- cbind( X2_0, x0[, pp-1 ] * meig0$sf_tc )
        id	<- c( id, rep( pp_id+pp, length( meig$ev_z[[2]] )))
      }
    }
    ntcv    <- nsv
  } else {
    ntcv    <- 0
  }


  if( is.null( Xconst )&is.null( X1 )){
    X0	<- as.matrix( rep(1,n) )
  } else {
    X0	<- as.matrix( cbind( 1, Xconst, X1 ) )
    #X0_0<- as.matrix( cbind( 1, xconst0, x0 ) )
  }
  X	  <- as.matrix( cbind( X0, X2 ) )
  #X_0	<- as.matrix( cbind( X0_0, X2_0 ) )

  X0_org<- X0
  X2_org<- X2
  X_org <- X
  X0    <- sqrt(weight) * X0_org
  X2    <- sqrt(weight) * X2_org
  X     <- sqrt(weight) * X_org

  Mo 	  <- crossprod( X0 )
  nx	  <- dim( X )[ 2 ] - dim( X2 )[ 2 ]
  M     <- crossprod( X )

  y0  <- y
  plim<- 10
  if( ( tr_num > 0 )|( y_nonneg == TRUE ) ){
    Moinv  <- solve2(Mo, tol=tol)
    m0   <- weight * X0_org
    cv_init<-lm_cw(y=y0, M=Mo, Minv=Moinv,m0=m0,k=tr_num,noconst_last=noconst_last,
                   y_nonneg=y_nonneg,jackup=jackup, weight = weight, plim=plim,#, bc_value=bc_value
                   y_added2=y_added2)
    y      <-cv_init$z
    tr_par <-cv_init$vpar[1:tr_num]
    tr_bpar<-cv_init$bc_par
    if( y_nonneg == TRUE ){
      if( jackup ){
        tr_bpar[2]<-abs(tr_bpar[2])
      } else {
        tr_bpar[2]<-0
      }
    }

    tr_npar<-length(unlist(c(tr_par,tr_bpar)))
    #if( jackup==FALSE ) tr_npar<- tr_npar - 1
    tr_comp<-cv_init$comp
    tr_par0<-NULL

  } else {
    tr_par <- NULL
    tr_bpar<- NULL
    tr_comp<- 0
    tr_npar<- 0
  }

  w_y <- weight * y
  mo	<- crossprod( X0_org, w_y )
  m	  <- crossprod( X_org , w_y )
  yy  <- sum( y * w_y )
  resid_naive<- y - X0_org %*% ( solve2( Mo, tol=tol ) %*% mo )
  parVmax_sq <- sqrt( sd( y ) / sd( resid_naive ) )

  if( penalty == "aic" ){
    pen	<- 2
  } else if( penalty == "bic" ){
    pen	<- log( n )
  }

  par0	  <- rep( 1, 2 * nsv )
  par0[1] <- parVmax_sq / 3
  par0_est<- c( 1, 1 )
  par0_id	<- rep( 1:nsv, 2 )
  null_dum<- rep( 0, nsv )

  if( is.null( xgroup ) == FALSE ){
    par0	<- c( par0, rep( parVmax_sq/3, ng ) )
    par0_id	<- c( par0_id , max(par0_id) + ( 1 : ng ) )
    null_dum<- c( null_dum, rep( 0, ng ) )
  }
  if( is.logical( nvc_xconst[ 1 ] ) == FALSE ){
    par0	  <- c( par0, rep( parVmax_sq/10^4.5, nnxf ) )#10^5
    par0_id	<- c( par0_id , max(par0_id) + ( 1:nnxf ) )
    null_dum0<- rep( 0, nnxf )
    null_dum0[np_xxconst==0]<-1
    null_dum<- c( null_dum, null_dum0 )
  }
  if( is.logical( nvc_x[ 1 ] ) == FALSE ){
    par0	<- c( par0, rep( parVmax_sq/10^4.5, nnsv ) )#1
    par0_id	<- c( par0_id , max(par0_id) + ( 1:nnsv ) )
    null_dum0<- rep( 0, nnsv )
    null_dum0[np_xx==0]<-1
    null_dum<- c( null_dum, null_dum0 )
  }
  if( ntv > 0 ){
    par0_sub   <-rep( 1, 2*ntv )
    par0_sub[1]<- parVmax_sq / 3
    par0	  <- c( par0, par0_sub )
    par0_id	<- c( par0_id , max(par0_id) + rep( 1:ntv, 2 ) )
    null_dum<- c( null_dum,  rep( 0, ntv ) )
  }
  if( ntcv >0 ){
    par0_sub   <-rep( 1, 2*ntcv )
    par0_sub[1]<- parVmax_sq / 3
    par0	  <- c( par0, par0_sub )
    par0_id	<- c( par0_id , max(par0_id) + rep( 1:ntcv, 2 ) )
    null_dum<- c( null_dum,  rep( 0, ntcv ) )
  }

  id_exist<-unique(id)[unique(id)>0]
  par0[(par0_id %in% id_exist) ==FALSE]<-0
  Par    <-par0
}

  if(nsv>=7){
    message("Note: At most 5 covariates in x are recommended for stable estimation.")
    message("      Consider using xconst")
  }

  miniter<- ifelse( is.null( miniter ), 3, miniter)
  obj	   <- sd( y )
  LL	   <- NULL
  res_old<- score<-Inf
  iter	 <- 1
  gmess  <- 1
  warn   <- FALSE
  while( (( obj > 0.0001 ) & ( iter <= maxiter ))|( iter <= miniter ) ){#sd( y ) * 0.0001

    if( !is.null(x) ) print( paste( "-------  Iteration ", iter, "  -------", sep = "" ) )

    LL0	<- LL
    for( par0_sel in (1:( nsv+ ng + nnxf + nnsv + ntv + ntcv ))[id_exist]){
      null_dum2<- null_dum
      null_dum2[ par0_sel ] <- 0
      null_dum3_0<- c( rep( 0, nx ), null_dum2[ id ]) == 0

      out_evSqrt<-proc_evSqrt(nsv=nsv,ng=ng,nnxf=nnxf,nnsv=nnsv,ntv=ntv,ntcv=ntcv,
                              M=M,par0=par0,id=id,meig=meig,par0_sel=par0_sel,
                              null_dum=NULL)
      evSqrt    <-out_evSqrt$evSqrt
      M0        <-out_evSqrt$M0
      omit_list <-out_evSqrt$omit_list
      omit_id   <-c(rep(TRUE,nx),unlist(omit_list))#[null_dum3_0]

      null_dum3 <-null_dum3_0 & omit_id
      #null_dum3[omit_id]<-null_dum3_0
      MM	<- M [ null_dum3, null_dum3 ]#[omit_id,omit_id]
      MM0	<- M0[ null_dum3, null_dum3 ]#[omit_id,omit_id]
      mm	<- m [ null_dum3 ]#[omit_id]
      idd	<- id[ null_dum3 ]#[omit_id]

      try1	<- try( M0inv	<- solve2( MM0,  tol = tol ), silent = TRUE )
      try2	<- try( Mdet0	<- Mdet_f0( M = MM, M0 = MM0, id = idd,
                                     par0_sel = par0_sel, emet = method ), silent = TRUE)
      er_dum	<- inherits( try1, "try-error" )|inherits( try2, "try-error" )
      if( er_dum == TRUE ){
        #stop("Singular fit. Simplify the model")
        loglik	<- ( - 1 / 2 ) * res_old
        LL0	<- c( LL0, loglik )
        if( !is.null(x) ) print( paste( par0_sel, "/", nsv + ng + nnxf + nnsv + ntv + ntcv, sep = "" ) )
        next
      }

      term2	  <- Mdet0$term2
      term3_0	<- Mdet0$term3_0

      if( min( par0[ par0_id == par0_sel ] ) >= 1e-5 ){
        par00	<- par0[ par0_id == par0_sel ]
      } else {
        if(par0_sel<=nsv){
          par00_id<- max( which( c(apply(Par[,par0_id == par0_sel], 1, min)) != 0 ) )
        } else if(par0_sel <=  nsv + ng + nnxf + nnsv ){
          par00_id<- max( which( Par[,par0_id == par0_sel] != 0 ) )
        } else {
          par00_id<- max( which( c(apply(Par[,par0_id == par0_sel], 1, min)) != 0 ) )
        }

        par00	<- Par[ par00_id, par0_id == par0_sel ]
      }

      M0inv_01<- M0inv[ 		, idd == par0_sel ]
      M0inv_00<- M0inv[ idd ==par0_sel, idd == par0_sel ]
      b_01	<- M0inv %*% mm
      b_02	<- t( M0inv_01 ) %*% mm

      if( par0_sel == 1){
        llim    <- c( parVmax_sq / 1000, 1e-03)
        ulim    <- c( parVmax_sq, 5 )
        omethod <- "L-BFGS-B"

        ################# changed !!!!!!! 2024/11/14
        res    	<- optim( fn = lik_resf_vc, par00, par0 = par0, M = MM, M0inv = M0inv,
                          meig=meig, omit_list=omit_list,
                          M0inv_01 = M0inv_01, M0inv_00 = M0inv_00, b_01 = b_01, b_02 = b_02,
                          term2 = term2, term3_0 = term3_0, m = mm, yy = yy,nnxf=nnxf,nnsv=nnsv,
                          n = n, nx = nx, nsv = nsv, ntv = ntv, ntcv = ntcv, ng = ng, emet = method,
                          id = idd, par0_sel = par0_sel, par0_id = par0_id, null_dum2 = null_dum2,
                          tr_comp=tr_comp,max_alpha=10 )#, max_sigma=5*parVmax_sq,max_alpha=5
        #lower = llim, upper = ulim, method = omethod,
        res_int <- res

      } else if( par0_sel <= nsv ){
        res    	<- optim( fn = lik_resf_vc, par00, par0 = par0, M = MM, M0inv = M0inv,
                          meig=meig, omit_list=omit_list,
                          M0inv_01 = M0inv_01, M0inv_00 = M0inv_00, b_01 = b_01, b_02 = b_02,
                          term2 = term2, term3_0 = term3_0, m = mm, yy = yy,nnxf=nnxf,nnsv=nnsv,
                          n = n, nx = nx, nsv = nsv, ntv = ntv, ntcv = ntcv, ng = ng, emet = method,
                          id = idd, par0_sel = par0_sel, par0_id = par0_id, null_dum2 = null_dum2,
                          tr_comp=tr_comp, max_alpha=10 )
        res$par[2]<-min(10,res$par[2])

      } else if( par0_sel <= nsv+ng ){
        test<-try(res    	<- optimize( f = lik_resf_vc, lower = parVmax_sq/1000, upper = parVmax_sq, par00,
                                       par0 = par0, M = MM, M0inv = M0inv, meig=meig, omit_list=omit_list,
                                       M0inv_01 = M0inv_01, M0inv_00 = M0inv_00, b_01 = b_01, b_02 = b_02,
                                       term2 = term2, term3_0 = term3_0, m = mm, yy = yy,nnxf=nnxf,nnsv=nnsv,
                                       n = n, nx = nx, nsv = nsv, ntv = ntv, ntcv = ntcv, ng = ng, emet = method,
                                       id = idd, par0_sel = par0_sel, par0_id = par0_id, null_dum2 = null_dum2, tr_comp=tr_comp ),silent = TRUE)
        if( inherits( test, "try-error") ){
          res$value	<- Inf
          res$par	  <- 0
        } else {
          res$value	<- c( res$objective )
          res$par	  <- c( res$minimum )
        }

      } else if( par0_sel <= nsv+ng+nnsv+nnxf ) { ############ +nnxf added

        res    	  <- optimize( f = lik_resf_vc, lower = parVmax_sq/10^5, upper = 10^10, par00, par0 = par0,#10^10
                               M = MM, M0inv = M0inv, meig=meig, omit_list=omit_list,
                               M0inv_01 = M0inv_01, M0inv_00 = M0inv_00, b_01 = b_01, b_02 = b_02,
                               term2 = term2, term3_0 = term3_0, m = mm, yy = yy,nnxf=nnxf,nnsv=nnsv,
                               n = n, nx = nx, nsv = nsv, ntv = ntv, ntcv = ntcv, ng = ng, emet = method, id = idd,
                               par0_sel = par0_sel, par0_id = par0_id, null_dum2 = null_dum2, tr_comp=tr_comp )
        res$value	<- c( res$objective )
        res$par	  <- c( res$minimum )
      } else {

        #par0_est<-par00; par0 = par0; M = MM;m = mm;emet = method;id = idd
        #par0_est[2]  <-100
        res    	  <- optim( fn = lik_resf_vc, par00, par0 = par0, M = MM, M0inv = M0inv,
                            meig=meig, omit_list=omit_list, M0inv_01 = M0inv_01, M0inv_00 = M0inv_00,
                            b_01 = b_01, b_02 = b_02, term2 = term2, term3_0 = term3_0, m = mm, yy = yy,
                            nnxf=nnxf,nnsv=nnsv,n = n, nx = nx, nsv = nsv, ntv = ntv, ntcv = ntcv, ng = ng,
                            emet = method,id = idd, par0_sel = par0_sel, par0_id = par0_id,
                            null_dum2 = null_dum2, tr_comp=tr_comp, max_alpha=10 )
        res$par[2]<-min(10,res$par[2])
      }

      if( res$value > res_old ){#( iter > 3 ) & (
        loglik	<- ( - 1 / 2 ) * res_old
      } else {
        if( par0_sel != 1 ){
          MM_ref 	 <- MM [ idd != par0_sel, idd != par0_sel ]
          mm_ref	 <- mm [ idd != par0_sel ]
          null_dum4<- null_dum2
          null_dum4[ par0_sel ] <- 1
          res_ref  <- lik_resf_vc0( par0, M = MM_ref, m = mm_ref, yy = yy,nnxf=nnxf,nnsv=nnsv,
                                    meig=meig, omit_list=omit_list,
                                    n = n, nx = nx, nsv = nsv, ntv = ntv, ntcv = ntcv,
                                    ng = ng, emet = method, null_dum4 = null_dum4, tr_comp=tr_comp )
        } else {
          res_ref	<- Inf
        }

        np_add       <- length(par00)
        flag         <- 0
        if( par0_sel <= nsv ){
          if( allsvc_flag == 0 ){
            allvc  <- allsvc
          } else {
            flag   <-1
            allvc  <- c(TRUE, allsvc_id)[ par0_sel ]
          }
        } else if( ( nsv < par0_sel )&( par0_sel <= nsv + ng )){
          allvc  <- TRUE
        } else if( ( nsv +ng < par0_sel )&( par0_sel <= nsv + ng + nnxf )){
          if( allnvc_xconst_flag == 0 ){
            allvc  <- allnvc_xconst
          } else {
            flag   <-1
            allvc  <- allnvc_xconst_id[ par0_sel - nsv - ng ]
          }
        } else if( ( nsv +ng +nnxf < par0_sel )&( par0_sel <= nsv + ng + nnxf + nnsv) ){
          if( allnvc_x_flag == 0 ){
            allvc  <- allnvc_x
          } else {
            flag   <-1
            allvc  <- allnvc_x_id[ par0_sel - nsv - ng - nnxf ]
          }
        } else if( ( nsv+ng +nnxf+nnsv < par0_sel )&( par0_sel <= nsv+ng+nnxf+nnsv+ntv) ){
          if( alltvc_flag == 0 ){
            allvc  <- alltvc
          } else {
            flag   <-1
            allvc  <- alltvc_id[ par0_sel - nsv - ng - nnxf - nnsv ]
          }
        } else {
          if( alltcvc_flag == 0 ){
            allvc  <- alltcvc
          } else {
            flag   <-1
            allvc  <- alltvc_id[ par0_sel - nsv - ng - nnxf - nnsv - ntv ]
          }
        }

        score_cand <- res$value + pen * (nx + sum(par0[ par0_id != par0_sel ]!=0)+2 + 1 + tr_npar) - sum( log( weight ))#### 241205
        if(score_cand < score){#### 241205
          if( ((flag==0)&( res_ref < res$value + pen * np_add )&(allvc[1]==FALSE))|((flag==1)&(allvc[1]==FALSE))){### reject
            par0[ par0_id == par0_sel ] 	<- 0
            null_dum[ par0_sel ]		<- 1

            loglik 	<- ( -1 / 2 ) * c( res_ref )
            res_old	<- res_ref
            score   <- res_ref+ pen * (nx+sum(par0!=0)+1+tr_npar) - sum( log( weight ))
          } else {      ####### accept
            par0[ par0_id == par0_sel ]	<- res$par
            null_dum[ par0_sel ]		<- 0
            res_old	<- res$value
            loglik 	<- ( -1 / 2 ) * res$value
            score   <- res$value + pen * (nx + sum(par0[ par0_id != par0_sel ]!=0)+2 + 1 + tr_npar) - sum( log( weight ))
          }
        }#### 241205
      }

      LL0	<- c( LL0, loglik )
      if( !is.null(x) ) print( paste( par0_sel, "/", nsv + ng + nnxf + nnsv + ntv + ntcv, sep = "" ) )
      #print(paste0(round(score,3),"***",round(loglik,3),"***",round(tr_comp,3) ))
    }

    ########################################### TWGP: main
    tr_do <- TRUE
    if((tr_num != 0)|( y_nonneg )){ ####### with transformation
      out_evSqrt<-proc_evSqrt(nsv=nsv,ng=ng,nnxf=nnxf,nnsv=nnsv,ntv=ntv,ntcv=ntcv,
                              M=M,par0=par0,id=id,meig=meig,par0_sel=Inf)#,
                              #null_dum=null_dum)
      M0        <-out_evSqrt$M0
      omit_list <-out_evSqrt$omit_list

      null_dum3_0<- c( rep( 0, nx ), null_dum[ id ]) == 0
      omit_id    <- c(rep(TRUE,nx),unlist(omit_list))#[null_dum3]
      null_dum3  <- null_dum3_0 & omit_id

      evSqrt    <-out_evSqrt$evSqrt[ null_dum3[-(1:nx)] ]#[ omit_id[-(1:nx)] ]
      MM	<- M [ null_dum3, null_dum3 ]#[omit_id,omit_id]
      MM0	<- M0[ null_dum3, null_dum3 ]#[omit_id,omit_id]
      mm	<- m [ null_dum3 ]#[omit_id]
      idd	<- id[ null_dum3 ]#[omit_id]
      id_omit1<- c( diff( id ), 1)
      id_omit2_0<- which( id_omit1 != 0 )                           ### changed
      id_omit2<-id_omit2_0[2:(nsv + 1)]                             ### changed
      id_omit1[ id_omit2_0[(id_omit2_0 %in% id_omit2)==FALSE] ] <- 0### changed

      try1	<- try( M0inv	<- solve2( MM0,  tol = tol ), silent = TRUE )
      if( inherits( try1, "try-error") == TRUE ){
        #stop("Singular fit. Simplify the model")
        loglik	<- ( - 1 / 2 ) * res_old
        LL0	<- c( LL0, loglik )
        if( !is.null(x) ) print( paste( par0_sel, "/", nsv + ng + nnxf + nnsv + ntv + ntcv, sep = "" ) )
        tr_do<-FALSE
      }

      if(tr_do){
        ev	    <- ev[   1:sum( id == 1 ) ]
        meig$sf	<- meig$sf[ , 1:sum( id == 1 ) ]

        MM[ -( 1:nx ), -( 1:nx ) ]	<- t(MM[ -( 1:nx ), -( 1:nx ) ] * evSqrt ) * evSqrt
        MM[ -( 1:nx ),    1:nx   ]	<-   MM[ -( 1:nx ),    1:nx   ] * evSqrt
        MM[    1:nx  , -( 1:nx ) ]	<- t(MM[ -( 1:nx ),    1:nx   ] )
        MM0		<- MM
        #diag( MM[ -( 1:nx ), -( 1:nx ) ] ) <- diag( MM[ -( 1:nx ), -( 1:nx ) ] ) + 1
        diag( MM )[ -( 1:nx ) ] <- diag( MM )[ -( 1:nx ) ] + 1
        MMinv	<- solve2( MM, tol = tol )
        if( method == "reml" ){
          term1	<- determinant( MM )$modulus
        } else if( method == "ml" ){
          term1	<- determinant( as.matrix( MM[ -( 1:nx ), -( 1:nx ) ] ) )$modulus
        }
      }
    }

    if(tr_do){
      if(tr_num != 0){
        if(is.null(tr_par0)){
          tr_par0  <-rep(c(1,0,0,1),tr_num)
          if(noconst_last) tr_par0<-tr_par0[1:(4*tr_num-2)]
          if(y_nonneg)    tr_par0<-c(tr_par0,1,0.01)
        }

        res      <-optim(par=tr_par0,fn=lik_resf_vc_tr, k=tr_num, y_nonneg=y_nonneg,noconst_last=noconst_last,
                         evSqrt=evSqrt, y0=y0, n=n, nx=nx, emet=method, X=X,M0=MM0,Minv=MMinv,term1=term1,
                         null_dum3=null_dum3,jackup=jackup, plim=plim, y_type=y_type,# omit_id=omit_id,
                         y_added2=y_added2,weight=weight, control=list(maxit=100))#,
        #method="L-BFGS-B",lower=lower,upper=upper)

        score_tr <- res$value + pen * (nx + sum(par0!=0) + 1 + tr_npar) - sum( log( weight ))

        if(score_tr< score){
          res_old	 <- res$value
          loglik	 <- (-1/2)*c( res$value )
          score    <-score_tr
          tr_est0  <-res$par
          tr_est0[tr_est0 >  plim] <-  plim
          tr_est0[tr_est0 < -plim] <- -plim

          tr_par   <-list(NULL)
          for(kk in 1:tr_num){
            if(noconst_last & (kk==tr_num)){
              tr_est0[(4*(kk-1)+1)]<- abs( tr_est0[(4*(kk-1)+1)] )
              tr_par[[kk]]         <- tr_est0[(4*(kk-1)+1):(4*(kk-1)+2)]
            } else {
              tr_est0[(4*(kk-1)+1)]<- abs( tr_est0[(4*(kk-1)+1)] )
              tr_est0[(4*(kk-1)+4)]<- abs( tr_est0[(4*(kk-1)+4)] )
              tr_par[[kk]]         <- tr_est0[(4*(kk-1)+1):(4*(kk-1)+4)]
            }
          }

          if(y_nonneg){
            np_b      <-length(tr_est0)
            tr_bpar   <- tr_est0[(np_b - 1):np_b]
            if( jackup  ){
              tr_bpar[2]<- abs(tr_bpar[2])
            } else {
              tr_bpar[2]<- 0
            }
            if( !is.null( bc_value ) ){
              tr_bpar[1]<- ifelse( bc_value < -5, -5, ifelse( bc_value > 5, 5, bc_value) )
            } else {
              tr_bpar[1]<- ifelse( tr_bpar[1] < -5, -5, ifelse( tr_bpar[1] > 5, 5, tr_bpar[1]) )
            }

          } else {
            tr_bpar   <- NULL
          }

          z0       <- sal_k(par=tr_par,y=y0,k=tr_num,noconst_last=noconst_last,bc_par=tr_bpar,jackup=jackup,y_added2=y_added2)
          y        <- z0$y
          z_ms     <- z0$z_ms
          y_ms     <- z0$y_ms

          d_abs   <- abs( d_sal_k(par=tr_par,y=y0,k=tr_num,noconst_last=noconst_last,
                                  bc_par=tr_bpar,y_ms=y_ms,z_ms=z_ms,jackup=jackup) )
          tr_comp <- sum( log( d_abs ) )
          if( !is.null(weight) ){
            w_y <- weight * y
            mo	<- crossprod( X0_org, w_y )
            m	  <- crossprod( X_org, w_y )
            yy  <- sum( y * w_y )
            parVmax_sq <- sqrt( sd( y ) / sd( y - X0_org %*% ( solve2( Mo, tol=tol ) %*% mo ) ) )
          } else {
            mo	<- crossprod( X0, y )
            m	  <- crossprod( X, y )
            yy  <- sum( y ^ 2 )
            parVmax_sq <- sqrt( sd( y ) / sd( y - X0 %*% ( solve2( Mo, tol=tol ) %*% mo ) ) )
          }

          #print(tr_comp)
        }

      } else if(y_nonneg){ ####### without sal layer (only boxcox)

        if(!is.null(bc_value)){
          bc0_par<-c(bc_value,0.01)
        } else {
          bc0_par<-c(1,0.01)
        }
        res     <-optim(par=bc0_par,fn=lik_resf_vc_tr, k=tr_num, y_nonneg=y_nonneg,noconst_last=noconst_last,
                        evSqrt=evSqrt, y0=y0, n=n, nx=nx, emet=method, X=X,M0=MM0,Minv=MMinv,term1=term1,
                        null_dum3=null_dum3,jackup=jackup, weight=weight, plim=plim,y_added2=y_added2,# omit_id=omit_id,
                        y_type=y_type,lower = c(-2,0.0001), upper =  c(2,100), method="L-BFGS-B")#interval=c(-5, 5),
        if( !is.null( bc_value ) ){
          res$par[1]<- ifelse( bc_value < -5, -5, ifelse( bc_value > 5, 5, bc_value) )
        } else {
          res$par[1]<- ifelse( res$par[1] < -5, -5, ifelse( res$par[1] > 5, 5, res$par[1]) )
        }

        score_tr <- res$value + pen * (nx + sum(par0!=0) + 1 + tr_npar) - sum( log( weight ))

        if(score_tr < score){
          res_old	<- res$value
          loglik	<- (-1/2)*c( res$value )
          score   <- score_tr
          tr_bpar <- res$par#minimum
          if(jackup){
            tr_bpar[2]<-abs(tr_bpar[2])
          } else {
            tr_bpar[2]<-0
          }

          y       <- bc(par=tr_bpar, y=y0,jackup=jackup ) + y_added2
          y_ms    <- c(mean(y),sd(y))
          y       <- (y-y_ms[1])/y_ms[2]
          d_abs   <- abs( d_bc(par=tr_bpar,y0,jackup=jackup)*d_sc(par=y_ms,y) )
          tr_comp <- sum( log( d_abs ) )

          w_y <- weight * y
          mo	<- crossprod( X0_org, w_y )
          m	  <- crossprod( X_org, w_y )
          yy  <- sum( y * w_y )
          parVmax_sq <- sqrt( sd( y ) / sd( y - X0_org %*% ( solve2( Mo, tol=tol ) %*% mo ) ) )
        }
      }# else {
       #
       #res_old	  <- res$value
       #loglik 	  <- ( -1 / 2 ) * res_old
       #score     <- res_old + pen * (nx+sum(par0!=0)+1+tr_npar) - sum( log( weight ))
      #}
    }
    #print(paste0(round(score,3),"***",round(loglik,3),"***",round(tr_comp,3) ))

    ###############################################

    if( iter > 1 ){
      #if( sum( n_omit ) == 0 ){
        obj 	<- abs( loglik - loglik0 )
      #} else {
      #  obj	<- Inf
      #}
    } else {
      obj   <- Inf
    }

    Par	    <- rbind( Par, par0 )
    loglik0	<- loglik
    LL	    <- c( LL, loglik )
    #if( !is.null(x) ) print( round( loglik, 3 ) )
    if( !is.null(x) ) print( paste( toupper( penalty ), ": ", round( score, 3 ), sep = "" ) )
    iter	  <- iter + 1
  }

  {
    par2   	 <- par0 ^ 2
    null_dum3_0<- rep(0, nx)
    evSqrt2	 <- omit_list<-list(NULL)
    for( i in 1:nsv ){
      evv  	<- ev ^ par2[ nsv + i ] * sum( ev ) / sum( ev ^ par2[ nsv + i ] )
      omit_list[[i]]  <- evv/max(evv) >= erat_min
      evSqrt2[[ i ]]	<- par2[ i ] * sqrt( evv[ omit_list[[i]] ] )
      null_dum3_0     <- c(null_dum3_0, rep(null_dum[ i ], sum( omit_list[[i]] ) ) )
    }

    if( ng != 0 ){
      for( i2 in 1:ng ){
        xgg	<- rep( 1, sum( id == ( nsv + i2 ) ) )
        omit_list[[nsv + i2]] <- rep(TRUE, length(xgg))
        evSqrt2[[ nsv + i2 ]]	<- par2[ 2 * nsv + i2 ] * xgg
        null_dum3_0           <- c(null_dum3_0, rep(null_dum[ nsv + i2 ], sum( omit_list[[nsv + i2]] ) ) )
      }
    }

    if( nnxf != 0 ){
      for( i2 in 1:nnxf ){
        xgg	<- rep( 1, sum( id == ( nsv + ng + i2 ) ) )
        omit_list[[nsv + ng + i2]] <- rep(TRUE, length(xgg))
        evSqrt2[[ nsv + ng + i2 ]] <- par2[ 2 * nsv + ng + i2 ] * xgg
        null_dum3_0                <- c(null_dum3_0, rep(null_dum[ nsv + ng + i2 ], sum( omit_list[[nsv + ng + i2]] ) ) )
      }
      id_nxf0    <- nsv + ng
      id_nxf     <- (id_nxf0+1):(id_nxf0 +  nsv)
    }

    id_nsv<-NULL
    if( nnsv != 0 ){
      for( i2 in 1:nnsv ){
        xgg	<- rep( 1, sum( id == ( nsv + ng + nnxf + i2 ) ) )
        omit_list[[nsv + ng + nnxf + i2]] <- rep(TRUE, length(xgg))
        evSqrt2[[ nsv + ng + nnxf + i2 ]]	<- par2[ 2 * nsv + ng + nnxf + i2 ] * xgg
        null_dum3_0 <- c(null_dum3_0, rep(null_dum[ nsv + ng + nnxf + i2 ], sum( omit_list[[nsv + ng + nnxf + i2]] ) ) )
      }
      id_nsv0    <- nsv + ng + nnxf
      id_nsv     <- c(NA,(id_nsv0+1):(id_nsv0 + nnsv))
    }

    id_ntv<-NULL
    if( ntv != 0 ){
      for( i2 in 1:ntv ){
        n_base  <- 2 * nsv + ng + nnxf + nnsv
        evv  	  <- meig$ev_z[[1]] ^ par2[ n_base + ntv + i2 ] * sum( meig$ev_z[[1]] ) / sum( meig$ev_z[[1]] ^ par2[ n_base+ntv+i2 ] )
        omit_list[[nsv + ng + nnxf + nnsv + i2 ]] <- evv/max(evv) >= erat_min
        evSqrt2[[ nsv + ng + nnxf + nnsv + i2  ]]	<- par2[  n_base + i2 ] * sqrt( evv[ omit_list[[nsv + ng + nnxf + nnsv + i2 ]] ] )
        null_dum3_0 <- c(null_dum3_0, rep(null_dum[ nsv + ng + nnxf + nnsv + i2 ], sum( omit_list[[nsv + ng + nnxf + nnsv + i2]] ) ) )
      }
      id_ntv0        <- nsv + ng + nnxf + nnsv
      id_ntv         <- (id_ntv0+1):(id_ntv0 + ntv)#nsv + ng + nnxf + nnsv
      id_ntv_interact<- NULL
      if(meig$other$interact) id_ntv_interact<- (id_ntv0 + ntv + ntcv + 1):(id_ntv0 + ntv + ntcv + ntv)

    } else {
      id_ntv         <- NULL
      id_ntv_interact<- NULL
    }

    if( ntcv != 0 ){
      for( i2 in 1:ntcv ){
        n_base  <- 2 * nsv + ng + nnxf + nnsv + 2*ntv
        evv  	  <- meig$ev_z[[2]] ^ par2[ n_base + ntcv + i2 ] * sum( meig$ev_z[[2]] ) / sum( meig$ev_z[[2]] ^ par2[ n_base+ntcv+i2 ] )
        omit_list[[nsv + ng + nnxf + nnsv + ntv + i2 ]] <- evv/max(evv) >= erat_min
        evSqrt2[[  nsv + ng + nnxf + nnsv + ntv + i2 ]]	<- par2[  n_base + i2 ] * sqrt( evv[ omit_list[[nsv + ng + nnxf + nnsv + ntv + i2 ]] ] )
        null_dum3_0<- c(null_dum3_0, rep(null_dum[ nsv + ng + nnxf + nnsv + ntv + i2 ], sum( omit_list[[nsv + ng + nnxf + nnsv + ntv + i2]] ) ) )
      }

      id_ntcv0        <- nsv + ng + nnxf + nnsv + ntv
      id_ntcv         <- (id_ntcv0+1):(id_ntcv0 + ntcv)
      id_ntcv_interact<- NULL
      if(meig$other$interact) id_ntcv_interact<- (id_ntcv0 + ntcv + ntv + 1):(id_ntcv0 + ntcv + ntv + ntcv)

    } else {
      id_ntcv         <- NULL
      id_ntcv_interact<- NULL
    }
    evSqrt    <-unlist( evSqrt2 )
    eevSqrt		<- evSqrt[ null_dum3_0[ -( 1:nx )]==0 ]#[ omit_id[ -( 1:nx )] ]

    #null_dum2<- null_dum
    #null_dum3<- c( rep( 0, nx ), null_dum2[ id[id > 0] ]) == 0
    #null_dum4<- null_dum3[c(rep(TRUE,nx),unlist(omit_list))]
    weight_lik<-weight

    omit_id   <- c(rep(TRUE,nx),unlist(omit_list))#[ null_dum3 ]
    null_dum3 <- omit_id
    null_dum3[omit_id]<- null_dum3_0==0

    MM		    <- M[null_dum3,null_dum3]
    MM[ -( 1:nx ), -( 1:nx ) ]	<- t( MM[ -( 1:nx ), -( 1:nx ) ] * eevSqrt ) * eevSqrt
    MM[ -( 1:nx ),    1:nx   ]	<-    MM[ -( 1:nx ),    1:nx   ] * eevSqrt
    MM[    1:nx  , -( 1:nx ) ]	<- t( MM[ -( 1:nx ),    1:nx   ] )

    diag( MM )[ -( 1:nx )] <- diag( MM)[-(1:nx)] + 1
    MMinv	    <- solve2( MM, tol = tol )

    M0        <- M[null_dum3,null_dum3]
    diag(M0)[ -(1:nx)]<- diag(M[null_dum3,null_dum3])[ -(1:nx)] + 1/eevSqrt^2#1/evSqrt[ null_dum3[ -( 1:nx )] ]^2
    MM0	      <- M0
    idd	      <- id[null_dum3]

    if( method == "reml" ){
      term1	<- determinant( MM )$modulus
    } else if( method == "ml" ){
      term1	<- determinant( as.matrix( MM[ -( 1:nx ), -( 1:nx ) ] ) )$modulus
    }

    mm		          <- m [null_dum3]
    mm[ -( 1:nx ) ]	<- mm[ -( 1:nx ) ] * eevSqrt

    #X_org           <- X_org[,null_dum3][,omit_id]
    #X               <- X[,null_dum3][,omit_id]
    #null_dum3       <- null_dum3[omit_id]
  }

  ############################################# inserted (NEW version)

  par0_nonInt     <- par0
  ii_sel_all      <- omit_id_t_list<-NULL
  sf_z_list       <- ev_z_list<-ev_sel_list<-sf_z_a_list<-list(NULL)
  if( meig$other$interact ){
    n_evsqrt2     <- length(evSqrt2)
    st_threshold0 <- 0.25
    x_sel         <- NULL
    evz_num       <- NULL
    if(ntv > 0){
      evz_num     <- c(evz_num, rep(1,ntv))
      x_sel       <- c(x_sel, (1:ntv) -1 )
      ev_z_all    <- meig$ev_z[[ 1 ]] %x% meig$ev
      ev_z_max    <- max(ev_z_all)
      ev_min_val  <- st_threshold0*ev_z_max
      if( sum( ev_z_all > ev_min_val ) > meig$other$interact_max_dim ){
        ev_min_val <- ev_z_all[order(ev_z_all,decreasing=TRUE)][meig$other$interact_max_dim]
      }

      ev_sel      <- ev_z_all >= ev_min_val
      ev_z_id     <- rep(1:length(meig$ev_z[[ 1 ]]), each=length(meig$ev))[ev_sel]
      ev_id       <- rep(1:length(meig$ev), length(meig$ev_z[[ 1 ]]))[ev_sel]
      ev_sel_list[[1]]<-ev_sel
      ev_z_list[[1]]<- ev_z_all[ev_sel]
      sf_z_list[[1]]<- meig$sf[ , ev_id ]*meig$sf_z[[ 1 ]][ , ev_z_id ]
    }
    if(ntcv > 0){
      evz_num     <- c(evz_num, rep(2,ntcv))
      x_sel       <- c(x_sel, (1:ntcv) -1 )

      ev_z_all     <- meig$ev_z[[ 2 ]] %x% meig$ev
      ev_z_max     <- max( ev_z_all )
      ev_min_val   <- st_threshold0*ev_z_max
      if( sum( ev_z_all > ev_min_val ) > meig$other$interact_max_dim ){
        ev_min_val <- ev_z_all[order(ev_z_all,decreasing=TRUE)][meig$other$interact_max_dim]
      }

      ev_sel       <- ev_z_all > ev_min_val
      ev_z_id      <- rep(1:length(meig$ev_z[[ 2 ]]), each=length(meig$ev))[ev_sel]
      ev_id        <- rep(1:length(meig$ev), length(meig$ev_z[[ 2 ]]))[ev_sel]
      ev_sel_list[[2]]<-ev_sel
      ev_z_list[[2]]<- ev_z_all[ev_sel]
      sf_z_list[[2]]<- meig$sf[ , ev_id ]*meig$sf_z[[ 2 ]][ , ev_z_id ]
    }

    M_t  <- M0_t   <- M_x_t <-m_t <- eev_z<- list(NULL)
    if( !is.null(evz_num) ){
      ii_list      <- 1:length(evz_num)
      for( ii in ii_list ){

        x_sf_z     <- sf_z_list[[ evz_num[ii] ]]
        if( x_sel[[ii]] >  0 ) x_sf_z  <- X1[, x_sel[ii]]*x_sf_z
        M_t[[ii]]  <- M0_t[[ii]] <-crossprod( sqrt(weight)*x_sf_z )  #crossprod_ag(E_a=sf_z_a,wei=st_wei,id_a=id_a)#
        M_x_t[[ii]]<- crossprod( X[,null_dum3], sqrt(weight)*x_sf_z )#   crossprod_ag(Z=NULL,Z_dt=Z_dt_list[[ x_sel[ii] + 1 ]],E_a=sf_z_a,wei=st_wei,id_a=id_a)#
        m_t[[ii]]  <- crossprod( x_sf_z, w_y )    #crossprod( sf_z_a[id_a,], st_wei*y )#
        eev_z[[ii]]<- ev_z_list[[ evz_num[ii] ]]
      }
      sf_z      <- NULL
    }

    omit_id_t_list <- lapply(eev_z,function(x) rep(TRUE,length(x)))
    omit_id_all    <- rep(TRUE,sum(null_dum3))#rep(TRUE,sum(null_dum5==0))
    score_pre      <- score
    loglik_pre     <- loglik
    #par0_sel_t  <- max(id) + 1
    ii_run         <- TRUE
    ii_sel_all     <- NULL
    ii_use         <- rep(TRUE,length(ii_list))
    ii_iter        <- 1
    while( ii_run ){
      #print("------ iter -------")
      if( !is.null(x) ) print( paste( "-------  Iteration (interaction term) ", ii_iter, "  -------", sep = "" ) )
      if( ii_iter == 1 ){
        M_old      <- M[ null_dum3, null_dum3 ]
        m_old		   <- m[ null_dum3 ]
        idd_old    <- idd
        evSqrt_sub <- eevSqrt#evSqrt[ null_dum3[ -( 1:nx )] ]
        M_old2		 <- M_old
        for( j in 1:max( idd_old ) ){
          if( sum( idd_old == j )==1){
            M_old2[ idd_old == j, idd_old == j ]       <- M_old2[ idd_old == j, idd_old == j ] + 1/evSqrt_sub[ idd_old[ -c( 1:nx ) ] == j ] ^ 2
          } else {
            diag( M_old2[ idd_old == j, idd_old == j ])<- diag( M_old2[ idd_old == j, idd_old == j ] ) + 1/evSqrt_sub[ idd_old[ -c( 1:nx ) ] == j ] ^ 2
          }
        }

      } else {
        M_old        <-M_new
        m_old        <-m_new
        idd_old       <-idd_new
        evSqrt_sub   <-evSqrt_sub_new
      }

      ############################ trial
      if(ii_iter==1){
        M_old2_inv  <- solve(M_old2,tol=tol)
        M_old2_det  <- determinant(M_old2)$modulus[1]
      } else {
        ii_pre           <- ii_sel_all[length(ii_sel_all)]
        M_t_evSqrt       <- M_t[[ii_pre]]

        diag(M_t_evSqrt) <- diag(M_t[[ii_pre]]) + 1/evSqrt_sub_t^2
        shur0            <- solve_shur(x11inv=M_old2_inv,x12=M_x_t[[ii_pre]][1:n_old,],x22=M_t_evSqrt, tol = tol)
        M_old2_inv       <- shur0$M_inv
        M_old2_det       <- M_old2_det + determinant(shur0$Q_shur)$modulus[1]#1/determinant(M_old2_inv)$modulus[1]
      }

      n_old              <- ncol(M_old)
      SS_tv<-score<-M_new<- M0_new <- ii_sel <- score_pre0<-NULL
      for(ii in ii_list[ ii_use ]){
        n_t          <- ncol(M_t[[ii]])
        MM_t         <- matrix(NA, n_old + n_t,n_old + n_t)
        MM_t[1:n_old, 1:n_old]<-M_old
        MM_t[1:n_old, (n_old+1):(n_old+n_t)]<- M_x_t[[ii]]
        MM_t[(n_old+1):(n_old+n_t), 1:n_old]<-t(M_x_t[[ii]])
        MM_t[(n_old+1):(n_old+n_t), (n_old+1):(n_old+n_t)]<- M_t[[ii]]

        par0_alpha_init<-log(erat_min)/( log( min(eev_z[[ii]])/max(eev_z[[ii]]) ))
        par0_est <- c(1,par0_alpha_init)
        max_alpha<- max(10,par0_alpha_init+1)
        #ev_z=eev_z[[ii]]
        #M_t=M_t[[ii]];M_x_t=M_x_t[[ii]];m_t=m_t[[ii]];emet = method
        res     <-optim(fn=lik_resf_vc_int_Schur, par0_est, ev_z=eev_z[[ii]],
                        M_t=M_t[[ii]],MM_t=MM_t,M_old=M_old,M_x_t=M_x_t[[ii]],
                        M_old2_inv=M_old2_inv,M_old2_det=M_old2_det,
                        evSqrt_sub=evSqrt_sub, m_t=m_t[[ii]],m_old=m_old,
                        yy = yy, n = n, nx = nx, emet = method,tr_comp=tr_comp,
                        max_alpha=max_alpha,control = list(maxit = 50) )
        res$par[2]<- min(res$par[2], max_alpha)
        par0_t    <- c(par0, res$par)
        loglik0   <- ( -1/2 ) * res$value
        score0    <- res$value + pen * (nx+sum(par0_t!=0)+1+tr_npar) - sum( log( weight ))
        if(score0 <  score_pre){
          ii_sel  <- ii
          M_new   <- MM_t
          score   <- score0
          score_pre<-score0
          loglik  <- loglik0
          par0_sub<- res$par ^ 2
        }
        #print(score0)
        SS_tv<-c(SS_tv,score0)
      }

      if( !is.null(score) ){
        evv_sub_t0    <- eev_z[[ii_sel]] ^ par0_sub[ 2 ]
        evv_sub_t0_max<- max(evv_sub_t0[is.finite(evv_sub_t0)])
        evv_sub_t1    <- exp( par0_sub[ 2 ]*log( eev_z[[ii_sel]] )-log( evv_sub_t0_max ) )
        evv_sub_t  	  <- evv_sub_t1 * sum( eev_z[[ii_sel]] ) / sum( evv_sub_t1 )# to avoid overflow
        #evv_sub_t  	  <- eev_z[[ii_sel]] ^ par0_sub[ 2 ] * sum( eev_z[[ii_sel]] ) / sum( eev_z[[ii_sel]] ^ par0_sub[ 2 ] )
        omit_id_t     <- omit_id_t_list[[ii_sel]] <- evv_sub_t/max(evv_sub_t) >= erat_min
        omit_id_all   <- c( omit_id_all, omit_id_t )

        evSqrt_sub_t  <- c(par0_sub[ 1 ] * sqrt( evv_sub_t[omit_id_t] ))
        evSqrt_sub_new<- c(evSqrt_sub, evSqrt_sub_t)
        M_new         <- M_new[c(1:n_old, n_old+which(omit_id_t)),c(1:n_old, n_old+which(omit_id_t))]
        m_new         <- c(m_old,m_t[[ii_sel]][omit_id_t])
        idd_new        <- c(idd_old,rep(max(idd_old+1), sum(omit_id_t)))#length(m_t[[ii_sel]][omit_id_t] )
        ii_sel_all    <- c( ii_sel_all, ii_sel )
        evSqrt2[[n_evsqrt2+ii_sel]]<- evSqrt_sub_t
        par0          <- c(par0, res$par^2) ####################

        ###### update M_t
        M_t[[ii_sel]] <- as.matrix( M_t[[ii_sel]][omit_id_t,omit_id_t] )

        ###### update M_x_t
        sf_z_sel      <- sf_z_list[[ evz_num[ii_sel] ]][,omit_id_t]
        if( x_sel[[ii_sel]] >  0 ) sf_z_sel  <- X1[, x_sel[ii_sel]]*sf_z_sel

        ii_use[ii_sel] <- FALSE
        for( ii in ii_list ){
          sf_z         <- sf_z_list[[ evz_num[ii] ]][, omit_id_t_list[[ ii ]] ]
          if( x_sel[[ii]] >  0 ) sf_z  <- X1[, x_sel[ii]]*sf_z
          if( ii == ii_sel){
            M_x_t[[ii]] <- rbind( as.matrix( M_x_t[[ii]][, omit_id_t_list[[ ii ]] ] ), crossprod( sf_z_sel, weight * sf_z ) )#sqrt(weight)
          } else {
            M_x_t[[ii]] <- rbind( as.matrix( M_x_t[[ii]] ), crossprod( sf_z_sel, weight * sf_z ) )#sqrt(weight)
          }

        }
        omit_id_all  <- rep( TRUE, sum(omit_id_all) )
        ii_iter      <- ii_iter + 1
        if( !is.null(x) ) print( paste( toupper( penalty ), ": ", round( score, 3 ), sep = "" ) )

        if( sum(ii_use)==0 ){
          M_out      <- M_new
          m_out      <- m_new
          idd_out    <- idd_new
          evSqrt_out <- evSqrt_sub_new
          break
        }

      } else {
        M_out        <- M_old
        m_out        <- m_old
        idd_out      <- idd_old
        evSqrt_out   <- evSqrt_sub
        ii_run       <- FALSE
      }
    }
    #if(parallel) stopCluster(ncores)
  }

  if( meig$other$interact ){
    id_nonInt_max<-Inf
    if( ii_iter > 0 ){
      for( ii in ii_sel_all ){
        if( x_sel[[ii]] == 0 ) X_org <- cbind(X_org, sf_z_list[[ evz_num[ii] ]] )
        if( x_sel[[ii]] >  0 ) X_org <- cbind(X_org, X1[, x_sel[ii]]*sf_z_list[[ evz_num[ii] ]] )
        null_dum3<- c( null_dum3, omit_id_t_list[[ii]] )
      } #sum(null_dum3)+length(nnn_dum3)=385(288+97)
      X        <- sqrt(weight) * X_org

      M        <- M_out#[omit_id,omit_id]     # 443x443
      m        <- m_out#[omit_id]     # 443
      id_nonInt_max<-max(idd)
      idd      <- idd_out    # 458
      eevSqrt  <- evSqrt_out#[omit_id[-(1:nx)]]
      MM		   <- M
      MM[ -( 1:nx ), -( 1:nx ) ]	<- t( MM[ -( 1:nx ), -( 1:nx ) ] * eevSqrt ) * eevSqrt
      MM[ -( 1:nx ),    1:nx   ]	<-    MM[ -( 1:nx ),    1:nx   ] * eevSqrt
      MM[    1:nx  , -( 1:nx ) ]	<- t( MM[ -( 1:nx ),    1:nx   ] )
      diag( MM )[ -( 1:nx ) ] <- diag( MM )[ -( 1:nx ) ] + 1

      M0      <- M
      diag(M0)[ -(1:nx)]<- diag(M)[ -(1:nx)] + 1/eevSqrt^2
      MMinv	  <- solve2( MM, tol = tol )
      MM0	    <- M0

      if( method == "reml" ){
        term1	<- determinant( MM )$modulus
      } else if( method == "ml" ){
        term1	<- determinant( as.matrix( MM[ -( 1:nx ), -( 1:nx ) ] ) )$modulus
      }

      mm		          <- m#[ null_dum3 ]
      mm[ -( 1:nx ) ]	<- mm[ -( 1:nx ) ] * eevSqrt

      null_dum     <-c(null_dum, ifelse( ii_list %in% ii_sel_all,0,1))
    }

    idd          <- idd[idd <= id_nonInt_max]
    id_interact0 <- min( c(id_ntv_interact, id_ntcv_interact), na.rm=TRUE) - 1
    for( iii in ii_sel_all){
      idd        <- c(idd, rep(id_interact0+iii,  sum(omit_id_t_list[[iii]]) ) )
    }
  }
  omit_list       <- c(omit_list, omit_id_t_list)

  par0_int        <- par0_int_dat<-NULL
  if( length( par0_nonInt ) < length(par0) ){
    par0_int<- par0[-(1:length(par0_nonInt))]
    par0_int_dat<-data.frame(x_sel=x_sel[ii_sel_all]+1,int_sel=evz_num[ii_sel_all],
                             random_SE=par0_int[seq(1,length(par0_int),2)],
                             alpha=par0_int[seq(2,length(par0_int),2)])
  }

  b               <- MMinv %*% mm
  b[ -( 1:nx ) ]	<- b[ -( 1:nx ) ] * eevSqrt
  b_sub           <- b[ -( 1:nx ) ][eevSqrt>0]/eevSqrt[eevSqrt>0]

  X3              <- X
  #X3[,null_dum3][,omit_id][,-( 1:nx )]<- t(t(X[,null_dum3][,omit_id][,-( 1:nx )])* eevSqrt)
  X3[,null_dum3][,-( 1:nx )]<- t(t(X[,null_dum3][,-( 1:nx )])* eevSqrt)
  nxx             <- nx

  ############################################# Poisson estimator
  #ulim      <- NULL
  if( y_type == "count" ){
    pred0_lm   <- X_org[,null_dum3] %*% b ########### [,omit_id] added #########[ , null_dum3 ][,omit_id]

    if( tr_num > 0 ){
      z0_p     <- sal_k(par=tr_par,y=y0,k=tr_num,noconst_last=noconst_last,
                        bc_par=tr_bpar,jackup=jackup,y_added2=y_added2)#
      z_ms     <- z0_p$z_ms
      y_ms     <- z0_p$y_ms
      #ulim     <- max(z0_p$y)
      pred0_ex <- i_sal_k(par=tr_par,y=pred0_lm,k=tr_num,noconst_last=noconst_last,
                          bc_par=tr_bpar,y_ms=y_ms,z_ms=z_ms,jackup=jackup)#y_added2=y_added2,ulim=ulim
    } else {
      pred0_ex  <-pred0_lm
      #ulim      <- max(bc(par=tr_bpar,y=y0,jackup=jackup))
      #pred0_ex  <- i_bc(par=tr_bpar,y=pred0_lm,jackup=jackup)#,ulim=ulim
    }

    pred0_ex   <- exp(pred0_ex)
    if( !is.null( offset ) ) pred0_ex <- offset*pred0_ex

    #########
    #weight00    <- c(  weight0 * pred0_ex )
    #w_scale_test<- n/sum(weight00)
    #weight_test <- c( weight00*w_scale_test )

    #y_p        <- sqrt( weight_test ) * (pred0_lm + (y_org - pred0_ex)/pred0_ex)
    #X_test     <- sqrt( weight_test ) * X_org
    #M_test     <- crossprod( X_test[,null_dum3] )
    #m_test	   <- crossprod( X_test[,null_dum3] , y_p )

    #MM_test		  <- M_test
    #MM_test[ -( 1:nx ), -( 1:nx ) ]	<- t( MM_test[ -( 1:nx ), -( 1:nx ) ] * eevSqrt ) * eevSqrt
    #MM_test[ -( 1:nx ),    1:nx   ]	<-    MM_test[ -( 1:nx ),    1:nx   ] * eevSqrt
    #MM_test[    1:nx  , -( 1:nx ) ]	<- t( MM_test[ -( 1:nx ),    1:nx   ] )
    #diag( MM_test [ -( 1:nx ), -( 1:nx ) ] ) <- diag( MM_test[ -( 1:nx ), -( 1:nx ) ] ) + 1
    #MMinv_test	  <- solve2( MM_test, tol = 1e-30 )
    #mm_test		            <- m_test
    #mm_test[ -( 1:nx ) ]	<- mm_test[ -( 1:nx ) ] * eevSqrt
    #b_test                <- MMinv_test %*% mm_test
    #b_test[ -( 1:nx ) ]	  <- b_test[ -( 1:nx ) ] * eevSqrt

    #X3_test        <- X_test[,null_dum3]
    #X3_test[,-( 1:nx )]<- t(t(X_test[,null_dum3][,-( 1:nx )])* eevSqrt)#[, null_dum3 ][, null_dum3 ]

    #pred0_lm_test  <- X_org[,null_dum3] %*% b_test##############[,omit_id] added#[ , null_dum3 ][,omit_id]

    #if( tr_num > 0 ){
    #  pred0_ex_test <- i_sal_k(par=tr_par,y=pred0_lm_test,k=tr_num,noconst_last=noconst_last,
    #                           bc_par=tr_bpar,y_ms=y_ms,z_ms=z_ms,jackup=jackup)#y_added2=y_added2,ulim=ulim
    #} else {
    #  pred0_ex_test <- pred0_lm_test
    #}

    #pred0_ex_test   <- exp(pred0_ex_test)
    #if( !is.null( offset ) ) pred0_ex_test <- offset*pred0_ex_test

    ######## devance
    nz_dum     <- (y_org > 0) & ( pred0_ex[,1] > 0)# ( pred0_ex_test[,1] > 0) &
    dev1        <- 2*sum(y_org[nz_dum]*log(y_org[nz_dum]/pred0_ex[nz_dum])) - 2*sum(y_org - pred0_ex)
    #dev1_test   <- 2*sum(y_org[nz_dum]*log(y_org[nz_dum]/pred0_ex_test[nz_dum])) - 2*sum(y_org - pred0_ex_test)
    #if(dev1 >= dev1_test){
    #  X    <- X_test
    #  X3   <- X3_test
    #  M    <- m_test
    #  m    <- m_test
    #  MM   <- MM_test
    #  mm   <- mm_test
    #  MMinv<- MMinv_test
    #  weight<-weight_test
    #  w_scale<-w_scale_test
    #  b    <- b_test
    #  pred0_ex<-pred0_ex_test
    #  dev1 <- dev1_test
    #}

    ######## null deviance
    y_mean    <- mean( y_org[nz_dum] )
    dev0      <- 2*sum(y_org[nz_dum]*log(y_org[nz_dum]/y_mean )) - 2*sum(y_org - y_mean )
    dev_rat0  <- (sum( dev0 ) - sum( dev1 ))/sum( dev0 )
    dev_rat   <- ifelse( dev_rat0 < 0, 0, dev_rat0)

    ######## mean squared error
    rmse      <- sqrt( mean((y_org - pred0_ex)^2) )
  }

  if( y_type == "continuous"){
      pred0	  <- X_org[,null_dum3 ] %*% b

      wsq_y   <- sqrt(weight)*y
      resid_w <- wsq_y - X[,null_dum3 ] %*% b
      SSE     <- sum( resid_w ^ 2 )
      SSY		  <- sum( ( wsq_y - mean( wsq_y ) ) ^ 2 )
      sig		  <- SSE /( n - nxx )

      resid	  <- y - pred0
      SSE     <- sum( resid ^ 2 )
      SSY		  <- sum( ( y - mean( y ) ) ^ 2 )
      sig_org <- SSE /( n - nxx )

      b_cov		<- sig * MMinv
      b_cov2  <- b_cov[ b!=0, b!=0 ]
      pred0_se00<-colSums( t( X3[,null_dum3 ][ ,b!=0 ] ) * ( b_cov2 %*% t( X3[,null_dum3 ][ ,b!=0 ] ) ) )
      pred0_se00[pred0_se00<0] <- 0
      pred0_se0<-sqrt(pred0_se00)
      if( meig$other$fast ){
        pred0_se0<- pred0_se0 + bse_vc_adj(meig=meig,bse_vc=pred0_se0 )$add
        pred0_se <- sqrt(pred0_se0^2 + sig)
      } else {
        pred0_se <- sqrt( pred0_se0^2 + sig )
      }
      pred0_se<- pred0_se / sqrt( weight )

  } else if( y_type == "count" ){
    pred0	  <- X_org[,null_dum3 ] %*% b######### [,omit_id] added #####################################[ , null_dum3 ][,omit_id]
    resid	  <- y - pred0
    sig_org <- sum( ( y_org - pred0_ex )^2 / pred0_ex )/( n - nxx )

    wsq_y   <- sqrt(weight)*y
    resid_w <- wsq_y - X[,null_dum3 ]  %*% b#
    SSE     <- sum( resid_w ^ 2 )
    SSY		  <- sum( ( wsq_y - mean( wsq_y ) ) ^ 2 )
    sig     <- SSE /( n - nxx )

    b_cov		<- sig * MMinv
    b_cov2  <- b_cov[ b!=0, b!=0 ]

    ######## (Gaussian approx)
    pred0_se0<- sqrt( colSums( t( X3[,null_dum3 ][ ,b!=0 ] ) *( b_cov2 %*% t( X3[,null_dum3 ][ ,b!=0 ] ) ) ) )
    if( meig$other$fast ){
      pred0_se0<- pred0_se0 + bse_vc_adj(meig=meig,bse_vc=pred0_se0 )$add
      pred0_se <- sqrt(pred0_se0^2 + sig)
    } else {
      pred0_se <- sqrt( pred0_se0^2 + sig )
    }

  }

  bse		  <- sqrt( diag( b_cov ) )
  nev		  <- length( ev )

  b_vc_s0	  <- b_vc_t0	<- b_vc_n0  <- NULL
  B_vc_s	  <- B_vc_t	  <- B_vc_n  	<- NULL
  moran_vc  <- rep( 0, nsv )
  b_vc		  <- matrix( 0, nrow = n, ncol = nsv )
  bse_vc		<- matrix( 0, nrow = n, ncol = nsv )
  bb		    <- bb_cov		<- list(NULL)

  evSqrts		   <- evSqrts_n <- evSqrts_t <- evSqrts_tc <- evSqrts_t_int <- evSqrts_tc_int<- list(NULL)
  evSqrts[1:max(c(nsv,id_ntv,id_ntcv,id_ntv_interact,id_ntcv_interact),na.rm=TRUE)]<- NA
  evSqrts_n[1:nsv] <- evSqrts_t[1:nsv] <- evSqrts_tc[1:nsv] <- evSqrts_t_int[1:nsv] <- evSqrts_tc_int[1:nsv]<- NA

  if( nnsv == 0 ){
    j0		<- 1
    for( j in 1:nsv ){
      bid0		<- ifelse( j == 1, 1, nxf + j )
      dum_s   <- null_dum[ j ] == 0
      dum_t   <- dum_t_interact   <- FALSE
      dum_tc  <- dum_tc_interact  <- FALSE
      if( !is.null( id_ntv ) ) dum_t <- null_dum[ id_ntv[j]  ] == 0
      if( !is.null( id_ntcv) ) dum_tc<- null_dum[ id_ntcv[j] ] == 0
      if( !is.null( id_ntv_interact  ) ) dum_t_interact <- null_dum[ id_ntv_interact[j]  ] == 0
      if( !is.null( id_ntcv_interact ) ) dum_tc_interact<- null_dum[ id_ntcv_interact[j] ] == 0

      j_t              <- j_t_interact <-NULL
      j_tc             <- j_tc_interact<-NULL
      if( dum_t )  j_t <- id_ntv[j]
      if( dum_tc ) j_tc<- id_ntcv[j]
      if( dum_t_interact ) j_t_interact <-id_ntv_interact[j]
      if( dum_tc_interact) j_tc_interact<-id_ntcv_interact[j]

      if( dum_s|dum_t|dum_tc|dum_t_interact|dum_tc_interact ){
        bid            <-bid0
        bb0            <-b[ bid0 ]
        b_vc[ , j ]    <-b[ bid0 ]
        x_sf           <-matrix(1,nrow=n,ncol=1)
        if( dum_s ){
          bid0_vc_s	   <- which( idd == j )
          bid <- bid_s <- c( bid, bid0_vc_s )
          bb0          <- c(bb0, b[ bid0_vc_s ] )
          meig_sel     <- as.matrix( meig$sf[, omit_list[[j]]] )
          b_vc[ , j ]	 <- b_vc[ , j ] + meig_sel %*% b[ bid0_vc_s ]
          x_sf		     <- cbind(x_sf, t( t( meig_sel ) * evSqrt2[[ j ]] ))
          evSqrts[[ j ]]<- evSqrt2[[ j ]]
          moran_vc[ j ] <- sum(b[ bid0_vc_s ]^2*meig$ev[ omit_list[[j]] ])/(meig$ev[1]*sum(b[ bid0_vc_s ]^2))
        } else {
          evSqrts[[ j ]]<- NA
        }
        if( dum_t ){
          bid0_vc_t	  <- which( idd == j_t )
          bid         <- c( bid, bid0_vc_t )
          bb0         <- c(bb0, b[ bid0_vc_t ] )
          meig_sel    <- as.matrix( meig$sf_z[[1]][, omit_list[[j_t]] ] )
          b_vc[ , j ]	<- b_vc[ , j ] + meig_sel %*% b[ bid0_vc_t ]
          x_sf		    <- cbind(x_sf, t( t( meig_sel ) * evSqrt2[[ j_t ]] ))
          evSqrts[[ j_t ]]<- evSqrts_t[[ j ]]<-evSqrt2[[ j_t ]]
        } else {
          evSqrts_t[[ j ]]<-NA
          if( !is.null(j_t) ){
            evSqrts[[ j_t ]]<- NA
          }
        }
        if( dum_tc ){
          bid0_vc_tc	<- which( idd == j_tc )
          bid         <- c( bid, bid0_vc_tc )
          bb0         <- c(bb0, b[ bid0_vc_tc ] )
          meig_sel    <- as.matrix(meig$sf_z[[2]][, omit_list[[j_tc]] ])
          b_vc[ , j ]	<- b_vc[ , j ] + meig_sel %*% b[ bid0_vc_tc ]
          x_sf		    <- cbind(x_sf, t( t( meig_sel ) * evSqrt2[[ j_tc ]] ))
          evSqrts[[ j_tc ]]<- evSqrts_tc[[ j ]]<-evSqrt2[[ j_tc ]]
        } else {
          evSqrts_tc[[ j ]]<-NA
          if( !is.null(j_tc) ){
            evSqrts[[ j_tc ]]<- NA
          }
        }
        if( dum_t_interact ){
          bid0_vc_t_interact<- which( idd == j_t_interact )
          bid         <- c( bid, bid0_vc_t_interact )
          bb0         <- c(bb0, b[ bid0_vc_t_interact ] )
          meig_sel    <- as.matrix( sf_z_list[[1]][, omit_list[[j_t_interact]] ] )
          b_vc[ , j ]	<- b_vc[ , j ] + meig_sel %*% b[ bid0_vc_t_interact ]
          x_sf		    <- cbind(x_sf, t( t( meig_sel ) * evSqrt2[[ j_t_interact ]] ))
          evSqrts[[ j_t_interact ]]<- evSqrts_t_int[[ j ]]<-evSqrt2[[ j_t_interact ]]
        } else {
          evSqrts_t_int[[ j ]]<-NA
          if( !is.null(j_t_interact) ){
            evSqrts[[ j_t_interact ]]<- NA
          }

        }
        if( dum_tc_interact ){
          bid0_vc_tc_interact<- which( idd == j_tc_interact )
          bid         <- c( bid, bid0_vc_tc_interact )
          bb0         <- c(bb0, b[ bid0_vc_tc_interact ] )
          meig_sel    <- as.matrix( sf_z_list[[2]][, omit_list[[j_tc_interact]] ] )
          b_vc[ , j ]	<- b_vc[ , j ] + meig_sel %*% b[ bid0_vc_tc_interact ]
          x_sf		    <- cbind(x_sf, t( t( meig_sel ) * evSqrt2[[ j_tc_interact ]] ))
          evSqrts[[ j_tc_interact ]]<- evSqrts_tc_int[[ j ]]<-evSqrt2[[ j_tc_interact ]]
        } else {
          evSqrts_tc_int[[ j ]]<-NA
          if( !is.null(j_tc_interact) ){
            evSqrts[[ j_tc_interact ]]<- NA
          }

        }

        x_sf		      <- as.matrix( x_sf )
        bse_vc0       <- colSums( t( x_sf ) * ( b_cov[ bid, bid ] %*% t( x_sf ) ) )
        bse_vc0[ bse_vc0<=0 ] <- min( bse_vc0[bse_vc0>0] )
        bse_vc[ , j ]	<- sqrt( bse_vc0 )
        bb[[ j ]]	    <- bb0
        bb_cov[[ j ]]	<- b_cov[ bid, bid ]

        #if( dum_s & meig$other$fast ){
        #  #bse_vc_s	   <- sqrt( colSums( t( x_sf[,bid %in% bid_s] ) * ( b_cov[ bid_s, bid_s ] %*% t( x_sf[,bid %in% bid_s] ) ) ) )
        #  bse_vc[ , j ]<- bse_vc[,j]+bse_vc_adj( meig=meig,bse_vc=bse_vc[,j] )$add
        #}
        j0		        <- j0 + 1
      } else {
        b_vc[ , j ]	  <- b[ bid0 ]
        moran_vc[ j ] <- NA
        bse_vc[ , j ]	<- sqrt( b_cov[ bid0, bid0 ] )
        bb[[ j ]]	    <- b[ bid0 ]
        bb_cov[[ j ]]	<- b_cov[ bid0, bid0 ]
        evSqrts[[ j ]]<- NA #######################################
      }
    }
  } else {
    b_vc_s0		<- matrix( 0, nrow = n, ncol = nsv )
    b_vc_n0	  <- matrix( 0, nrow = n, ncol = nsv )
    bse_vc_s0	<- matrix( 0, nrow = n, ncol = nsv )
    bse_vc_n0	<- matrix( 0, nrow = n, ncol = nsv )
    #b_vc		  <- matrix( 0, nrow = n, ncol = nsv )
    #bse_vc		<- matrix( 0, nrow = n, ncol = nsv )

    j_nsv     <- 1
    for( j in 1:nsv ){
      bid0		<- ifelse( j == 1, 1, nxf + j )
      dum_s   <- null_dum[ j ] == 0
      dum_nsv <- FALSE
      dum_t   <- dum_t_interact   <- FALSE
      dum_tc  <- dum_tc_interact  <- FALSE
      if( !is.null( id_nsv )&(j>1) ) dum_nsv<- null_dum[ id_nsv[j]  ] == 0
      if( !is.null( id_ntv ) ) dum_t  <- null_dum[ id_ntv[j]  ] == 0
      if( !is.null( id_ntcv) ) dum_tc <- null_dum[ id_ntcv[j] ] == 0
      if( !is.null( id_ntv_interact  ) ) dum_t_interact <- null_dum[ id_ntv_interact[j]  ] == 0
      if( !is.null( id_ntcv_interact ) ) dum_tc_interact<- null_dum[ id_ntcv_interact[j] ] == 0

      j_nsv <- j_t <- j_t_interact <-j_tc<- j_tc_interact<-NULL
      if( dum_nsv ) j_nsv<-id_nsv[j]
      if( dum_t )  j_t <- id_ntv[j]
      if( dum_tc ) j_tc<- id_ntcv[j]
      if( dum_t_interact ) j_t_interact <-id_ntv_interact[j]
      if( dum_tc_interact) j_tc_interact<-id_ntcv_interact[j]

      bid             <- bid0
      bb0             <- b[ bid0 ]
      b_vc[ , j ]     <- b_vc_s0[ , j ]  <- b_vc_n0[ , j ]  <- b[ bid0 ]
      x_sf            <- matrix(1,nrow=n,ncol=1)
      bse_vc[ , j ]	  <- bse_vc_n0[ , j ]<- bse_vc_s0[ , j ]<- sqrt( b_cov[ bid0, bid0 ] )
      evSqrts[[ j ]]  <- evSqrts_n[[ j ]]<- evSqrts_t[[ j ]]<- evSqrts_tc[[ j ]]<- evSqrts_t_int[[ j ]]<- evSqrts_tc_int[[ j ]]<- NA
      if( dum_s|dum_nsv|dum_t|dum_tc|dum_t_interact|dum_tc_interact ){
        if( dum_s ){
          bid0_vc		    <- which( idd == j )
          bid <- bid_s  <- c( bid, bid0_vc )
          bb0           <- c(bb0, b[ bid0_vc ] )
          meig_sel      <- as.matrix( meig$sf[, omit_list[[j]]] )
          b_vc_s0_0     <- meig_sel %*% b[ bid0_vc ]
          b_vc_s0[ , j ]<- b_vc_s0[ , j ] + b_vc_s0_0
          b_vc[ , j ]	  <- b_vc[ , j ] + b_vc_s0_0

          moran_vc[ j ] <- sum(b[ bid0_vc ]^2*meig$ev[ omit_list[[j]] ])/(meig$ev[1]*sum(b[ bid0_vc ]^2))
          x_sf_ss       <- as.matrix( cbind( 1, t( t( meig_sel ) * evSqrt2[[ j ]] ) ))
          x_sf          <- cbind(x_sf, x_sf_ss[,-1] )

          bse_vc_s0[,j]<- sqrt( colSums( t( x_sf_ss ) * (  b_cov[ bid_s, bid_s ] %*% t( x_sf_ss ) ) ) )
          #if( meig$other$fast ){
          #  bse_vc_s0[ , j ]<-bse_vc_s0[,j]+bse_vc_adj( meig=meig,bse_vc=bse_vc_s0[,j] )$add
          #}
          evSqrts[[ j ]]<- evSqrt2[[ j ]]
        } else {
          evSqrts[[ j ]]<- NA
        }
        if( dum_nsv ){#(( 1:nsv ) %in% ( nvc_x + 1 ))[ j ]
          bid0_nvc      <- which(idd == j_nsv )#nsv + ng + nnxf + j - 1
          bid           <- c( bid, bid0_nvc)
          bb0           <- c(bb0, b[ bid0_nvc ] )
          b_vc_n0_0     <- XX1_0[[ j - 1 ]] %*% b[ bid0_nvc ]
          b_vc_n0[ , j ]<- b_vc_n0[ , j ] + b_vc_n0_0
          b_vc[ , j ]	  <- b_vc[ , j ] + b_vc_n0_0

          x_sf_nn       <- as.matrix( cbind( 1, t( t( XX1_0[[ j - 1 ]] ) * evSqrt2[[ j_nsv ]] ) ))
          x_sf          <- cbind(x_sf, x_sf_nn[,-1] )
          b_cov_sub_n0  <- b_cov[ c(bid0, bid0_nvc), c(bid0, bid0_nvc) ]
          bse_vc_n0[ , j ]<- sqrt( colSums( t( x_sf_nn ) * ( b_cov_sub_n0 %*% t( x_sf_nn ) ) ) )
          evSqrts[[j_nsv]]<-evSqrts_n[[ j ]]<- evSqrt2[[ j_nsv ]]
        } else {
          evSqrts_n[[ j ]]<- NA
        }
        if( dum_t ){
          bid0_vc_t	  <- which( idd == j_t )
          bid         <- c( bid, bid0_vc_t )
          bb0         <- c(bb0, b[ bid0_vc_t ] )
          meig_sel    <- as.matrix( meig$sf_z[[1]][, omit_list[[j_t]] ] )
          b_vc[ , j ]	<- b_vc[ , j ] + meig_sel %*% b[ bid0_vc_t ]
          x_sf		    <- cbind(x_sf, t( t( meig_sel ) * evSqrt2[[ j_t ]] ))
          evSqrts[[ j_t ]]<- evSqrts_t[[ j ]]<-evSqrt2[[ j_t ]]
        } else {
          evSqrts_t[[ j ]]<- NA
        }

        if( dum_tc ){
          bid0_vc_tc	<- which( idd == j_tc )
          bid         <- c( bid, bid0_vc_tc )
          bb0         <- c(bb0, b[ bid0_vc_tc ] )
          meig_sel    <- as.matrix(meig$sf_z[[2]][, omit_list[[j_tc]] ])
          b_vc[ , j ]	<- b_vc[ , j ] + meig_sel %*% b[ bid0_vc_tc ]
          x_sf		    <- cbind(x_sf, t( t( meig_sel ) * evSqrt2[[ j_tc ]] ))
          evSqrts[[ j_tc ]]<- evSqrts_tc[[ j ]]<-evSqrt2[[ j_tc ]]
        } else {
          evSqrts_tc[[ j ]]<- NA
        }
        if( dum_t_interact ){
          bid0_vc_t_interact<- which( idd == j_t_interact )
          bid         <- c( bid, bid0_vc_t_interact )
          bb0         <- c(bb0, b[ bid0_vc_t_interact ] )
          meig_sel    <- as.matrix( sf_z_list[[1]][, omit_list[[j_t_interact]] ] )
          b_vc[ , j ]	<- b_vc[ , j ] + meig_sel %*% b[ bid0_vc_t_interact ]
          x_sf		    <- cbind(x_sf, t( t( meig_sel ) * evSqrt2[[ j_t_interact ]] ))
          evSqrts[[ j_t_interact ]]<- evSqrts_t_int[[ j ]]<-evSqrt2[[ j_t_interact ]]
        } else {
          evSqrts_t_int[[ j ]]<- NA
        }
        if( dum_tc_interact ){
          bid0_vc_tc_interact<- which( idd == j_tc_interact )
          bid         <- c( bid, bid0_vc_tc_interact )
          bb0         <- c(bb0, b[ bid0_vc_tc_interact ] )
          meig_sel    <- as.matrix( sf_z_list[[2]][, omit_list[[j_tc_interact]] ] )
          b_vc[ , j ]	<- b_vc[ , j ] + meig_sel %*% b[ bid0_vc_tc_interact ]
          x_sf		    <- cbind(x_sf, t( t( meig_sel ) * evSqrt2[[ j_tc_interact ]] ))
          evSqrts[[ j_tc_interact ]]<- evSqrts_tc_int[[ j ]]<- evSqrt2[[ j_tc_interact ]]
        } else {
          evSqrts_tc_int[[ j ]]<- NA
        }

        x_sf		      <- as.matrix( x_sf )#cbind( 1, sf2 )
        bse_vc0       <- colSums( t( x_sf ) * ( b_cov[ bid, bid ] %*% t( x_sf ) ) )
        bse_vc0[ bse_vc0<=0 ] <- min( bse_vc0[bse_vc0>0] )
        bse_vc[ , j ]	<- sqrt( bse_vc0 )
        bb[[ j ]]	    <- bb0
        bb_cov[[ j ]]	<- b_cov[ bid, bid ]

        #if( meig$other$fast & dum_s ){
        #  bse_vc[ , j ]<-bse_vc[,j] + bse_vc_adj( meig = meig, bse_vc = bse_vc[,j] )$add
        #}
      } else {
        moran_vc[ j ] <- NA
        bb[[ j ]]	    <- bb0
        bb_cov[[ j ]]	<- b_cov[ bid0, bid0 ]
        evSqrts[[ j ]]<- NA ##################################
      }
    }
  }

  if( nnxf != 0 ){
    bf_vc		  <- matrix( 0, nrow = n, ncol = nxf )
    bfse_vc		<- matrix( 0, nrow = n, ncol = nxf )
    bf_s		  <- list(NULL)
    bf_covs		<- list(NULL)
    evSqrts_c	<- list(NULL)
    evSqrts_c[1:nxf]<-NA
    for( j in which( (1:nxf) %in% nvc_xconst ) ){
      #bid0		        <- 1 + j
      if( (null_dum[ nsv + ng + j ] == 0)&( sum(idd == nsv + ng + j)>0 ) ){
        bid0_nvc      <- which( idd == nsv + ng + j )
        bf_vc[ , j ]  <- b[ 1 + j ] + XXconst_0[[ j ]] %*% b[ bid0_nvc ]
        x_sf_n        <- t( t( XXconst_0[[ j ]] ) * evSqrt2[[ nsv + ng + j ]] )
        bid           <-c(1 + j, bid0_nvc)
        b_cov_sub	    <- b_cov[ bid, bid ]
        x_sf		      <- as.matrix( cbind( 1,x_sf_n ))
        bfse_vc[ , j ]<- sqrt( colSums( t( x_sf ) * ( b_cov_sub %*% t( x_sf ) ) ) )

      } else {
        bid            <- 1 + j#bid0
        bf_vc  [ , j ] <- b[ bid ]
        b_cov_sub      <- sqrt( b_cov[ bid, bid ])
        bfse_vc[ , j ] <- b_cov_sub
        bid0_nvc       <- NULL
        x_sf_n         <- NULL
      }

      bf_s[[ j ]]	    <- b[ bid ]
      bf_covs[[ j ]]	<- b_cov_sub
      evSqrts_c[[ j ]]<- evSqrt2[[ nsv + ng  + j ]]
      if( null_dum[ nsv + ng + j ] == 1 ) evSqrts_c[[ j ]]<- NA# 0 #####################################
    }
  } else {
    bf_vc		  <- NULL
    bfse_vc		<- NULL
    bf_s		  <- NULL
    bf_covs		<- NULL
    evSqrts_c	<- NA
  }

  parV		<- par2[ 1:nsv  ]
  parR		<- par2[ (nsv + 1):( 2 * nsv ) ]
  parV_t  <-parR_t<-parV_tc<-parR_tc<-NULL
  if(ntv > 0){
    parV_t	<- par2[ 2*nsv + ng +nnxf + nnsv + 1:ntv  ]
    parR_t	<- par2[ 2*nsv + ng +nnxf + nnsv + (ntv + 1):( 2 * ntv ) ]
  }
  if(ntcv > 0){
    parV_tc	<- par2[ 2*nsv + ng +nnxf + nnsv + 2*ntv +1:ntcv  ]
    parR_tc	<- par2[ 2*nsv + ng +nnxf + nnsv + 2*ntv +(ntcv + 1):( 2 * ntcv ) ]
  }

  Bias  <-0
  if( ng != 0 ){####### Check id2!!! (2023/4)
    id2   <- idd[( idd %in% which(parV!=0) )|( idd > nsv )|( idd == 0 ) ]
    b_g		<- b  [ id2 > nsv ]
    bse_g <- bse[ id2 > nsv ]
    bt_g	<- b_g / bse_g
    bpar_g		<- data.frame( Estimate = b_g, SE = bse_g, t_value = bt_g )

    glist <-unique(id2[id2 > nsv])
    bpar_g2<-list(NULL)
    nulldat<-data.frame(Estimate = 0, SE = NA, t_value = NA)
    for(ggid in 1:ng){
      bpar_g2_0		<- bpar_g[id2[id2 > nsv] == glist[ggid],]
      xg_idid_0   <- xg_idid2[ggid]
      if(xg_idid_0 == 1){
        bpar_g2_0		<- rbind(nulldat,bpar_g2_0)
      } else if( xg_idid_0 == dim(bpar_g2_0)[1] ){
        bpar_g2_0		<- rbind(bpar_g2_0,nulldat)
      } else {
        bpar_g2_0		<- rbind(bpar_g2_0[1:(xg_idid_0-1),],nulldat,bpar_g2_0[-(1:(xg_idid_0-1)),])
      }
      bias  <- mean(bpar_g2_0[xg_link_id[[ggid]],1])
      bpar_g2_0[,1]<- bpar_g2_0[,1] - bias
      Bias  <- Bias + bias

      bpar_g2_0[is.na(bpar_g2_0[,2])==FALSE,3]<- bpar_g2_0[is.na(bpar_g2_0[,2])==FALSE,1]/bpar_g2_0[is.na(bpar_g2_0[,2])==FALSE,2]
      xg00	<- factor( xgroup[ , ggid ] )
      rownames(bpar_g2_0) <- paste( names( as.data.frame(xgroup) )[ ggid ], "_", levels( xg00 ), sep = "" )
      bpar_g2[[ggid]] <- bpar_g2_0
    }

  } else {
    bpar_g2  <- NULL
  }

  if( ng != 0 ){
    parG	<- par2[ (2 * nsv + 1 ):( 2 * nsv + ng ) ]
  } else {
    parG  <- NULL
  }
  if(nnxf != 0 ){
    parNxf<- par2[ (2 * nsv +ng + 1 ):( 2 * nsv + ng + nnxf) ]
  } else {
    parNxf<- NULL
  }
  if(nnsv != 0 ){
    parNsv<- par2[ (2 * nsv + ng + nnxf + 1 ):( 2 * nsv + ng + nnxf + nnsv ) ]
  } else {
    parNsv<- NULL
  }

  nsv2		<- sum( parV != 0 )
  nnxf2   <- sum( parNxf != 0 )
  nnsv2   <- sum( parNsv != 0 )
  ntv2		<- sum( parV_t != 0 )
  ntcv2		<- sum( parV_tc != 0 )
  Xm		  <- X[,null_dum3]
  Xm[ , -( 1:nx ) ] <- t( t( Xm[ , -( 1:nx ) ] ) * eevSqrt )
  np		  <- nxx + 2*nsv2 + ng + nnxf2 + nnsv2 + 2*ntv2 + 2*ntcv2 + 1 + tr_npar

  loglik<-loglik + 0.5*sum( log( weight_lik ))########## check!!
  Xm_org    <- X_org[,null_dum3 ]
  Xm_org[ , -( 1:nx ) ] <- t( t( Xm_org[ , -( 1:nx ) ] ) * eevSqrt )

  AIC		  <- -2 * loglik + np * 2
  BIC		  <- -2 * loglik + np * log( n )
  r2_0		<- 1 - SSE / SSY
  r2		  <- 1 - ( 1- r2_0 ) * ( n - 1 ) / ( n - np - 1 )

  if( is.null( x ) ){
    b_vc[,1]<- b_vc[,1] + Bias
    b_c		  <- b[ 1:( nxf + 1 ) ]
    b_c[1]  <- b_c[1] + Bias
    r       <- b[( nx + 1 ):( nx + nev) ]
    bse_c	  <- bse[ 1:( nxf + 1 ) ]
    bt_c	  <- b_c / bse_c
    df	    <- sum( t( Xm_org ) * ( MMinv %*% t( weight * Xm_org  ) ) )

    bp_c	  <- 2 - 2 * pt( abs( bt_c ), df = n - df )
    b_par		<- data.frame( Estimate = b_c, SE = bse_c, t_value = bt_c, p_value = bp_c )
    rownames( b_par )<- c("(Intercept)", xfname)

  } else {
    if( nxf != 0 ) {
      b_c		<- b  [ 2:( nxf + 1 ) ]
      bse_c	<- bse[ 2:( nxf + 1 ) ]
      bt_c	<- b_c / bse_c
      df	  <- sum( t( Xm_org ) * ( MMinv %*% t( weight * Xm_org  ) ) )
      bp_c	<- 2 - 2 * pt( abs( bt_c ), df = n - df )
      b_par	<- data.frame( Estimate = b_c, SE = bse_c, t_value = bt_c, p_value = bp_c )
      rownames( b_par )<- xfname
    } else {
      df	  <- sum( t( Xm_org ) * ( MMinv %*% t( weight * Xm_org  ) ) ) ######################### time consuming
      b_par	<- NULL
    }

    b_vc[,1]<- b_vc[,1] + Bias
    r       <- NULL
  }

  bt_vc		<- as.matrix( b_vc / bse_vc )
  bp_vc		<- 2 - 2 * pt( abs( bt_vc ), df = n - df )
  b_vc		<- data.frame( b_vc )
  bse_vc		<- data.frame( bse_vc )
  bt_vc		<- data.frame( bt_vc )
  bp_vc		<- data.frame( bp_vc )
  bp_vc_adj<- data.frame( apply(bp_vc, 2, function(x) p.adjust(x, method="BY" ) ) )
  names( b_vc )	  <- c( "(Intercept)", xname )
  names( bse_vc )	<- c( "(Intercept)", xname )
  names( bt_vc )	<- c( "(Intercept)", xname )
  names( bp_vc )	<- c( "(Intercept)", xname )
  names( bp_vc_adj )<- c( "(Intercept)", xname )

  if(is.null(b_vc_s0) ==FALSE ){### the same for b_vc_t0 and b_vc_tc0 (2023/4)
    bt_vc_s0		    <- as.matrix( b_vc_s0 / bse_vc_s0 )
    bp_vc_s0        <- 2 - 2 * pt( abs( bt_vc_s0 ), df = n - df )

    b_vc_s0		      <- data.frame( b_vc_s0 )
    bse_vc_s0       <- data.frame( bse_vc_s0 )
    bt_vc_s0		    <- data.frame( bt_vc_s0 )
    #bp_vc_s0		    <- data.frame( bp_vc_s0 )
    bp_vc_s0    <- data.frame( apply(bp_vc_s0, 2, function(x) p.adjust(x, method="BY" ) ) )

    names( b_vc_s0 )	   <- c( "(Intercept)", xname )
    names( bse_vc_s0 )   <- c( "(Intercept)", xname )
    names( bt_vc_s0 )	   <- c( "(Intercept)", xname )
    names( bp_vc_s0 )	   <- c( "(Intercept)", xname )

    B_vc_s <-list(b_vc_s0, bse_vc_s0, bt_vc_s0, bp_vc_s0)
  }
  if(is.null(b_vc_n0) ==FALSE ){
    bt_vc_n0		    <- as.matrix( b_vc_n0 / bse_vc_n0 )
    bp_vc_n0        <- 2 - 2 * pt( abs( bt_vc_n0 ), df = n - df )

    b_vc_n0		      <- data.frame( b_vc_n0 )
    bse_vc_n0       <- data.frame( bse_vc_n0 )
    bt_vc_n0		    <- data.frame( bt_vc_n0 )
    #bp_vc_n0		    <- data.frame( bp_vc_n0 )
    bp_vc_n0        <- data.frame( apply(bp_vc_n0, 2, function(x) p.adjust(x, method="BY" ) ) )

    names( b_vc_n0 )	<- c( "(Intercept)", xname )
    names( bse_vc_n0 )<- c( "(Intercept)", xname )
    names( bt_vc_n0 )	<- c( "(Intercept)", xname )
    names( bp_vc_n0 )	<- c( "(Intercept)", xname )

    B_vc_n <-list(b_vc_n0, bse_vc_n0, bt_vc_n0, bp_vc_n0)
  }

  if( nnxf != 0 ) {
    bft_vc		<- as.matrix( bf_vc / bfse_vc )
    bfp_vc		<- 2 - 2 * pt( abs( bft_vc ), df = n - df )
    bf_vc		<- data.frame( bf_vc )
    bfse_vc		<- data.frame( bfse_vc )
    bft_vc		<- data.frame( bft_vc )
    bfp_vc		<- data.frame( bfp_vc )
    names( bf_vc )	  <- xxfname
    names( bfse_vc )	<- xxfname
    names( bft_vc )	<- xxfname
    names( bfp_vc )	<- xxfname
  } else {
    bf_vc    <- NULL
    bfse_vc  <- NULL
    bft_vc   <- NULL
    bfp_vc   <- NULL
  }

  parV		        <- parV * sqrt( sig )
  sf_par		      <- data.frame( rbind( parV, moran_vc ) )
  if( is.null( parV_t )&is.null( parV_tc ) ){
    sf_name         <- c( "random SD: Spatial", "Moran.I/max(Moran.I)" )
  } else {
    sf_name         <- c( "random SD: Spatial", "--Moran.I/max(Moran.I)" )
  }

  if( !is.null( parV_t ) ){
    sf_par    <- rbind(sf_par, parV_t * sqrt( sig ))
    sf_name   <- c(sf_name, "random SD: coords_z[,1]")
  }
  if( !is.null( parV_tc ) ){
    sf_par   <- rbind(sf_par, parV_tc * sqrt( sig ))
    sf_name<- c(sf_name, "random SD: coords_z[,2]")
  }
  names( sf_par )	  <- c( "(Intercept)", xname )
  rownames( sf_par )<- sf_name

  s_g             <- NULL
  if( ng != 0 ){
    parG		<- t( parG * sqrt( sig ) )
    bg_par	<- data.frame( parG )
    rownames( bg_par )	<- "random SD"
    names( bg_par )	<- names( as.data.frame(xgroup) )
    s_g     <- bg_par
  }

  vc		          <- data.frame(ifelse( sf_par[1,] ==0, 0, 1) )
  names( vc )	    <- names( sf_par )
  rownames( vc )	<- "Spatial"

  vc_row          <- 2
  if( !is.null( parV_t ) ){
    vc  <- rbind( vc, ifelse( parV_t == 0, 0, 1 ))
    rownames( vc )[vc_row]	<- "coords_z[,1]"
    vc_row        <-vc_row + 1
  }
  if( !is.null( parV_tc ) ){
    vc  <- rbind( vc, ifelse( parV_tc == 0, 0, 1 ))
    rownames( vc )[vc_row]	<- "coords_z[,2]"
    vc_row        <-vc_row + 1
  }

  if(!is.null(id_ntv_interact)){
    vc     <- rbind( vc, as.integer(id_ntv_interact%in%unique(idd)))
    rownames( vc )[vc_row]	<- "Spatial x coords_z[,1]"
    vc_row        <-vc_row + 1

    sf_par <- rbind(sf_par,0)
    sf_name<- c(sf_name, "random SD: Spatial x coords_z[,1]")
    rownames( sf_par )<- sf_name
    if( sum( vc[nrow(vc),] ) > 0 ){
      for(vv in which(par0_int_dat$int_sel==1)){
        sf_par[nrow(sf_par),par0_int_dat$x_sel[vv]]<-par0_int_dat$random_SE[vv]
      }
    }
  }

  if(!is.null(id_ntcv_interact)){
    vc     <- rbind( vc, as.integer(id_ntcv_interact%in%unique(idd)))
    rownames( vc )[vc_row]	<- "Spatial x coords_z[,2]"
    vc_row        <-vc_row + 1

    sf_par <- rbind(sf_par,0)
    sf_name<- c(sf_name, "random SD: Spatial x coords_z[,2]")
    rownames( sf_par )<- sf_name
    if( sum( vc[nrow(vc),] ) > 0 ){
      for(vv in which(par0_int_dat$int_sel==2)){
        sf_par[nrow(sf_par),par0_int_dat$x_sel[vv]]<-par0_int_dat$random_SE[vv]
      }
    }
  }

  if( !is.null( parNsv ) ){
    vc  <- rbind( vc, c(0, ifelse( parNsv ==0, 0, 1 )))
    rownames( vc )[vc_row]	<- "Non-spatial/Non-linear"
    vc_row        <-vc_row + 1

    parN		<- t( parNsv * sqrt( sig ) )
    bn_par	<- data.frame( parN )
    rownames( bn_par )	<- "random SD"
    names( bn_par )	<- xxname
    s_n    <- bn_par

    totv   <-c(t(sf_par[1,]+c(0,s_n)))#as.vector(sf_par[1,]+c(0,s_n))
    s_rat  <-c(c(t(sf_par[1,]))/totv)#as.vector(sf_par[1,])/totv
    n_rat  <-c(c(0,t(s_n))/totv)#as.vector(c(0,s_n))/totv
    s_rat[totv==0]<-NA;n_rat[totv==0]<-NA
    sn_rat <-data.frame(rbind(s_rat,n_rat))
    rownames(sn_rat)<-c("Share (Spatial)","Share (Non-spatial/Non-linear)")
    names(sn_rat)[1]<-"(Intercept)"
    s       <- list(sf_par, s_n, sn_rat)### Temporal is not implemented (2023/4)
  } else {
    s       <- sf_par
  }

  s_nxconst<- NULL
  if( is.null( parNxf ) ==FALSE ){
    vc_n      <- as.data.frame( t(ifelse( parNxf ==0, 0, 1 )))
    names( vc_n )	  <- xxfname
    rownames( vc_n )<- "Non-spatial"
    vc        <- list( vc, vc_n )
    parNf		  <- t( parNxf * sqrt( sig ) )
    bnf_par		<- data.frame( parNf )
    rownames( bnf_par )	<- "random SE"
    names( bnf_par )	  <- xxfname
    s_nxconst<- bnf_par
  }



  ################ pred_quantile
  pquant        <- c(0.01, 0.025, 0.05, seq(0.1,0.9,0.1), 0.95, 0.975, 0.99)
  pquant2       <- seq(0.001,0.999,0.001)

  pq_dat0       <- NULL
  for(pq in pquant){
    pq_dat0     <- cbind(pq_dat0,qnorm(pq,pred0,pred0_se))
  }
  pq_dat0       <- as.data.frame(pq_dat0)
  names(pq_dat0)<- paste("q",pquant,sep="")

  if( y_nonneg==TRUE ){
    tr_bpar          <- data.frame(tr_bpar[1])
    names(tr_bpar)   <- "Estimates"
    rownames(tr_bpar)<- "Lambda (Box-Cox)"
  } else {
    tr_bpar          <- NULL
  }

  if( tr_num > 0 ){
    z0       <- sal_k(par=tr_par,y=y0,k=tr_num,noconst_last=noconst_last,
                      bc_par=tr_bpar$Estimates,jackup=jackup,y_added2=y_added2)#
    z_ms     <- z0$z_ms
    y_ms     <- z0$y_ms
    dif      <- 1/d_sal_k(par=tr_par,y=y0,k=tr_num,noconst_last=noconst_last,
                          bc_par=tr_bpar$Estimates,y_ms=y_ms,z_ms=z_ms,jackup=jackup)
    pred     <- i_sal_k(par=tr_par,y=pred0,k=tr_num,noconst_last=noconst_last,
                        bc_par=tr_bpar$Estimates,y_ms=y_ms,z_ms=z_ms,jackup=jackup,y_added2=y_added2)
    # - y_added,ulim=ulim

    pq_dat   <- pq_dat0
    for(pq in 1:ncol(pq_dat0)){
      ptest<-try(pq_pred<- i_sal_k(par=tr_par,y=pq_dat0[,pq],k=tr_num,noconst_last=noconst_last,
                                   bc_par=tr_bpar$Estimates,y_ms=y_ms,z_ms=z_ms,jackup=jackup,
                                   y_added2=y_added2))#,  - y_added, ulim=ulim
      if( !inherits(ptest, "try-error") ){#class(ptest)!="try-error"
        pq_dat[,pq]     <- pq_pred
      } else {
        pq_dat[,pq]     <- NA
      }
    }

    min_s        <- min(y) - (max(y)-min(y))/10
    max_s        <- max(y) + (max(y)-min(y))/10
    if( y_type == "continuous"){
      snorm       <- seq(min_s, max_s, len=100000)
      snorm_d     <- dnorm( snorm )
      snorm_tr    <- i_sal_k(par=tr_par,y=snorm,k=tr_num,noconst_last=noconst_last,
                             bc_par=tr_bpar$Estimates,y_ms=y_ms,z_ms=z_ms,jackup=jackup) - y_added#, ulim=ulim,y_added2=y_added2
      y_probfun00 <- cbind(snorm_tr,snorm_d)
      y_probfun00 <- y_probfun00[is.finite(y_probfun00[,1]),]
      y_prob_sk      <- y_probfun00
      y_bin          <- seq(0.8*min(y0), 1.2*max(y0), len=1000)
      suppressWarnings(approx_fun     <- approxfun(x=y_probfun00[,1],y=y_probfun00[,2], rule=2))
      y_probfun      <- cbind( y_bin, approx_fun(y_bin))

      pred        <- data.frame(pred = pred, pred_transG = pred0, pred_transG_se = pred0_se)

    } else if( y_type =="count" ){
      snorm       <- seq(min_s, max_s, len=1000000)
      snorm_d     <- dnorm( snorm )
      snorm_d     <- c(pnorm( snorm )[1],diff(pnorm( snorm )))

      snorm_tr    <- i_sal_k(par=tr_par,y=snorm,k=tr_num,noconst_last=noconst_last,
                             bc_par=tr_bpar$Estimates,y_ms=y_ms,z_ms=z_ms,jackup=jackup)# - y_added#,y_added2=y_added2
      if( !is.null(offset)){
        snorm_tr_pred <- median( offset )*exp( snorm_tr )
        snorm_tr      <- snorm_tr_pred * exp(- (1+0.5*zrat)/( snorm_tr_pred + 0.5 ) )
      } else {
        snorm_tr_pred <- exp( snorm_tr )
        snorm_tr      <- snorm_tr_pred * exp(- (1+0.5*zrat)/( snorm_tr_pred + 0.5) )
      }
      y_probfun00 <- cbind( round( snorm_tr ), snorm_d )
      y_probfun00 <- y_probfun00[is.finite(y_probfun00[,1]),]

      y_bin          <- min( y_probfun00[,1] ):min(max(y_org*1.5),max( y_probfun00[,1] ))
      y_probfun00_id <- get.knnx(y_bin,y_probfun00[,1],k=1)$nn.index
      y_bin2         <- y_bin[y_probfun00_id]
      y_probfun      <- aggregate(y_probfun00[,2],list(y_bin2),sum)
      y_probfun[,2]  <- y_probfun[,2]/sum(y_probfun[,2])
      y_prob_sk      <- y_probfun
      y_probfun      <- y_probfun[( y_probfun[,1] >= 0.8*min(y_org) )&( y_probfun[,1] <= 1.2*max(y_org) ),]

      dif            <- dif*y_org
      pq_dat         <- exp(pq_dat)
      pred           <- data.frame(pred = exp(pred), pred_transG = pred0)
    }

  } else if( y_nonneg ==TRUE ){
    y        <- bc(par=tr_bpar$Estimates,y=y0,jackup=jackup)
    z_ms     <- NULL
    y_ms    <- c(mean(y),sd(y))

    dif      <- 1 / ( d_bc(par=tr_bpar$Estimates,y=y0,jackup=jackup)*d_sc(par=y_ms,y) )### changed 2022/08
    pred0    <- pred0*y_ms[2] + y_ms[1]
    pred     <- i_bc(par=tr_bpar$Estimates,y=pred0,jackup=jackup) - y_added# - y_added2,ulim=ulim

    pq_dat   <- pq_dat0
    for(pq in 1:ncol(pq_dat0)){
      ptest  <- try(suppressWarnings(pq_pred<- i_bc(par=tr_bpar$Estimates,
                    y=pq_dat0[,pq]*y_ms[2] + y_ms[1],jackup=jackup) - y_added))#,ulim=ulim
      if( !inherits(ptest, "try-error")){#class(ptest)!="try-error"
        pq_dat[,pq]       <-pq_pred
      } else {
        pq_dat[,pq]       <-NA
      }
    }

    min_s     <- min(y) - (max(y)-min(y))/10
    max_s     <- max(y) + (max(y)-min(y))/10
    snorm     <- seq(min_s, max_s, len=10000)
    snorm_d   <- dnorm(snorm,mean=mean(y),sd=sd(y))
    tr_bpar00 <- tr_bpar
    snorm_tr     <- i_bc(par=tr_bpar00$Estimates,y=snorm,jackup=jackup)#exp(  - y_added
    y_probfun00  <- cbind(snorm_tr,snorm_d)
    y_probfun00  <- y_probfun00[is.finite(y_probfun00[,1]),]
    y_prob_sk    <- y_probfun00

    #suppressWarnings(approx_fun<- approxfun(x=y_probfun00[,1],y=y_probfun00[,2], rule=2))
    #y_bin     <- seq(max(0, 0.8*min(y0)), 1.2*max(y0), len=1000)
    #y_probfun <- cbind( y_bin, approx_fun(y_bin))

    y_bin          <- seq(0.8*min(y0), 1.2*max(y0), len=1000)
    suppressWarnings(approx_fun     <- approxfun(x=y_probfun00[,1],y=y_probfun00[,2], rule=2))
    y_probfun      <- cbind( y_bin, approx_fun(y_bin))

    pred        <- data.frame(pred = pred, pred_transG = pred0, pred_transG_se = pred0_se)

  } else {
    dif          <- 1
    z_ms         <- NULL
    y_ms         <- NULL
    tr_par       <- NULL
    tr_bpar      <- NULL
    min_s        <- min(y) - (max(y)-min(y))/10
    max_s        <- max(y) + (max(y)-min(y))/10
    if( y_type == "continuous" ){
      pq_dat   <- pq_dat0
      pred     <- data.frame(pred = pred0, pred_se=pred0_se)
      snorm    <- seq(min_s, max_s, len=10000)
    } else if( y_type == "count" ){
      dif      <- dif*y_org
      pq_dat   <- exp( pq_dat0 )
      pred     <- data.frame(pred =exp( pred0 ), pred_transG = pred0 )
      snorm    <- seq(min_s, max_s, len=1000000)
    }

    snorm_d  <- dnorm(snorm,mean=mean(pred0),sd=sd(pred0))
    if( y_type=="count" ){

      if( !is.null(offset)){
        snorm_tr_pred <- median( offset )*exp( snorm )
        snorm_tr      <- snorm_tr_pred * exp(- (1+0.5*zrat)/( snorm_tr_pred + 0.5 ) )
      } else {
        snorm_tr_pred <- exp( snorm )
        snorm_tr      <- snorm_tr_pred * exp(- (1+0.5*zrat)/( snorm_tr_pred + 0.5) )
      }

      y_probfun   <- cbind( round( snorm_tr ), snorm_d )
      y_probfun    <- cbind(snorm_tr,snorm_d)
      y_probfun00  <- y_probfun[is.finite(y_probfun[,1]),]
      y_bin        <- min( y_probfun00[,1] ):max( y_probfun00[,1] )
      y_probfun00_id<-get.knnx(y_bin,y_probfun00[,1],k=1)$nn.index
      y_bin2       <- y_bin[y_probfun00_id]
      y_probfun    <- aggregate(y_probfun00[,2],list(y_bin2),sum)
      y_prob_sk    <- y_probfun

      y_probfun[,2]<- y_probfun[,2]/sum(y_probfun[,2])
      y_probfun    <- y_probfun[( y_probfun[,1] >= 0.8*min(y_org) )&( y_probfun[,1] <= 1.2*max(y_org) ),]

    } else {
      y_probfun00  <- cbind(snorm,snorm_d)
      y_probfun00  <- y_probfun00[is.finite(y_probfun00[,1]),]
      y_prob_sk    <- y_probfun00

      suppressWarnings(approx_fun   <- approxfun(x=y_probfun00[,1],y=y_probfun00[,2], rule=2))
      y_bin        <- seq(0.8*min(y0), 1.2*max(y0), len=1000)
      y_probfun    <- cbind( y_bin, approx_fun(y_bin))
    }
  }

  y_probfun        <- data.frame(y_probfun)
  names(y_probfun) <- c("Value","Intensity")

  cumrat<- cumsum(y_prob_sk[,2])/sum(y_prob_sk[,2])
  m1_est<- weighted.mean(x=y_prob_sk[,1],w=y_prob_sk[,2])
  m2_est<- weighted.mean(x=(y_prob_sk[,1] - m1_est)^2,w=y_prob_sk[,2])
  m3_est<- weighted.mean(x=(y_prob_sk[,1] - m1_est)^3,w=y_prob_sk[,2])
  m4_est<- weighted.mean(x=(y_prob_sk[,1] - m1_est)^4,w=y_prob_sk[,2])
  if( (( tr_num > 0 )|( y_nonneg ==TRUE ))&( y_type == "continuous" ) ){
    Skew  <- m3_est/m2_est^1.5
    Kurt  <- m4_est/m2_est^2 - 3
    skew_kurt<- data.frame(Estimates = c(Skew, Kurt))
    row.names( skew_kurt )<-c("skewness","excess kurtosis")
  } else if(y_type == "count"){
    Skew      <- m3_est/m2_est^1.5
    Kurt      <- m4_est/m2_est^2 - 3
    skew_kurt<- data.frame(Estimates= c(Skew, Kurt))
    row.names( skew_kurt )<-c("skewness","excess kurtosis")
  } else {
    Skew  <- 0
    Kurt  <- 0
    skew_kurt<- data.frame(Estimates= c(Skew, Kurt))
    row.names( skew_kurt )<-c("skewness","excess kurtosis")
  }

  if( y_type=="count" ){
    if( !is.null( offset )) {
      pred[,1] <-pred[,1]*offset
      pq_dat   <-round(pq_dat*offset)
    } else {
      pq_dat   <-round(pq_dat)
    }
  }

  if( y_type=="continuous" ){
    e_stat		 <- data.frame( stat = c( sqrt( sig_org ), r2, loglik, AIC, BIC ) )

    if( y_nonneg ==FALSE & tr_num ==0 ){
      rownames( e_stat ) <- c( "resid SE", "adjR2(cond)", lik_nam, "AIC", "BIC")
    } else {
      rownames( e_stat ) <- c( "resid SE(trans)", "adjR2(cond)", lik_nam, "AIC", "BIC")
    }

    ######################## Approximate inference for the NULL model
    if( length(weight)==1 ){
      mod_NULL_id0<-" )"
      weight      <- NULL#rep(1,n)
    } else {
      mod_NULL_id0<-", weights = weight )"
    }

    if( is.null( x ) & is.null( xconst ) ){
      mod_NULL   <- lm(y0 ~ 1, weights = weight )
      mod_NULL_id<- paste("lm( y ~ 1", mod_NULL_id0, sep="")
    } else if( is.null( x ) ){
      lm_dat     <- data.frame(y0, xconst)
      mod_NULL   <- lm(y0 ~ ., data=lm_dat, weights = weight )
      if( resf_flag ){
        mod_NULL_id<- paste("lm( y ~ x", mod_NULL_id0, sep="")
      } else {
        mod_NULL_id<- paste("lm( y ~ xconst", mod_NULL_id0, sep="")
      }
    } else if( is.null( xconst ) ){
      lm_dat     <- data.frame(y0, x)
      mod_NULL   <- lm(y0 ~ ., data=lm_dat, weights = weight )
      mod_NULL_id<- paste("lm( y ~ x", mod_NULL_id0, sep="")
    } else {
      lm_dat     <- data.frame(y0, x, xconst)
      mod_NULL   <- lm(y0 ~ ., data = lm_dat, weights = weight )
      mod_NULL_id<- paste("lm( y ~ x + xconst", mod_NULL_id0, sep="")
    }

    e_stat_N     <- data.frame( stat = c( logLik(mod_NULL), AIC(mod_NULL), BIC(mod_NULL) ) )
    rownames( e_stat_N ) <- c( lik_nam, "AIC", "BIC")
    e_stat_NULL  <- list(NULL)
    e_stat_NULL[[1]]<- e_stat_N
    e_stat_NULL[[2]]<- mod_NULL_id
    r2_devrat<-r2

  } else {
    lik_nam2  <- paste( "Gaussian ",lik_nam, " approximating the model",sep="")
    #if( tr_num==0 & y_nonneg==FALSE ){
    #  e_stat		<- data.frame( stat = c( sig_org, rmse, 100*dev_rat, loglik, AIC, BIC ) )
    #  rownames( e_stat ) <- c( "dispersion parameter","RMSE","deviance explained (%)", lik_nam2, "AIC", "BIC")
    #} else {
    e_stat		<- data.frame( stat = c( rmse, 100*dev_rat, loglik, AIC, BIC ) )
    rownames( e_stat ) <- c("RMSE","deviance explained (%)", lik_nam2, "AIC", "BIC")
    #}

    ######################## Approximate Gaussian likelihood for the NULL model
    if( is.null(offset) ){
      y_org2   <- log( y_org + 0.5) + y_added2
      mod_NULL_id0<-", family = poisson )"
    } else {
      y_org2<- log( ( y_org+0.5 )/offset ) +y_added2
      mod_NULL_id0<-", offset = log( offset ), family = poisson )"
    }

    if( length( weight0 ) > 1 ){
      mod_NULL_id0<- paste( ", weights = weight", mod_NULL_id0, sep="" )
    }

    if( is.null( x ) & is.null( xconst ) ){
      mod_NULL   <- lm(y_org2 ~ 1, weights = weight_lik )
      mod_NULL_id<- paste("glm( y ~ 1", mod_NULL_id0, sep="")
    } else if( is.null( x ) ){
      lm_dat     <- data.frame(y_org2, xconst)
      mod_NULL   <- lm(y_org2 ~ ., data=lm_dat, weights = weight_lik )
      #mod_NULL   <- lm(y_org2 ~ as.matrix( xconst ), weights = weight_lik )
      if( resf_flag ){
        mod_NULL_id<- paste("glm( y ~ x", mod_NULL_id0, sep="")
      } else {
        mod_NULL_id<- paste("glm( y ~ xconst", mod_NULL_id0, sep="")
      }
    } else if( is.null( xconst ) ){
      lm_dat     <- data.frame(y_org2, x)
      mod_NULL   <- lm(y_org2 ~ ., data=lm_dat, weights = weight_lik )
      mod_NULL_id<- paste("glm( y ~ x", mod_NULL_id0, sep="")
    } else {
      lm_dat     <- data.frame(y_org2, x, xconst)
      mod_NULL   <- lm(y_org2 ~ ., data=lm_dat, weights = weight_lik )
      mod_NULL_id<- paste("glm( y ~ x + xconst", mod_NULL_id0, sep="")
    }

    e_stat_N     <- data.frame( stat = c( logLik(mod_NULL), AIC(mod_NULL), BIC(mod_NULL) ) )
    lik_nam2_NULL<- paste("Gaussian ",lik_nam, " approximating the NULL model",sep="")
    rownames( e_stat_N ) <- c( lik_nam2_NULL, "AIC", "BIC")
    e_stat_NULL  <- list(NULL)
    e_stat_NULL[[1]]<- e_stat_N
    e_stat_NULL[[2]]<- mod_NULL_id

    r2_devrat<-dev_rat0
  }

  messs   <- 0
  #if( sum( n_omit ) > 5 ) {
  #  if(y_type=="count" & r2_devrat <= -0.1 ){
  #    message( "Note: Singular fit. Simplify the model")
  #  } else {
  #    if( r2_devrat <= -0.1 ){
  #      #  message( "Note: The model is nearly singular. Consider simplifying the model" )
  #      #} else {
  #      message( "Note: Singular fit. Simplify the model")
  #    }
  #    messs <- 1
  #  }
  #
  #} else
  if( y_type=="continuous" & r2_devrat < -0.1 ){
    message("Note: Singular fit. Simplify the model")
    messs <- 1
  } else if( y_type=="count" & r2_devrat <= -0.1 ){
    message("Note: Singular fit. Simplify the model")
    messs <- 1
  }

  #if( (messs == 0)&( loglik < logLik( mod_NULL )) ){
  #  message( "Note: The model is nearly singular. Consider simplifying the model")
  #}


  other		<- list( res_int =res_int, r = r, sf_alpha = parR, x_id = x_id, nx=nx, nxf = nxf, xf_id = xf_id, df = df,
                  omit_list=omit_list,null_dum3=null_dum3, idd=idd, nev0=nev0,
                  b_s = bb, b_covs = bb_cov, B_covs = b_cov2, sig = sig, sig_org=sig_org ,xg_levels = xg_levels,
                  weight=weight,is_weight = length( weight )>1, penalty=penalty,
                  eevSqrt=eevSqrt, evSqrts = evSqrts, evSqrts_c = evSqrts_c, evSqrts_n = evSqrts_n,
                  evSqrts_t = evSqrts_t, evSqrts_tc = evSqrts_tc, evSqrts_t_int = evSqrts_t_int, evSqrts_tc_int = evSqrts_tc_int,
                  model = "resf_vc", b_c = bf_s, b_covs_c = bf_covs,
                  Bias=Bias, nvc_x=nvc_x, nvc_xconst=nvc_xconst, nvc_num = nvc_num, sel_basis_c = sel_basis_c, sel_basis_n = sel_basis_n,
                  x = x, xconst = xconst, coords = meig$other$coords, coords_z = meig$other$coords_z, dif=dif, y = y, y0 = y0,
                  tr_num=tr_num, y_nonneg = y_nonneg, y_type = y_type,
                  method=method, y_added = y_added, jackup=jackup, np=np, offset = offset, e_NULL = e_stat_NULL,
                  w_scale = w_scale,tr_comp=tr_comp,par0=par0,par0_int_dat=par0_int_dat,
                  b_sub=b_sub,interact_sel_all=ii_sel_all, int_ev_sel_list=ev_sel_list,
                  id_nsv=id_nsv, id_ntv=id_ntv, id_ntcv=id_ntcv, id_ntv_interact=id_ntv_interact,id_ntcv_interact=id_ntcv_interact,
                  Par=Par)

  result    <- list( b_vc = b_vc, bse_vc = bse_vc, t_vc = bt_vc, p_vc = bp_vc_adj, p_vc_naive=bp_vc, B_vc_s = B_vc_s, B_vc_n = B_vc_n,
                     c = b_par, c_vc = bf_vc, cse_vc = bfse_vc, ct_vc = bft_vc, cp_vc = bfp_vc, b_g = bpar_g2,
                     s = s, s_c=s_nxconst, s_g = s_g, vc = vc, e = e_stat, pred = pred, pred_quantile=pq_dat,
                     tr_par=tr_par,tr_bpar=tr_bpar,tr_y=y,
                     resid = resid, pdf=y_probfun, skew_kurt = skew_kurt, other = other,
                     call = match.call() )
  class( result) <- "resf_vc"
  return( result )
}

print.resf_vc <- function(x, ...)
{
  cat("Call:\n")
  print(x$call)
  sss  <-x$s
  if( is( x$s )[1] == "list" ) sss<-x$s[[1]]

  if( dim( sss )[ 1 ]==2 ){
    if( is( x$s )[1] != "list" ){
      cat("\n----Spatially varying coefficients on x (summary)----\n")
    } else {
      cat("\n----Spatially varying coefficients (S&NVC) on x (summary)----\n")
    }
  } else {
    if( is( x$s )[1] != "list" ){
      cat("\n----Spatio-temporally varying coefficients on x (summary)----\n")
    } else {
      cat("\n----Spatio-temporally varying coefficients (ST&NVC) on x (summary)----\n")
    }
  }

  cat("\nCoefficient estimates:\n")
  print( summary( x$b_vc ) )
  cat("\nStatistical significance (adjusted by Benjamini-Yekutieli method):\n")
  p01<-apply(x$p_vc,2,function(x) sum(x<0.01))
  p05<-apply(x$p_vc,2,function(x) sum(x<0.05)) - p01
  p10<-apply(x$p_vc,2,function(x) sum(x<0.10)) - p01 - p05
  p90<-length(x$p_vc[,1]) - p01 - p05 - p10
  pv <-data.frame(rbind( p90, p10, p05, p01))
  names(pv)[1]  <- "Intercept"
  row.names(pv) <- c("Not significant", "Significant (10% level)",
                     "Significant ( 5% level)","Significant ( 1% level)")
  print(pv)
  if( is.null(x$s_c) & !is.null(x$c) ){
    cat("\n----Constant coefficients on xconst----------------------------\n")
    print( x$c )

  } else if( !is.null(x$c_vc) ){
    cat("\n----Non-linear varying coefficients on xconst (summary)----\n")
    cat("\nCoefficient estimates:\n")
    print( summary( x$c_vc ) )
    cat("\nStatistical significance (adjusted by Benjamini-Yekutieli method):\n")
    cp01<-apply(x$cp_vc,2,function(x) sum(x<0.01))
    cp05<-apply(x$cp_vc,2,function(x) sum(x<0.05)) - cp01
    cp10<-apply(x$cp_vc,2,function(x) sum(x<0.10)) - cp01 - cp05
    cp90<-length(x$cp_vc[,1]) - cp01 - cp05 - cp10
    cpv <-data.frame(rbind( cp90, cp10, cp05, cp01))
    row.names(cpv) <- c("Not significant", "Significant (10% level)",
                        "Significant ( 5% level)","Significant ( 1% level)")
    print(cpv)
  }
  cat("\n----Variance parameters----------------------------------\n")
  if( dim(sss)[1]>2 ){
    cat("\nSpatio-temporal effects (coefficients on x):\n")
  } else {
    cat("\nSpatial effects (coefficients on x):\n")
  }
  print( sss )

  if( is( x$s )[1] == "list" ){
    cat("\nNon-linear effects wrt x (coefficients on x):\n")
    print(x$s[[2]])
  }
  if( !is.null(x$s_c) ){
    cat("\nNon-linear effects wrt xconst (coefficients on xconst):\n")
    print(x$s_c)
  }
  if( !is.null(x$s_g) ){
    cat("\nGroup effects:\n")
    print(x$s_g)
  }

  if( !is.null(x$skew_kurt)|!is.null(x$tr_bpar) ){
    cat("\n----Estimated probability distribution of y--------------\n")
    if( !is.null(x$skew_kurt) ) print(x$skew_kurt)
    if( !is.null(x$tr_bpar) & x$other$y_type =="continuous" ){
      cat( paste("(Box-Cox parameter: ", format(x$tr_bpar[1], digits=7),")\n",sep="") )
    } else if( x$other$y_type =="count" ){
      cat( paste("(dispersion parameter: ", format(x$other$sig_org, digits=7),")\n",sep="") )
    }
  }

  cat("\n----Error statistics-------------------------------------\n")
  print(x$e)

  loglik_NULL<- x$other$e_NULL[[1]][1,1]
  AIC_NULL   <- x$other$e_NULL[[1]][2,1]
  BIC_NULL   <- x$other$e_NULL[[1]][3,1]
  mod_NULL   <- x$other$e_NULL[[2]]
  if( x$other$y_type =="continuous" ){
    ml_name    <- ifelse( x$other$method=="reml", "(r)loglik: ", "loglik: " )
    cat( paste('\nNULL model: ', mod_NULL, sep="") )
    cat(paste("\n   ",ml_name,format(loglik_NULL,digits=7),sep=""))
    cat(paste(" ( AIC: ",
              format(AIC_NULL,digits=7), ",  BIC: ", format(BIC_NULL,digits=7)," )\n",sep=""))

  } else if( x$other$y_type =="count" ){
    ml_name0 <- ifelse( x$other$method=="reml", "(r)loglik", "loglik" )
    ml_name  <- paste("\n   Gaussian ",ml_name0, " approximating the model: ",sep="")
    cat( paste('\nNULL model: ', mod_NULL, sep="") )
    cat(paste(ml_name,format(loglik_NULL,digits=7), "\n",sep=""))
    cat(paste("   ( AIC: ",
              format(AIC_NULL,digits=7), ",  BIC: ", format(BIC_NULL,digits=7)," )\n",sep=""))
  }

  if( (x$other$method=="reml")&(x$other$y_type=="continuous") ){
    cat('\nNote: AIC and BIC are based on the restricted/marginal likelihood.')
    cat('\n      Use method="ml" for comparison of models with different fixed effects (x and xconst)\n')
  }
  #if(x$other$y_type=="count"){
  #  cat('\nNote: Error statistics of the log-Gaussian model used to approximate')
  #  cat('\n      the Poisson model are available from other$e_stat_logG\n')
  #}
  invisible(x)
}

Try the spmoran package in your browser

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

spmoran documentation built on April 4, 2025, 12:20 a.m.