Nothing
besf_vc <- function( y, x, xconst = NULL, coords, s_id = NULL,
x_nvc = FALSE, xconst_nvc = FALSE,
x_sel = TRUE, x_nvc_sel = TRUE, xconst_nvc_sel = TRUE, nvc_num=5,
method = "reml", penalty = "bic", maxiter = 30, tol=1e-30,
covmodel="exp",enum = 200, bsize = 4000, ncores=NULL ){
n <- length( y )
nvc_x <- x_nvc
nvc_xconst <- xconst_nvc
xgroup <- NULL
ng <- 0
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( sel_xgroup ) ){
# allg<- !sel_xgroup
#} else {
# allg<- sel_xgroup
#}
if( x_nvc[1] == FALSE ) x_nvc_sel = FALSE
if( xconst_nvc[1] == FALSE ) xconst_nvc_sel = FALSE
sel_basis_vif <-function(test0, vif_max=15){
tres3 <-TRUE
rem_ord <-NULL
while(tres3){
suppressWarnings(cor_test0<-cor(test0))
testt2<-try(vif <- diag(solve(cor_test0, tol=tol)), silent=TRUE)
#testt2<-try(vif <- diag(solve(cor(test0))), silent=TRUE)
if( inherits(testt2, "try-error") ){#class(testt2)=="try-error"
suppressWarnings(corr_0<-cor(test0))
#corr_0 <- cor(test0)
diag(corr_0)<- 0
cor_max <- apply(corr_0,2,max)
rem_ord0 <- which(cor_max == max(cor_max))[1]
test0 <- test0[, -rem_ord0]
rem_ord <- c( rem_ord, rem_ord0 )
} else {
if( max(vif)>vif_max ){
rem_ord0<- which( vif == max( vif ))[ 1 ]
test0 <- test0[ ,-rem_ord0 ]
} else {
tres3<-FALSE
}
}
}
return(test0)#list(selvar=test0, rem_ord=rem_ord)
}
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 ]
term3_0 <- Msub_00 - Msub_01 %*% solve( 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
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)
}
lik_resf_vc <- function( par0, par0_est, par0_id, par0_sel, ev, M, M0inv, M0inv_01, M0inv_00,
m, yy, b_01, b_02, n, nx, nsv, nnxf, nnsv, ng, emet, term2, term3_0, null_dum2, id ){
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 ] ){
evv <- ev ^ par[ nsv + i ] * sum( ev ) / sum( ev ^ par[ nsv + i ] )
evSqrt <- c( evSqrt, par[ i ] * sqrt( evv ) )
}
if( ng != 0 ){
for( j in (1:ng)[ null_dum2[ (( nsv + 1 ):( nsv + ng )) ] == 0 ] ){#null_dum2[ -(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_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( 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 )
}
}
Mdet <- Mdet_f( id = id, par0_sel=par0_sel, term2 = term2, term3_0 = term3_0, evSqrt = evSqrt )
M2 <- M
for( j in 1:max( id ) ){
diag( M[ id == j, id == j ])<-diag( M[ id == j, id == j ] ) + 1/evSqrt[ id[ -c( 1:nx ) ] == j ] ^ 2
}
diag(M0inv_00) <- diag( M0inv_00 ) + evSqrt[ id[ id != 0 ] == par0_sel ] ^ 2
b_02_b <- solve( M0inv_00, tol=tol ) %*% b_02
b_02 <- M0inv_01 %*% b_02_b
b <- b_01 - b_02
sse <- yy - 2 * t( b ) %*% m + t( b ) %*% M2 %*% b
dd <- sse + sum( ( b[ -( 1:nx ) ] / evSqrt ) ^ 2 )
if( emet == "reml" ){
loglik <- Mdet + ( n - nx ) * ( 1 + log( 2 * pi * dd / ( n - nx ) ) )
} else if( emet == "ml" ){
loglik <- Mdet + n * ( 1 + log( 2 * pi * dd / n ) )
}
return( loglik )
}
lik_resf_vc0 <- function( par0, ev, M, m, yy, n, nx, nsv, ng, nnxf, nnsv, emet, null_dum4 ){
par <- par0 ^ 2
evSqrt <- NULL
for( i in ( 1:nsv )[ null_dum4[ 1:nsv ] == 0 ] ){
evv <- ev ^ par[ nsv + i ] * sum( ev ) / sum( ev ^ 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 )
}
}
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
m[-(1:nx)] <- m[ -( 1:nx ) ] * evSqrt
test <-try(Minv <- solve( M, tol=tol ))
#if(class(test)[1]=="try-error"){
if( inherits(test,"try-error") ){
loglik <- Inf
} else {
b <- Minv %*% m
sse <- yy - 2 * t( b ) %*% m + t( b ) %*% M0 %*% b
dd <- 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
}
return( loglik[ 1 ] )
}
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( 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
}
#XX1_0 <- list(NULL)
#XX1 <- NULL
#nvc_x_use <- apply( X1, 2, function( x ) length( unique( x ))) >= 3
if( is.logical( nvc_x[ 1 ] )&( nvc_x[ 1 ] == TRUE) ) nvc_x <- 1:nx0
pmod_x <- list(NULL)
nsamp <- 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 )#chganged
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)]
#if( n > 2000000 ){############################# to be implemented #################
# if( !is.null(nsamp) ) nsamp <- sample( 2000000 )
# testt<-try(XX1_00_b<- ns( X1_nvc[ nsamp,ii ], knots=knots ), silent=TRUE)
#} else {
testt<-try(XX1_00<- ns( X1_nvc[ ,ii ], knots = knots ), silent=TRUE)
#}
test <- inherits(testt, "try-error")#class(testt)[1] == "try-error"
iiii <- iiii+1
}
#if( n > 2000000 ){
# testt2<-try(XX1_00 <- sel_basis_vif( cbind( X1_nvc[nsamp,ii], XX1_00_b), vif_max=15),silent=TRUE)
#} else{
testt2<-try(XX1_00 <- sel_basis_vif( cbind( X1_nvc[,ii], XX1_00), vif_max=15),silent=TRUE )
#}
if( !inherits( testt2, "try-error" ) ){#class(testt2)[1] != "try-error"
np_xx[ ii ] <- dim(XX1_00)[2]
} else {
np_xx[ ii ] <- 0
}
if( np_xx[ ii ] > 2 ){
pmod_x[[ii]]<- scale( XX1_00 )
} else {
pmod_x[[ii]]<- 0
np_xx[ ii ] <- 0
}
}
} else {
nnsv <- 0
np_xx <- NULL
}
if( is.logical( nvc_xconst[ 1 ] )&( nvc_xconst[ 1 ] == TRUE) ) nvc_xconst<- 1:nxf
pmod_xconst <- 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 )
rem_ord_XXc<-list(NULL)
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)]
#if( n > 2000000 ){############################# to be implemented #################
# if( !is.null(nsamp) ) nsamp <- sample( 2000000 )
# testt<-try(XXconst_00_b<- ns( Xconst_nvc[ nsamp,ii ], knots = knots ), silent=TRUE)
#} else {
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
}
#if( n > 2000000 ){
# testt2<-try(XXconst_00 <- sel_basis_vif( cbind( Xconst_nvc[nsamp,ii], XXconst_00_b ), vif_max=15),silent=TRUE)
#} else{
testt2<-try(XXconst_00 <- sel_basis_vif( cbind( Xconst_nvc[,ii], XXconst_00 ), vif_max=15),silent=TRUE)
#}
if( !inherits(testt2, "try-error") ){#class(testt2)[1] != "try-error"
np_xxconst[ ii ] <- dim(XXconst_00)[2]
}else{
np_xxconst[ ii ] <- 0
}
if( np_xxconst[ ii ] > 2 ){
pmod_xconst[[ii]]<- scale( XXconst_00 )
} else {
pmod_xconst[[ii]]<- 0
np_xxconst[ ii ] <- 0
}
}
} else {
nnxf <- 0
np_xxconst <- NULL
}
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" )
}
}
#X0 <- as.matrix( cbind( 1, Xconst, X1 ) )
#Mo <- crossprod( X0 )
#mo <- crossprod( X0, y )
#parVmax_sq <- sqrt( sd( y ) / sd( y - X0 %*% ( solve( Mo ) %*% mo ) ) )
nsv <- ifelse( is.null( X1 ), 1, dim( X1 )[ 2 ] + 1 )
if(nsv>=7){
message("Note: At most 5 covariates in x are recommended for stable estimation.")
message(" Consider using xconst")
}
if( !is.null( s_id )[ 1 ] ){
coords_gx<-tapply(coords[,1],s_id,mean)
coords_gy<-tapply(coords[,2],s_id,mean)
coords_g <-as.matrix(cbind(coords_gx, coords_gy))
#s_id2 <-data.frame( s_id=rownames(coords_x), s_id_num = 1:length(coords_x) )
#s_id_dat<-data.frame( s_id = s_id, id = 1:length(s_id))
#s_id_dat2<-merge(s_id_dat, s_id2, by="s_id", all.x=TRUE)
#s_id_dat2<-s_id_dat2[order(s_id_dat2$id),c("s_id", "s_id_num")]
enum <- min( enum, length( coords_g ) - 1 )
suppressWarnings(coordk <- kmeans( coords, centers = enum + 1 )$centers)
} else {
suppressWarnings(coordk <- kmeans( coords, centers = enum + 1 )$centers)
}
D <- rdist( coordk )
h <- max( spantree( D )$dist )
if( covmodel == "exp" ){
C <- exp( -D / h )
} else if( covmodel == "gau" ){
C <- exp( -( D / h) ^ 2 )
} else if( covmodel == "sph" ){
C <- ifelse( D < h , 1 - 1.5 * (D / h ) + 0.5 * ( D / h ) ^ 3, 0 )
} else {
stop( "covmodel is not specified appropriately" )
}
Cmean <- apply( C, 1, mean )
MCM <- t( C - Cmean ) - Cmean + mean( Cmean )
eigenC <- eigen( MCM )
sf0 <- eigenC$vectors[ , eigenC$values > 1e-08 ]
ev0 <- eigenC$values [ eigenC$values > 1e-08 ]
sf_ap0 <- sf0 %*% diag( 1 / ev0 )
ev_full <- ev0 * ( n / enum ) - 1
ev <- ev_full[ ev_full > 1e-08 ]
nev <- length( ev )
if( is.null( Xconst )&is.null( X1 )){
X0 <- as.matrix( rep(1,n) )
} else {
X0 <- as.matrix( cbind( 1, Xconst, X1 ) )
}
#X0 <- as.matrix( cbind( 1, Xconst, X1 ) )
nx <- dim( X0 )[ 2 ]
Mo <- crossprod( X0 )
mo <- crossprod( X0, y )
parVmax_sq <- sqrt( sd( y ) / sd( y - X0 %*% ( solve( Mo, tol=tol ) %*% mo ) ) )
nblock <-max(1,round((n*(nxf + length(ev)*nsv)/bsize^2)))
ids <- unique(c(round(seq(1,n,len=nblock+1)),n))
if(is.null(ncores)) {
ncores00<- detectCores()
ncores0 <- ifelse(ncores00<=2,ncores00, ifelse(ncores00<=5, ncores00-1,ncores00-2))
ncores <- makeCluster(ncores0,setup_strategy = "sequential")
} else {
ncores0<- ncores
ncores <- makeCluster(ncores0,setup_strategy = "sequential")
}
registerDoParallel(ncores)
message(paste0("Inner product evaluation is parallelized with ", ncores0," cores"))
im <- NULL
sfsf <- foreach(im = 1:(length(ids)-1), .packages=c("fields", "splines"),.combine="+") %dopar% {
if(im==1){
ids_sub<-ids[im]:ids[im+1]
} else {
ids_sub<-(ids[im]+1):ids[im+1]
}
if( covmodel == "exp" ){
ccc0 <- exp( -rdist( coordk, coords[ids_sub,] ) / h )
} else if( covmodel == "gau" ){
C <- exp( -( D / h) ^ 2 )
ccc0 <- exp( -( rdist( coordk, coords[ids_sub,] ) / h )^2 )
} else if( covmodel == "sph" ){
ddd0 <- rdist( coordk, coords[ids_sub,] )
ccc0 <- ifelse( ddd0 < h , 1 - 1.5 * ( ddd0 / h ) + 0.5 * ( ddd0 / h ) ^ 3, 0 )
}
sf <- t( ccc0 - Cmean ) %*% sf_ap0[ , ev_full > 1e-08 ]
X2 <- cbind(X0[ids_sub, ], sf)
cSum_sub0<- rep( 0, nx )
cSum_sub1 <- colSums( sf )
if( nsv >= 2 ){
for( im2 in 1:( nsv - 1 ) ){
X2 <- cbind( X2, X1[ids_sub, im2 ] * sf )
cSum_sub1 <- c( cSum_sub1, colSums( sf ) )
}
}
cSum_sub2 <- NULL
if( is.logical( nvc_xconst[ 1 ] ) == FALSE ){
for( ii in 1:length( np_xxconst )){
if(np_xxconst[ii]>0){
basis <- as.matrix( pmod_xconst[[ ii ]][ids_sub, ] )
X2 <- cbind( X2, Xconst_nvc[ids_sub, ii ]*basis )
cSum_sub2 <- c( cSum_sub2, colSums( basis ) )
}
}
}
cSum_sub3 <- NULL
if( is.logical( nvc_x[ 1 ] ) == FALSE ){
for( ii in 1:length( np_xx ) ){
if(np_xx[ii]>0){
basis <- as.matrix( pmod_x[[ii]][ids_sub, ] )
X2 <- cbind( X2, X1_nvc[ids_sub, ii ]*basis )
cSum_sub3 <- c( cSum_sub3, colSums( basis ) )
}
}
}
cSum <- c(cSum_sub0, cSum_sub1, cSum_sub2, cSum_sub3)
sfsf1<-cbind(cSum, crossprod( X2, y[ids_sub] ),crossprod( X2 ))
#sfsf1<-cbind(rep(0,length(sfsf0[,1])),sfsf0)
#clsf<-colSums(sf)
#sfsf1[1:length(clsf),1]<-clsf
sfsf1
}
sfmean <-sfsf[,1]/n
m0 <-c(sfsf[,2])
EX <-sfsf[, -(1:2)][ ,1:nx ]
id <- c( rep( 0, nx ), rep( 1:nsv, each = nev ) )
sf_id <- c(rep(1,nx),rep(1, nev))
sf_id2 <- 1
sf_id3 <- 1
if( nsv >=2 ){
sf_id <-c( sf_id, nxf + rep(2:nsv,each=nev) )
sf_id2 <- c(sf_id2, nxf + 2:nsv)
sf_id3 <- c(sf_id3, nxf + 2:nsv)
}
ii_count <- max( id )
if( is.logical( nvc_xconst[ 1 ] ) == FALSE ){
for( ii in 1:length( np_xxconst ) ){
ii_count <- ii_count + 1
if(np_xxconst[ii]>0){
id <- c( id, rep( ii_count, np_xxconst[ii] ))
sf_id <- c(sf_id, rep( 1+ii, np_xxconst[ii] ) )
sf_id2 <- c( sf_id2, 1+ii )
} else {
sf_id2 <- c( sf_id2, NA )
}
sf_id3 <- c( sf_id3, 1+ii )
}
}
if( is.logical( nvc_x[ 1 ] ) == FALSE ){
for( ii in 1:length( np_xx ) ){
ii_count <- ii_count + 1
if(np_xx[ii]>0){
id <- c( id, rep( ii_count, np_xx[ii] ))
sf_id <-c(sf_id, rep( 1+nxf+ii, np_xx[ii] ) )
sf_id2 <- c( sf_id2, 1+nxf+ii )
} else {
sf_id2 <- c( sf_id2, NA )
}
sf_id3 <- c( sf_id3, 1+nxf+ii )
}
}
XX2 <-Mo[sf_id,sf_id]*outer(sfmean,sfmean)
EX2 <-t(t(as.matrix(EX)[,sf_id])*sfmean)
M <-sfsf[, -(1:2)] - EX2 - t(EX2) + XX2
m <-m0 - sfmean*mo[sf_id]
yy <- sum( y ^ 2 )
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 ) )## changed
# par0_id <- c( par0_id , max(par0_id) + ( 1 : ng ) )
# null_dum<- c( null_dum, rep( 0, ng ) )
#}
#pp_id <-max( id )
if( is.logical( nvc_xconst[ 1 ] ) == FALSE ){#nvc_xconst[ 1 ] != FALSE
par0 <- c( par0, rep( 1, nnxf ) )
par0_id <- c( par0_id , max(par0_id) + ( 1:nnxf ) )
null_dum0<- rep( 0, nnxf ) # added
null_dum0[np_xxconst==0]<-1# added
null_dum<- c( null_dum, null_dum0 )
}
if( is.logical( nvc_x[ 1 ] ) == FALSE ){#nvc_x[ 1 ] != FALSE
par0 <- c( par0, rep( 1, nnsv ) )
par0_id <- c( par0_id , max(par0_id) + ( 1:nnsv ) )
null_dum0<- rep( 0, nnsv ) # added
null_dum0[np_xx==0]<-1 # added
null_dum<- c( null_dum, null_dum0 )
}
id_exist<-unique(id)[unique(id)>0]
par0[(par0_id %in% id_exist) ==FALSE]<-0
Par <-par0
obj <- sd( y )
LL <- NULL
res_old<- Inf
n_omit <- 0
iter <- 1
gmess <- 1
warn <- FALSE
while( (( obj > 0.0001 ) & ( iter <= maxiter ))|( iter <= 4 ) ){#sd( y ) *
if( !is.null(x) ) print( paste( "------- Iteration ", iter, " -------", sep = "" ) )
LL0 <- LL
n_omit <- 0
for( par0_sel in (1:( nsv + nnxf+ ng + nnsv ))[id_exist]){
evSqrt <- NULL
par0_sq <- par0 ^ 2
for( i in 1:nsv ){
evv <- ev ^ par0_sq[ nsv + i ] * sum( ev ) / sum( ev ^ par0_sq[ nsv + i ] )
evSqrt <- c( evSqrt, par0_sq[ i ] * sqrt( evv ) )
}
M0 <- M
for( j in 1:nsv ){
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 ){
for( i2 in 1:ng ){
xgg <- rep( 1, sum( id == ( nsv + i2 ) ) )
evSqrt <- c( evSqrt, par0_sq[ 2 * nsv + i2 ] * xgg )
}
for( j2 in 1:ng ){
if( j2 != ( par0_sel - nsv ) ){
id_sub <- ( id == ( j2 + nsv ) )
diag( M0[ id_sub, id_sub ] )<-
diag( M0[ id_sub, id_sub ] ) + 1/evSqrt[ id[ - ( 1:nx ) ] == (j2 + nsv) ] ^ 2
} else {############################### added
id_sub <- ( id == ( j2 + nsv ) )
g_add <- (1/evSqrt[ id[ - ( 1:nx ) ] == (j2 + nsv) ] ^ 2)/100
diag( M0[ id_sub, id_sub ] )<-
diag( M0[ id_sub, id_sub ] ) + g_add
} ############################### added
}
}
if( nnxf != 0 ){
for( i2 in 1:nnxf ){
xgg <- rep( 1, sum( id == ( nsv+ ng + i2 ) ) )
evSqrt <- c( evSqrt, par0_sq[ (2 * nsv + ng ) + i2 ] * xgg )
}
for( j2 in 1:nnxf ){
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 ){
for( i2 in 1:nnsv ){
xgg <- rep( 1, sum( id == ( nsv + ng + nnxf + i2 ) ) )
evSqrt <- c( evSqrt, par0_sq[ (2 * nsv + ng + nnxf ) + i2 ] * xgg )
}
for( j2 in 1:nnsv ){
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
}
}
}
null_dum2<- null_dum
null_dum2[ par0_sel ] <- 0
null_dum3<- c( rep( 0, nx ), null_dum2[ id ]) == 0
MM <- M [ null_dum3, null_dum3 ]
MM0 <- M0[ null_dum3, null_dum3 ]
mm <- m [ null_dum3 ]
idd <- id[ null_dum3 ]
id_omit1<- c( diff( id ), 1)
#if( ng != 0 ) id_omit1[ id > nsv ] <- 0
#if( ng != 0 ) id_omit1[( nsv < id ) & ( id <= nsv + ng )] <- 0
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
sstop <-FALSE
if( ( ng != 0 )&( nsv < par0_sel )&( par0_sel <= nsv + ng ) ){
if( rcond( MM0 ) < 1e-30 ){
loglik <- ( - 1 / 2 ) * res_old
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)
#if(gmess == 1 ){
message( paste( "Note:", "group effect", par0_sel - nsv, "is omitted to stablize the estimates", sep = " " ) )
#}
gmess <- gmess + 1
sstop <-TRUE
}
}
if( sstop==FALSE ){
er_dum <- TRUE
n_omit0 <- 0
while( er_dum == TRUE ){
try1 <- try( M0inv <- solve( 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 <- ( class(try1)[1] =="try-error" ) | ( class(try2)[1] =="try-error" )
er_dum <- inherits( try1, "try-error" )|inherits( try2, "try-error" )
if( er_dum == TRUE ){
M <- M [ id_omit1 == 0, id_omit1 == 0 ]
M0 <- M0[ id_omit1 == 0, id_omit1 == 0 ]
m <- m [ id_omit1 == 0 ]
id <- id[ id_omit1 == 0 ]
#X <- X [, id_omit1 == 0 ] ### difference
null_dum3<- null_dum3[ id_omit1 == 0 ]
MM <- M [ null_dum3, null_dum3 ]
MM0 <- M0[ null_dum3, null_dum3 ]
mm <- m [ null_dum3 ]
idd <- id[ null_dum3 ]
id_omit1<- c( diff( id ), 1)
#if( ng != 0 ) id_omit1[ id > nsv ] <- 0
#if( ng != 0 ) id_omit1[( nsv < id ) & ( id <= nsv + ng )] <- 0 ## deleted
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
n_omit0 <- n_omit0 + 1
nev <-nev-1 ## added
warn <- TRUE
}
if( n_omit0 >= 1000){
stop( " Singular fit. Simplify the model" )
}
}
if( n_omit0 > 0 ){
message( paste( "Note: ", n_omit0, " eigenvectors are omitted to stablize the estimates", sep = "" ) )
}
n_omit <- c( n_omit, n_omit0 )
ev <- ev [ 1:sum( id == 1 ) ]
#meig$sf <- meig$sf[ , 1:sum( id == 1 ) ] ##################3 changed
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 ) )#fixed
} else {
par00_id<- max( which( Par[,par0_id == par0_sel] != 0 ) )## fixed
}
par00 <- Par[ par00_id, par0_id == par0_sel ]
}
if( n_omit0 > 0 ){
res_old <- lik_resf_vc0( par0, ev = ev, M = MM, m = mm, yy = yy,nnxf=nnxf,nnsv=nnsv,
n = n, nx = nx, nsv = nsv, ng = ng, emet = method, null_dum4 = null_dum2 )
}
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, 4 )
omethod <- "L-BFGS-B"##### error length of b[-(1:nx)] and evSqrt are inconsistent
#M = MM; M0inv = M0inv; m = mm; emet = method; id = idd
#lower = llim; upper = ulim; method = omethod
res <- optim( fn = lik_resf_vc, par00, par0 = par0, ev = ev, M = MM, M0inv = M0inv,
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, ng = ng, emet = method, id = idd,
par0_sel = par0_sel, par0_id = par0_id, null_dum2 = null_dum2,
lower = llim, upper = ulim, method = omethod )
res_int <- res
} else if( par0_sel <= nsv ){
res <- optim( fn = lik_resf_vc, par00, par0 = par0, ev = ev, M = MM, M0inv = M0inv,
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, ng = ng, emet = method, id = idd,
par0_sel = par0_sel, par0_id = par0_id, null_dum2 = null_dum2 )
} else if( par0_sel <= nsv+ng ){#### changed
res <- optimize( f = lik_resf_vc, lower = parVmax_sq/1000, upper = parVmax_sq, par00,
par0 = par0, ev = ev, M = MM, M0inv = M0inv,
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, ng = ng, emet = method, id = idd,
par0_sel = par0_sel, par0_id = par0_id, null_dum2 = null_dum2 )
res$value <- c( res$objective )
res$par <- c( res$minimum )
} else {
res <- optimize( f = lik_resf_vc, lower = parVmax_sq/10^5, upper = 10^10, par00, par0 = par0, ev = ev,
M = MM, M0inv = M0inv,
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, ng = ng, emet = method, id = idd,
par0_sel = par0_sel, par0_id = par0_id, null_dum2 = null_dum2 )
res$value <- c( res$objective )
res$par <- c( res$minimum )
}
if( ( iter > 3 ) & ( res$value > res_old ) ){
loglik <- ( - 1 / 2 ) * res_old
} else {
if( par0_sel != 1 ){#### ) & ( par0_sel <= nsv ))|( par0_sel > nsv + ng )
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, ev = ev, M = MM_ref, m = mm_ref, yy = yy,nnxf=nnxf,nnsv=nnsv,
n = n, nx = nx, nsv = nsv, ng = ng, emet = method, null_dum4 = null_dum4 )
} else {
res_ref <- Inf
}
np_add <- length(par00)
flag <- 0
if( par0_sel <= nsv ){
if( allsvc_flag == 0 ){ ####### added
allvc <- allsvc
} else {
flag <-1
allvc <- c(TRUE, allsvc_id)[ par0_sel ] # including intercept
}
} else if( ( nsv < par0_sel )&( par0_sel <= nsv + ng )){
allvc <- TRUE#allg
} else if( ( nsv +ng < par0_sel )&( par0_sel <= nsv + ng + nnxf )){
if( allnvc_xconst_flag == 0 ){ ####### added
allvc <- allnvc_xconst
} else {
flag <-1
allvc <- allnvc_xconst_id[ par0_sel - nsv - ng ]
}
} else {
if( allnvc_x_flag == 0 ){ ####### added
allvc <- allnvc_x
} else {
flag <-1
allvc <- allnvc_x_id[ par0_sel - nsv - ng - nnxf ]
}
}
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)
} 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_old + pen * (nx+sum(par0!=0)+1)
}
}
}
LL0 <- c( LL0, loglik )
if( !is.null(x) ) print( paste( par0_sel, "/", nsv + ng + nnxf + nnsv, sep = "" ) )
}
if( iter > 1 ){
if( sum( n_omit ) == 0 ){
obj <- abs( loglik - loglik0 )
} else {
obj <- Inf
}
}
Par <- rbind( Par, par0 )
loglik0 <- loglik
LL <- c( LL, loglik )
if( !is.null(x) ) print( paste( toupper( penalty ), ": ", round( score, 3 ), sep = "" ) )
iter <- iter + 1
}
par2 <- par0 ^ 2
evSqrt <- NULL
evSqrt2 <- list(NULL)
for( i in 1:nsv ){
evv <- ev ^ par2[ nsv + i ] * sum( ev ) / sum( ev ^ par2[ nsv + i ] )
evSqrt<- c( evSqrt, par2[ i ] * sqrt( evv ) )
evSqrt2[[ i ]] <- par2[ i ] * sqrt( evv )
}
if( ng != 0 ){
for( i2 in 1:ng ){
xgg <- rep( 1, sum( id == ( nsv + i2 ) ) )
evSqrt <- c( evSqrt, par2[ 2 * nsv + i2 ] * xgg )
evSqrt2[[ nsv + i2 ]] <- par2[ 2 * nsv + i2 ] * xgg
}
}
if( nnxf != 0 ){
for( i2 in 1:nnxf ){
xgg <- rep( 1, sum( id == ( nsv + ng + i2 ) ) )
evSqrt <- c( evSqrt, par2[ 2 * nsv + ng + i2 ] * xgg )
evSqrt2[[ nsv + ng + i2 ]] <- par2[ 2 * nsv + ng + i2 ] * xgg
}
}
if( nnsv != 0 ){
for( i2 in 1:nnsv ){
xgg <- rep( 1, sum( id == ( nsv + ng + nnxf + i2 ) ) )
evSqrt <- c( evSqrt, par2[ 2 * nsv + ng + nnxf + i2 ] * xgg )
evSqrt2[[ nsv + ng + nnxf + i2 ]] <- par2[ 2 * nsv + ng + nnxf + i2 ] * xgg
}
}
MM <- M [ null_dum3, null_dum3 ]
MM0 <- M0[ null_dum3, null_dum3 ]
mm <- m [ null_dum3 ]
idd <- id[ null_dum3 ]
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 ] )
diag( M[ -( 1:nx ), -( 1:nx ) ] ) <- diag( M[ -( 1:nx ), -( 1:nx ) ] ) + 1
MM <- M[ null_dum3, null_dum3 ]
MMinv <- solve( MM, tol=tol )
mm <- m[ null_dum3 ]
eevSqrt <- evSqrt[ null_dum3[ -( 1:nx )] ]
mm[ -( 1:nx ) ] <- mm[ -( 1:nx ) ] * eevSqrt
b <- MMinv %*% mm
b[ -( 1:nx ) ] <- b[ -( 1:nx ) ] * eevSqrt
null_dum3_id<- which(null_dum3)
sf_id4 <- c(nx,rep(nev,nsv),np_xxconst,np_xx) #fixed id_omit1
sf_id5 <- 1:nsv
eff_id <- rep(1,nsv)
if( nnxf != 0) {
sf_id5 <- c( sf_id5, nsv + 1:nnxf )
eff_id <- c( eff_id, rep(2, nnxf))
}
if( nnsv != 0){
sf_id5 <- c( sf_id5, 2:nsv )#not nsv
eff_id <- c( eff_id, rep(3, nnsv))#here is nnsv
}
sec_id <- 1:( nsv+ ng + nnxf + nnsv )
#b_s <- list(NULL)
#b_covs <- list(NULL)
bvc_all <- foreach(i = 1:(length(ids)-1), .packages=c("fields","splines"),.combine="rbind") %dopar% {
if(i==1){
ids_sub<-ids[i]:ids[i+1]
} else {
ids_sub<-(ids[i]+1):ids[i+1]
}
n_sub <- length(min(ids_sub):max(ids_sub) )
b_vc_s0 <- matrix(0, nrow= length(ids_sub), ncol= nsv );j_s<-1
b_vc_n0 <- matrix(0, nrow= length(ids_sub), ncol= nsv );j_n<-1
bse_vc_s0 <- matrix(0, nrow= length(ids_sub), ncol= nsv )
bse_vc_n0 <- matrix(0, nrow= length(ids_sub), ncol= nsv )
b_vc0 <- matrix(0, nrow= length(ids_sub), ncol= nsv )
bse_vc0 <- matrix(0, nrow= length(ids_sub), ncol= nsv )
if(nnxf>0) {
c_vc0 <- matrix(0, nrow= length(ids_sub), ncol= nnxf );j_c<-1
cse_vc0 <- matrix(0, nrow= length(ids_sub), ncol= nnxf )
} else { ## added
c_vc0 <-NULL
cse_vc0 <-NULL
}
for( j in unique( sf_id5 ) ){#( nsv+ ng + nnxf + nnsv )
bid0_NA <- sf_id2[ sf_id5 == j ]#ifelse( j == 1, 1, nxf + j )
bid0 <- sf_id3[ sf_id5 == j ]#ifelse( j == 1, 1, nxf + j )
null_dum_sub<- null_dum[ sf_id5 == j ]
b_full_id <- cumsum( sf_id4[ c(0, null_dum) == 0] )
if( (j %in% sf_id5[null_dum==0])==FALSE ){### if vc is not selected
if( j <= nsv ){
b_vc0 [,j]<- b_vc_s0 [,j]<- b_vc_n0 [,j]<- b[bid0[1]] ##### added
bse_vc0[,j]<- bse_vc_s0[,j]<- bse_vc_n0[,j]<- c(sqrt(MMinv[bid0[1], bid0[1]]))#### added
#b_vc0[ ,j ] <- b[ bid0[ 1 ] ]
#bse_vc0[ ,j ]<- c( sqrt( MMinv[ bid0[ 1 ], bid0[ 1 ] ] ) )
j_s <- j_s + 1
j_n <- j_n + 1
} else {
c_vc0[ ,j_c ] <- b[ bid0[ 1 ] ]
cse_vc0[ ,j_c ]<- c( sqrt( MMinv[ bid0[ 1 ], bid0[ 1 ] ] ) )
j_c <- j_c + 1
}
} else {
sf_id4_start<- b_full_id[ which(sf_id5[null_dum==0] == j) ] + 1
sf_id4_end <- b_full_id[ which(sf_id5[null_dum==0] == j) + 1 ]
eff_id_sub <- eff_id[null_dum==0][ sf_id5[null_dum==0] == j ]#(1,1,1,3,3)#(1=snvc, 3=nvc)
sec_id_sub <- sec_id[null_dum==0][ sf_id5[null_dum==0] == j ]#(1,2,3,4,5)
bid <- list( NULL )
sfmean_sub <- list(NULL)
bid[[ 1 ]] <- bid0[1]
bid_len <- length( sf_id4_start )
for(jjj in 1:bid_len){
sell <- ifelse(eff_id_sub[jjj]==1, 2, 3) #(1=snvc, 3=nvc)
bid[[sell]] <- sf_id4_start[ jjj ]:sf_id4_end[ jjj ]#selct column
sfmean_sub[[sell]]<-sfmean[ bid[[sell]] ]
}
if( 1 %in% eff_id_sub ){
if( (null_dum_sub[ 1 ] == 0)&( !is.na( bid0_NA[ 1 ] ) ) ){
if( covmodel == "exp" ){
ccc0 <- exp( -rdist( coordk, coords[ids_sub,] ) / h )
} else if( covmodel == "gau" ){
C <- exp( -( D / h) ^ 2 )
ccc0 <- exp( -( rdist( coordk, coords[ids_sub,] ) / h )^2 )
} else if( covmodel == "sph" ){
ddd0 <- rdist( coordk, coords[ids_sub,] )
ccc0 <- ifelse( ddd0 < h , 1 - 1.5 * ( ddd0 / h ) + 0.5 * ( ddd0 / h ) ^ 3, 0 )
}
basis_s <- t( ccc0 - Cmean ) %*% sf_ap0[ , 1:nev ]
basis_s <- t( t( basis_s ) - sfmean_sub[[ 2 ]] )
b_vc_s0[ ,j_s ]<- basis_s %*% b[ bid[[ 2 ]] ]
MMinv_sub_s <- MMinv[ unlist(bid[ 1:2 ]), unlist(bid[ 1:2 ])] ## Added below #########
basis_sel_s <- evSqrt2[which(sf_id5==j)][[1]]#[!is.null(basis_s)]
basis2_s <- t( t( basis_s ) * unlist( basis_sel_s ))## added
x_basis_s <- as.matrix( cbind( 1, basis2_s ) )###added
bse_vc_s0[ ,j ]<- sqrt( colSums( t( x_basis_s ) * ( MMinv_sub_s %*% t( x_basis_s ))))### Added above ##
} else {
basis_s <- NULL
bse_vc_s0[ ,j ]<- c( sqrt( MMinv[ bid0[ 1 ], bid0[ 1 ] ] ) )
}
j_s <- j_s + 1
} else if(1 %in% eff_id[ sf_id5==j ]){
basis_s <- NULL
bse_vc_s0[ ,j ]<- c( sqrt( MMinv[ bid0[ 1 ], bid0[ 1 ] ] ) )
j_s <- j_s + 1
}
if( 2 %in% eff_id_sub ){
if( (null_dum_sub[ 1 ] == 0)&( !is.na( bid0_NA[ 1 ] ) ) ){
jj <- sec_id_sub[ 1 ] - nsv - ng
basis_c<- as.matrix( pmod_xconst[[ jj ]][ids_sub, ] )
basis_c<- t( t(basis_c) - sfmean_sub[[ 3 ]])
c_vc0[ ,j_c ] <- b[ bid0[ 1 ] ] + basis_c %*% b[ bid[[ 3 ]] ]
MMinv_sub <-MMinv[ unlist(bid), unlist(bid)]
basis2 <- t( t( basis_c ) * unlist( evSqrt2[which(sf_id5==j)] ) )
x_basis <- as.matrix( cbind( 1, basis2 ) )###added
cse_vc0[ ,j_c ]<- sqrt( colSums( t( x_basis ) * ( MMinv_sub %*% t( x_basis ))))
} else {
basis_c <- NULL
c_vc0[ ,j_c ] <- c(b[ bid0[ 1 ] ])#rep( b[ bid0[ 1 ] ],n_sub )
cse_vc0[ ,j_c ]<- c( sqrt( MMinv[ bid0[ 1 ], bid0[ 1 ] ] ) )
#rep(sqrt( MMinv[ bid0[ 1 ], bid0[ 1 ] ] ),n_sub)
}
j_c <- j_c + 1
}
if( 3 %in% eff_id_sub ){
if( (null_dum_sub[ 2 ] == 0)&( !is.na( bid0_NA[ 2 ] ) ) ){
jj <- sec_id_sub[ length(sec_id_sub) ] - nsv - ng - nnxf
basis_n<- as.matrix( pmod_x[[ jj ]][ids_sub, ] )
#basis_n<- as.matrix( predict(pmod_x[[ jj ]],X1_nvc[ids_sub ,jj ]) )#[ ,-1 ]
basis_n<- t( t(basis_n) - sfmean_sub[[ 3 ]])
b_vc_n0[ ,j_n ]<- basis_n %*% b[ bid[[ 3 ]] ]
MMinv_sub_n <- MMinv[ unlist(bid[ c(1,3) ]), unlist(bid[ c(1,3) ])] ## Added below #########
basis_sel_n <- evSqrt2[which(sf_id5==j)][[2]]#[!is.null(basis_n)]
basis2_n <- t( t( basis_n ) * unlist( basis_sel_n ))## added
x_basis_n <- as.matrix( cbind( 1, basis2_n ) )###added
bse_vc_n0[ ,j ]<- sqrt( colSums( t( x_basis_n ) * ( MMinv_sub_n %*% t( x_basis_n ))))### Added above ##
} else {
basis_n <- NULL
bse_vc_n0[ ,j ]<- c( sqrt( MMinv[ bid0[ 1 ], bid0[ 1 ] ] ) )
}
j_n <- j_n + 1
} else if(1 %in% eff_id[sf_id5==j] ){#eff_id_sub
basis_n <- NULL
bse_vc_n0[ ,j ]<- c( sqrt( MMinv[ bid0[ 1 ], bid0[ 1 ] ] ) )
j_n <- j_n + 1
}
if( 1 %in% eff_id[sf_id5==j] ){#eff_id_sub
b_vc0[ ,j ] <- b[ bid0[ 1 ] ] + b_vc_s0[ ,j ] + b_vc_n0[ ,j ]
b_vc_s0[ ,j ]<- b[ bid0[ 1 ] ] + b_vc_s0[ ,j ] ### added
b_vc_n0[ ,j ]<- b[ bid0[ 1 ] ] + b_vc_n0[ ,j ] ### added
MMinv_sub <- MMinv[ unlist(bid), unlist(bid) ]
if( !is.null(basis_s)|!is.null(basis_n) ){
basis2 <- cbind( basis_s, basis_n )
basis_sel <-evSqrt2[which(sf_id5==j)][c(!is.null(basis_s), !is.null(basis_n))]# fixed
basis2 <- t( t( basis2 ) * unlist( basis_sel ) )
x_basis <- as.matrix( cbind( 1, basis2 ) )
bse_vc0[ ,j ]<- sqrt( colSums( t( x_basis ) * ( MMinv_sub %*% t( x_basis ))))
} else {
bse_vc0[ ,j ]<- c( sqrt( MMinv[ bid0[ 1 ], bid0[ 1 ] ] ) )#rep(sqrt( MMinv[ bid0[ 1 ], bid0[ 1 ] ] ),n_sub)
}
}
}
}
cbind(b_vc0, bse_vc0, b_vc_s0, b_vc_n0, bse_vc_s0, bse_vc_n0, c_vc0, cse_vc0)
}
stopCluster(ncores)
bvc_all <-as.matrix(bvc_all)
b_vc <-as.matrix(bvc_all[,1:nsv])
bse_vc <-as.matrix(bvc_all[,( nsv+1):(2*nsv)])
b_vc_s0 <-as.matrix(bvc_all[,(2*nsv+1):(3*nsv)])
b_vc_n0 <-as.matrix(bvc_all[,(3*nsv+1):(4*nsv)])
bse_vc_s0 <-as.matrix(bvc_all[,(4*nsv+1):(5*nsv)])
bse_vc_n0 <-as.matrix(bvc_all[,(5*nsv+1):(6*nsv)])
if( nnxf !=0 ){
c_vc <-as.matrix(bvc_all[,(6*nsv+1):(6*nsv + nnxf)])
cse_vc <-as.matrix(bvc_all[,(6*nsv + nnxf + 1):(6*nsv + 2*nnxf)])
} else {
c_vc <- NULL
cse_vc <- NULL
}
pred <- b_vc[,1]
if( nsv > 1 ) pred <- pred + rowSums( X1 * b_vc[, -1 ] )
if( nnxf !=0 ){
pred <- pred + rowSums( Xconst_nvc * c_vc )
} else {
c_vc <- NULL
if( nxf > 0){
pred <- pred + as.matrix( Xconst ) %*% c(b[ 2:( nxf + 1 ) ])
}
}
resid <- y - pred
SSE <- sum( resid ^ 2 )
SSY <- sum( ( y - mean( y ) ) ^ 2 )
sig <- SSE / ( n - nx )
bse_vc <-bse_vc * sqrt( sig ) #### added
moran_vc <- rep( 0, nsv )
j0 <- 1
for(j in 1:nsv){
bid0 <- ifelse( j == 1, 1, nxf + j )
if( null_dum[ j ] == 0 ){
bid0_vc <- ( nx + nev * ( j0 - 1 ) + 1 ):( nx + nev * j0 )
moran_vc[ j ] <- sum(b[ bid0_vc ]^2*ev)/(ev[1]*sum(b[ bid0_vc ]^2))
j0 <- j0 + 1
} else {
moran_vc[ j ] <- NA
}
}
parV <- par2[ 1:nsv ]
parR <- par2[ (nsv + 1):( 2 * nsv ) ]
if( nnsv !=0 ){
parNsv <- par2[ (2 * nsv + ng + nnxf + 1 ):( 2 * nsv + ng + nnxf + nnsv ) ]
} else {
parNsv <- NULL
}
if( nnxf !=0 ){
parNxf <- par2[ (2 * nsv +ng + 1 ):( 2 * nsv + ng + nnxf) ]
} else {
parNxf <- NULL
}
nsv2 <- sum( parV != 0 )
nnxf2 <- sum( parNxf != 0 )
nnsv2 <- sum( parNsv != 0 )
np <- nx + nsv2 * 2 + ng + nnxf2 + nnsv2 + 1#
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 )
b_cov <- sig * MMinv
bse <- sqrt( diag( b_cov ) )
if( is.null( x ) ){
#b_vc[,1]<- b_vc[,1]# + Bias
mean_bvc<- mean( b_vc[,1] )
b_vc[,1]<- b_vc[,1] - mean_bvc
b_c <- b[ 1:( nxf + 1 ) ]
b_c[1] <- mean_bvc
r <- b[( nx + 1 ):( nx + nev) ]
bse_c <- bse[ 1:( nxf + 1 ) ]
bz_c <- b_c / bse_c
bp_c <- pnorm(abs(bz_c),lower.tail=FALSE) * 2
b_par <- data.frame( Estimate = b_c, SE = bse_c, t_value = bz_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 ) ]
bz_c <- as.matrix( b_c / bse_c )
bp_c <- pnorm(abs(bz_c),lower.tail=FALSE) * 2
b_par <- data.frame( Estimate = b_c, SE = bse_c, z_value = bz_c, p_value = bp_c )
rownames( b_par )<- xfname
} else {
b_par <- NULL
}
r <- NULL
}
bz_vc <- b_vc / bse_vc
bp_vc <- pnorm(abs(bz_vc),lower.tail=FALSE) * 2
b_vc <- data.frame( b_vc )
bse_vc <- data.frame( bse_vc )
bz_vc <- data.frame( bz_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( bz_vc ) <- c( "(Intercept)", xname )
names( bp_vc ) <- c( "(Intercept)", xname )
names( bp_vc_adj )<- c( "(Intercept)", xname )
if(is.null(b_vc_s0) ==FALSE ){
bz_vc_s0<- as.matrix( b_vc_s0 / bse_vc_s0 )
bp_vc_s0<- pnorm(abs( bz_vc_s0 ),lower.tail=FALSE) * 2
b_vc_s0 <- data.frame( b_vc_s0 )
bse_vc_s0 <- data.frame( bse_vc_s0 )
bz_vc_s0 <- data.frame( bz_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( bz_vc_s0 ) <- c( "(Intercept)", xname )
names( bp_vc_s0 ) <- c( "(Intercept)", xname )
B_vc_s <-list(b_vc = b_vc_s0, bse_vc = bse_vc_s0, bz_vc = bz_vc_s0, bp_vc = bp_vc_s0)
}
if(is.null(b_vc_n0) ==FALSE ){
bz_vc_n0<- as.matrix( b_vc_n0 / bse_vc_n0 )
bp_vc_n0<- pnorm(abs( bz_vc_n0 ),lower.tail=FALSE) * 2
b_vc_n0 <- data.frame( b_vc_n0 )
bse_vc_n0 <- data.frame( bse_vc_n0 )
bz_vc_n0 <- data.frame( bz_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( bz_vc_n0 ) <- c( "(Intercept)", xname )
names( bp_vc_n0 ) <- c( "(Intercept)", xname )
B_vc_n <-list(b_vc = b_vc_n0, bse_vc = bse_vc_n0, bz_vc = bz_vc_n0, bp_vc = bp_vc_n0)
}
if( nnxf != 0 ) {
cse_vc <- cse_vc * sqrt( sig ) #### added
cz_vc <- as.matrix( c_vc / cse_vc )
cp_vc <- pnorm(abs(cz_vc),lower.tail=FALSE)*2
c_vc <- data.frame( c_vc )
cse_vc <- data.frame( cse_vc )
cz_vc <- data.frame( cz_vc )
cp_vc <- data.frame( cp_vc )
names( c_vc ) <- xfname
names( cse_vc ) <- xfname
names( cz_vc ) <- xfname
names( cp_vc ) <- xfname
} else {
c_vc <- NULL
cse_vc <- NULL
cz_vc <- NULL
cp_vc <- NULL
}
parV <- parV * sqrt( sig )
sf_par <- data.frame( rbind( parV, moran_vc ) )
names( sf_par ) <- c( "(Intercept)", xname )
rownames( sf_par )<- c( "random_SD", "Moran.I/max(Moran.I)" )
e_stat <- data.frame( stat = c( sqrt( sig ), r2, loglik, AIC, BIC ) )
rownames( e_stat ) <- c( "resid_SE", "adjR2(cond)", lik_nam, "AIC", "BIC")
vc <- data.frame(ifelse( sf_par[1,] ==0, 0, 1) )
names( vc ) <- names( sf_par )
rownames( vc ) <- "Spatial"
if( is.null( parNsv ) ==FALSE ){
vc <- rbind( vc, c(0, ifelse( parNsv ==0, 0, 1 )))
rownames( vc )[2] <- "Non-spatial"
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)")
names(sn_rat)[1]<-"(Intercept)"
s <- list(sf_par, s_n, sn_rat)
} 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_SD"
names( bnf_par ) <- xxfname
s_nxconst<- bnf_par
}
if( sum( n_omit ) >5 ) {
if( e_stat[2,1]> -0.2 ){
message( "Note: The model is nearly singular. Consider simplifying the model" )
} else {
message( "Note: Singular fit. Simplify the model")
}
} else if( e_stat[2,1]< -0.2 ){
message("Note: Singular fit. Simplify the model")
}
#evSqrt2 # evSqrts = evSqrts, evSqrts_n = evSqrts_n,### b_s = bb, b_covs = bb_cov,
other <- list( res_int =res_int, r = r, sf_alpha = parR, x_id = x_id, nxf = nxf,np=np, xf_id = xf_id,# df = df,
model = "resf_vc", nvc_x=nvc_x, nvc_xconst=nvc_xconst, nvc_num = nvc_num, dif=1,
method=method, penalty=penalty, y_nonneg =FALSE, w_scale = 1, ncores = ncores,
y = y, x = x, xconst = xconst, coords = coords, weight = NULL )#, Bias=bias
result <- list( b_vc = b_vc, bse_vc = bse_vc, z_vc = bz_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 = c_vc, cse_vc = cse_vc,
cz_vc = cz_vc, cp_vc = cp_vc, s = s, s_c=s_nxconst, vc = vc, e = e_stat,
pred = pred, resid = resid, other = other, call = match.call() )
class(result)<-"besf_vc"
return( result )
}
print.besf_vc <- function(x, ...)
{
cat("Call:\n")
print(x$call)
if( is( x$s )[1] == "list" ){
cat("\n----Spatially and non-spatially varying coefficients on x (summary)----\n")
} else {
cat("\n----Spatially varying coefficients 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-spatially 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")
cat("\nSpatial effects (coefficients on x):\n")
if( is( x$s )[1] != "list" ){
print( x$s )
} else {
print(x$s[[1]])
cat("\nNon-spatial effects (coefficients on x):\n")
print(x$s[[2]])
}
if( !is.null(x$s_c) ){
cat("\nNon-spatial effects (coefficients on xconst):\n")
print(x$s_c)
}
#if( !is.null(x$s_g) ){
# cat("\nGroup variation:\n")
# print(x$s_g)
#}
cat("\n----Error statistics-------------------------------------\n")
print(x$e)
if( x$other$method=="reml"){
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')
}
invisible(x)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.