R/back.search.R

In MRPP: Multiresponse permutation procedure and its variable importance and variable selection methods.

mrppBVS <-
function(y,permutedTrt, 
         importance=c('dp.dw','p.dd.dw'),
         alpha.in, #=if(match.arg(importance)=='dp.dw') 0 else 0.1, 
         alpha.del=0, size.in=0L, stepwise=FALSE, niter=Inf, verbose=TRUE, ...)
## y is a data matrix, with col's being variables and rows being observations
{
    if(!is.matrix(y))y=as.matrix(y)
    importance=match.arg(importance)
    if(missing(alpha.in)) alpha.in=if(importance=='dp.dw') 0 else 0.1
    N=nrow(y)
#    if(missing(cperm.mat)) cperm.mat=apply(permutedTrt,2,function(kk)(1:N)[-kk])
    B=nperms.permutedTrt(permutedTrt)
    importance=match.arg(importance)

    ans=vector('list'); attr(ans, 'parameter')=list(importance=importance, alpha.in=alpha.in, alpha.del=alpha.del, size.in=size.in, stepwise=stepwise, niter=niter)
    R=ncol(y)
    idx=1:R  # inclusion set
    xcl=integer(0L)
    i=1L
    imptnc.threshold=Inf
    imptnc=rep(-Inf, R)
    next.deleted.p=NA_real_
	ret.msg=c('All remaining variables meet importance requirements',
			  'Deleted variables fail to meet unimportance requirements if deletion continues', 
			  'Minimum number of selected variables reached',
			  'Maximum number of iterations reached')
    repeat{
        if(verbose && isTRUE(i%%verbose==0L)) {cat('iteration',i-1L,'...')
                    time0=proc.time()[3L]}
        idx=idx[imptnc<imptnc.threshold]
        if(length(idx)==0){
            ans[[i]]=list(iter=i-1L,var.idx=integer(0L), influence=numeric(0L), 
                        p.value=numeric(0L),deleted.p.value=ans[[1L]]$p.value)
            if(verbose)message('no variables left in the inclusion set')
			 attr(ans, 'status')=c(attr(ans, 'status'), ret.msg[3])
            return(ans)
        }
        dist0=dist(y[,idx,drop=FALSE])
        mrpp.rslt=mrpp.test.dist(dist0,permutedTrt=permutedTrt,...)
        mrpp.stats0=mrpp.rslt$all.statistics
        imptnc=if(importance=='dp.dw') 
                    get.dp.dw.kde(y[,idx,drop=FALSE],permutedTrt,distObj=dist0,mrpp.stats=mrpp.stats0, scale=B, ...) 
               else get.p.dd.dw(y[,idx,drop=FALSE],permutedTrt,...)

        var.rank=order(imptnc)
        ans[[i]]=list(iter=i-1L, var.idx=idx[var.rank], influence=imptnc[var.rank],
                      p.value=mrpp.rslt$p.value,
                      deleted.p.value=next.deleted.p)
        imptnc.threshold=if(stepwise) max(imptnc) else alpha.in
#        idx=idx[imptnc<max(imptnc)] else idx=idx[imptnc<alpha.in]
        if(alpha.del>=0) {
            xcl=c(xcl, idx[imptnc>=imptnc.threshold])
            dist.del=dist(y[,xcl,drop=FALSE])
            if(all(!is.na(dist.del)) && length(xcl)>0L)
#                ans[[i]]$deleted.p.value=mrpp.test.dist(dist.del, permutedTrt=permutedTrt)$p.value
                next.deleted.p=mrpp.test.dist(dist.del, permutedTrt=permutedTrt)$p.value
        }
        if(verbose && isTRUE(i%%verbose==0L)) {
          cat('\b\b\b:\t',length(idx),'variables left; mrpp.p =',ans[[i]]$p.value,';', 
                        'deleted.mrpp.p =',ans[[i]]$deleted.p.value,
                        ';', proc.time()[3L]-time0,'seconds passed;\n')
        }
        if( all(imptnc<=alpha.in) ) attr(ans, 'status') = c( attr(ans, 'status') , ret.msg[1L] )
		 if( i-1L>=niter ) attr(ans, 'status') = c( attr(ans, 'status') , ret.msg[4L] )
        if( isTRUE(next.deleted.p<=alpha.del) ) attr(ans, 'status') = c( attr(ans, 'status') , ret.msg[2L] )
        if( isTRUE(R-length(xcl)<size.in) ) attr(ans, 'status') = c( attr(ans, 'status') , ret.msg[3L] )
		 if (!is.null(attr(ans, 'status')))  return(ans)
        i=i+1L
#        if(stepwise) idx=idx[imptnc<max(imptnc)] else idx=idx[imptnc<alpha.in]
#        if(length(idx)==0) {
#            warning('not converged'); 
#            ans[[i]]=list(iter=i-1, var.idx=numeric(0), influence=numeric(0),
#                      p.value=numeric(0),
#                      deleted.p.value=ans[[1]]$p.value)
#            return(ans)}
    }
}