R/leaps.R

Defines functions leaps.from.biglm leaps.setup warn.extra leaps.exhaustive leaps.backward leaps.forward leaps.seqrep print.regsubsets print.summary.regsubsets summary.regsubsets regsubsets regsubsets.default regsubsets.biglm regsubsets.formula leaps vcov.regsubsets coef.regsubsets

Documented in coef.regsubsets leaps leaps.backward leaps.exhaustive leaps.forward leaps.seqrep leaps.setup print.regsubsets print.summary.regsubsets regsubsets regsubsets.biglm regsubsets.default regsubsets.formula summary.regsubsets vcov.regsubsets

# 
# R code for model selection using
# Alan Miller's FORTRAN routines
#


leaps.from.biglm<-function(x, force.in=NULL,
                      nvmax=8,nbest=1,warn.dep=TRUE){
    make.names<-function(np){
	if (np<27) letters[1:np] else as.character(1:np)
  }
  if(is.null(force.in) && attr(x$terms,"intercept")!=0) force.in<-1
  nn<-x$n
  np<-length(x$qr$D)
  index<-rep(0,np)
  names(index)<-x$names
  index[force.in]<--1
  force.in<-(index==-1) ## make force.in, force.out logical vectors
  ii<-order(index)
  force.in<-force.in[ii]
  ones<-rep(1,np)
  names(ones)<-colnames(x)
  first<-1+sum(ones[force.in])
  last<-np
  nvmax<-min(nvmax,np)

  vorder<-1:np
  il<-nvmax*(nvmax+1)/2
  nrbar<-np*(np-1)/2
  qrleaps<-x$qr
  tolset<-.Fortran("tolset",as.integer(np),
                   as.integer(nrbar),qrleaps$D,qrleaps$rbar,
                   tol=numeric(np),numeric(np),ier=as.integer(0), PACKAGE="leaps")
  if (tolset$ier!=0)
      warning(paste("TOLSET returned error code",tolset$ier))
  ss<-.Fortran("ssleaps",as.integer(np),qrleaps$D,
               qrleaps$thetab,qrleaps$ss,rss=numeric(np),
               ier=as.integer(0), PACKAGE="leaps")
  if (ss$ier!=0)
      warning(paste("SS returned error code",ss$ier))
 sing<-.Fortran("sing",np=as.integer(np),nrbar=as.integer(nrbar),
                d=qrleaps$D,rbar=qrleaps$rbar,thetab=qrleaps$thetab,
                sserr=qrleaps$ss,tol=tolset$tol,lindep=logical(np),
                work=numeric(np),ier=as.integer(0), PACKAGE="leaps")
   if (sing$ier>0)
       warning(paste("SING returned error code",sing$ier))
  sing$work<-NULL
  if(any(sing$lindep)) {
      if (warn.dep)
          warning(paste(sum(sing$lindep)," linear dependencies found"))
      if (any(sing$lindep[-1] & force.in)) stop("Linear dependency among variables forced in")
      rightorder<-sing$lindep 
      if (any((c(rightorder,1)-c(0,rightorder))<0)) {
          stop("Linear dependences in biglm fit: this can't happen")
      }
      lastsafe<-max((1:np)[!rightorder])
      if (lastsafe<min(nvmax,last)) {
          nvmax<-lastsafe
          if (warn.dep)
              warning(paste("nvmax reduced to ",nvmax))
      }
  }
  if (any(sing$lindep)){
      ss<-.Fortran("ssleaps",as.integer(np),sing$d,sing$thetab,
                   sing$sserr,rss=numeric(np),ier=as.integer(0),PACKAGE="leaps")
  	if (ss$ier!=0)
            warning(paste("SS returned error code",ss$ier))
  }
  initr<-.Fortran("initr",as.integer(np),as.integer(nvmax),as.integer(nbest),
                  bound=numeric(np),ress=numeric(nbest*nvmax),as.integer(nvmax),
                  lopt=integer(nbest*il),as.integer(il),vorder=as.integer(vorder),
                  ss$rss,ier=as.integer(0), PACKAGE="leaps")
  if (initr$ier!=0)
      warning(paste("INITR returned error code",initr$ier))
  nullrss<- ss$rss[1]
  rval<-c(sing,list(nn=nn,rss=ss$rss,bound=initr$bound,
                    ress=matrix(initr$ress,ncol=nbest),
                    lopt=matrix(initr$lopt,ncol=nbest),
                    nvmax=nvmax,nbest=nbest,nrbar=nrbar,il=il,
                    ir=nvmax,vorder=initr$vorder,
                    first=first,last=last,xnames=x$names,
                    force.in=(index==-1),force.out=(index==1),
                    intercept=TRUE,nullrss=nullrss))
  class(rval)<-"regsubsets"
  invisible(rval)
}


leaps.setup<-function(x,y,wt=rep(1,length(y)),force.in=NULL,
                      force.out=NULL,intercept=TRUE,
                      nvmax=8,nbest=1,warn.dep=TRUE){
  make.names<-function(np){
	if (np<27) letters[1:np] else as.character(1:np)
  }
  np<-NCOL(x)
  nn<-NROW(x)
  if (length(y)!=nn) stop("y and x different lengths")
  if (length(wt)!=nn) stop("wt and x different lengths")
  if (is.null(colnames(x))) colnames(x)<-make.names(np)
  index<-rep(0,np)
  names(index)<-colnames(x)
  index[force.in]<--1
  if (any(index[force.out]==-1)) stop("Can't force the same variable in and out")
  index[force.out]<-1
  force.in<-(index==-1) ## make force.in, force.out logical vectors
  force.out<-(index==1)
  ii<-order(index)
  xx<-x[,ii]
  force.in<-force.in[ii]
  force.out<-force.out[ii]
  ones<-rep(1,np)
  names(ones)<-colnames(x)
  first<-1+sum(ones[force.in])
  last<-np-sum(ones[force.out])
  nvmax<-min(nvmax,np)
  if (intercept){
    np<-np+1
    xnames<-c("(Intercept)",colnames(xx))    
    xx<-cbind(1,xx)
    colnames(xx)<-xnames
    first<-first+1
    last<-last+1
    nvmax<-nvmax+1
    index<-c(-1,index)
  }
  vorder<-1:np
  il<-nvmax*(nvmax+1)/2
  nrbar<-np*(np-1)/2
  qrleaps<-.Fortran("makeqr",np=as.integer(np),nn=as.integer(nn),
	wt=as.double(wt),tx=t(xx),y=as.double(y),d=numeric(np),
	rbar=numeric(nrbar),
        thetab=numeric(np),sserr=numeric(1),ier=as.integer(0),
        PACKAGE="leaps")
  if (qrleaps$ier!=0)
      warning(paste("MAKEQR returned error code",qrleaps$ier))
  qrleaps$tx<-NULL
  qrleaps$wt<-NULL
  tolset<-.Fortran("tolset",as.integer(np),
                   as.integer(nrbar),qrleaps$d,qrleaps$rbar,
                   tol=numeric(np),numeric(np),ier=as.integer(0), PACKAGE="leaps")
  if (tolset$ier!=0)
      warning(paste("TOLSET returned error code",tolset$ier))
  ss<-.Fortran("ssleaps",as.integer(np),qrleaps$d,
               qrleaps$thetab,qrleaps$sserr,rss=numeric(np),
               ier=as.integer(0), PACKAGE="leaps")
  if (ss$ier!=0)
      warning(paste("SS returned error code",ss$ier))
 sing<-.Fortran("sing",np=as.integer(qrleaps$np),nrbar=as.integer(nrbar),
                d=qrleaps$d,rbar=qrleaps$rbar,thetab=qrleaps$thetab,
                sserr=qrleaps$sserr,tol=tolset$tol,lindep=logical(qrleaps$np),
                work=numeric(qrleaps$np),ier=as.integer(0), PACKAGE="leaps")
   if (sing$ier>0)
       warning(paste("SING returned error code",sing$ier))
  sing$work<-NULL
  if(any(sing$lindep)) {
      if (warn.dep)
          warning(paste(sum(sing$lindep)," linear dependencies found"))
      if (any(sing$lindep[-1] & force.in)) stop("Linear dependency among variables forced in")
      rightorder<-sing$lindep | c(FALSE,force.out)
      if (any((c(rightorder,1)-c(0,rightorder))<0)) {
          if (warn.dep){
              cat("Reordering variables and trying again:\n")
          }
          reorder<-if(intercept) order(rightorder[-1]) else order(rightorder)
          rval<-leaps.setup(x[,ii[reorder],drop=FALSE],y,wt,
                            force.in[reorder],force.out[reorder],
                            intercept,nvmax,nbest,warn.dep=FALSE)
          rval$reorder<-ii[reorder]
          return(rval)
      }
      lastsafe<-max((1:np)[!rightorder])
      if (lastsafe<min(nvmax,last)) {
          nvmax<-lastsafe
          if (warn.dep)
              warning(paste("nvmax reduced to ",nvmax-intercept))
      }
  }
  if (any(sing$lindep)){
      ss<-.Fortran("ssleaps",as.integer(np),sing$d,sing$thetab,
                   sing$sserr,rss=numeric(np),ier=as.integer(0),PACKAGE="leaps")
  	if (ss$ier!=0)
            warning(paste("SS returned error code",ss$ier))
  }
  initr<-.Fortran("initr",as.integer(np),as.integer(nvmax),as.integer(nbest),
                  bound=numeric(np),ress=numeric(nbest*nvmax),as.integer(nvmax),
                  lopt=integer(nbest*il),as.integer(il),vorder=as.integer(vorder),
                  ss$rss,ier=as.integer(0), PACKAGE="leaps")
  if (initr$ier!=0)
      warning(paste("INITR returned error code",initr$ier))
  nullrss<-if (intercept) ss$rss[1] else sum(y^2)	
  rval<-c(sing,list(nn=qrleaps$nn,rss=ss$rss,bound=initr$bound,
                    ress=matrix(initr$ress,ncol=nbest),
                    lopt=matrix(initr$lopt,ncol=nbest),
                    nvmax=nvmax,nbest=nbest,nrbar=nrbar,il=il,
                    ir=nvmax,vorder=initr$vorder,
                    first=first,last=last,xnames=colnames(xx),
                    force.in=(index==-1),force.out=(index==1),
                    intercept=intercept,nullrss=nullrss))
  class(rval)<-"regsubsets"
  invisible(rval)
}

warn.extra<-function(obj){ 
	i<-which(obj$ress<obj$rss[seq_along(obj$ress)]*(1-1e-8),arr.ind=TRUE)[,1]
    if(length(i)) 
    	warning(paste0("model with initial (",paste(i-obj$int,collapse=","),") variables was better, and is reported"))
    invisible(i)	
}



leaps.exhaustive<-function(leaps.obj,really.big=FALSE){
    if (!inherits(leaps.obj,"regsubsets")){
        stop("Not a regsubsets object -- must run leaps.setup")
    }
    nbest<-leaps.obj$nbest
    if (!really.big & (leaps.obj$np>50 || leaps.obj$nbest>40)) {
        stop("Exhaustive search will be S L O W, must specify really.big=T")
    }
    dimwk<-3*leaps.obj$last
    dimiwk<-leaps.obj$nvmax
    rval<-.Fortran("xhaust",
                   np=as.integer(leaps.obj$np),
                   nrbar=as.integer(leaps.obj$nrbar),
                   d=leaps.obj$d,
                   rbar=leaps.obj$rbar,
                   thetab=leaps.obj$thetab,
                   first=as.integer(leaps.obj$first),
                   last=as.integer(leaps.obj$last),
                   vorder=as.integer(leaps.obj$vorder),
                   tol=leaps.obj$tol,
                   rss=leaps.obj$rss,
                   bound=leaps.obj$bound,
                   nvmax=as.integer(leaps.obj$nvmax),
                   ress=leaps.obj$ress,
                   ir=as.integer(leaps.obj$ir),
                   nbest=as.integer(leaps.obj$nbest),
                   lopt=matrix(as.integer(leaps.obj$lopt),ncol=nbest),
                   il=as.integer(leaps.obj$il),
                   wk=numeric(dimwk),
                   dimwk=as.integer(dimwk),
                   iwk=integer(dimiwk),
                   dimiwk=as.integer(dimiwk),
                   ier=as.integer(0), PACKAGE="leaps")
    rval$dimwk<-rval$dimiwk<-rval$iwk<-rval$wk<-NULL
    rval$xnames<-leaps.obj$xnames
    rval$method<-c("exhaustive",leaps.obj$method)
    rval$force.in<-leaps.obj$force.in
    rval$force.out<-leaps.obj$force.out
    rval$sserr<-leaps.obj$sserr
    rval$intercept<-leaps.obj$intercept
    rval$lindep<-leaps.obj$lindep
    rval$reorder<-leaps.obj$reorder
    rval$nullrss<-leaps.obj$nullrss
    rval$nn<-leaps.obj$nn
    class(rval)<-"regsubsets"
    if(rval$ier!=0) warning(paste("XHAUST returned error code",rval$ier))
    rval
}


leaps.backward<-function(leaps.obj,nested){
  if (!inherits(leaps.obj,"regsubsets")){
      stop("Not a regsubsets object -- must run leaps.setup")
  }
  nbest<-leaps.obj$nbest
  if(nested) leaps.obj$nbest<-0
  dimwk<-2*leaps.obj$last
  rval<-.Fortran("bakwrd",np=as.integer(leaps.obj$np),
                 nrbar=as.integer(leaps.obj$nrbar),d=leaps.obj$d,
                 rbar=leaps.obj$rbar,thetab=leaps.obj$thetab,
                 first=as.integer(leaps.obj$first),
                 last=as.integer(leaps.obj$last),
                 vorder=as.integer(leaps.obj$vorder),tol=leaps.obj$tol,
                 rss=leaps.obj$rss,bound=leaps.obj$bound,
                 nvmax=as.integer(leaps.obj$nvmax),
                 ress=leaps.obj$ress,ir=as.integer(leaps.obj$ir),
                 nbest=as.integer(leaps.obj$nbest),
                 lopt=matrix(as.integer(leaps.obj$lopt),ncol=nbest),
                 il=as.integer(leaps.obj$il),wk=numeric(dimwk),
                 dimwk=as.integer(dimwk),ier=as.integer(0), PACKAGE="leaps")
  rval$dimwk<-rval$wk<-NULL
  rval$xnames<-leaps.obj$xnames
  rval$method<-c("backward",leaps.obj$method)
  rval$force.in<-leaps.obj$force.in
  rval$force.out<-leaps.obj$force.out
  rval$sserr<-leaps.obj$sserr
  rval$intercept<-leaps.obj$intercept
  rval$lindep<-leaps.obj$lindep
  rval$reorder<-leaps.obj$reorder
  rval$nullrss<-leaps.obj$nullrss
  rval$nn<-leaps.obj$nn
  class(rval)<-"regsubsets"
  if(nested) {
      rval$nbest<-1
      rval$lopt[]<-rval$vorder[rval$lopt]
      rval$ress<-matrix(rval$rss,ncol=1)
  }
  if(rval$ier!=0)
      warning(paste("BAKWRD returned error code",rval$ier))
  warn.extra(rval)
  rval
}


leaps.forward<-function(leaps.obj,nested){
  if (!inherits(leaps.obj,"regsubsets")){
      stop("Not a regsubsets object -- must run leaps.setup")
  }
  nbest<-leaps.obj$nbest
  if(nested) leaps.obj$nbest<-0
  dimwk<-3*leaps.obj$last
  rval<-.Fortran("forwrd",np=as.integer(leaps.obj$np),
                 nrbar=as.integer(leaps.obj$nrbar),
                 d=leaps.obj$d,rbar=leaps.obj$rbar,
                 thetab=leaps.obj$thetab,first=as.integer(leaps.obj$first),
                 last=as.integer(leaps.obj$last),
                 vorder=as.integer(leaps.obj$vorder),
                 tol=leaps.obj$tol,rss=leaps.obj$rss,
                 bound=leaps.obj$bound,nvmax=as.integer(leaps.obj$nvmax),
                 ress=leaps.obj$ress,ir=as.integer(leaps.obj$ir),
                 nbest=as.integer(leaps.obj$nbest),
                 lopt=matrix(as.integer(leaps.obj$lopt),ncol=nbest),
                 il=as.integer(leaps.obj$il),wk=numeric(dimwk),
                 dimwk=as.integer(dimwk),ier=as.integer(0), PACKAGE="leaps")
  rval$dimwk<-rval$wk<-NULL
  rval$xnames<-leaps.obj$xnames
  rval$method<-c("forward",leaps.obj$method)
  rval$force.in<-leaps.obj$force.in
  rval$force.out<-leaps.obj$force.out
  rval$sserr<-leaps.obj$sserr
  rval$intercept<-leaps.obj$intercept
  rval$lindep<-leaps.obj$lindep
  rval$reorder<-leaps.obj$reorder
  rval$nullrss<-leaps.obj$nullrss
  rval$nn<-leaps.obj$nn
  if(nested) {
      rval$nbest<-1
      rval$lopt[]<-rval$vorder[rval$lopt]
      rval$ress<-matrix(rval$rss,ncol=1)
  }
  class(rval)<-"regsubsets"
  if(rval$ier!=0)
      warning(paste("FORWARD returned error code",rval$ier))
  warn.extra(rval)
  rval
}


leaps.seqrep<-function(leaps.obj){
    if (!inherits(leaps.obj,"regsubsets")){
        stop("Not a regsubsets object -- must run leaps.setup")
    }
    nbest<-leaps.obj$nbest
    dimwk<-3*leaps.obj$last
    rval<-.Fortran("seqrep",np=as.integer(leaps.obj$np),
                   nrbar=as.integer(leaps.obj$nrbar),
                   d=leaps.obj$d,rbar=leaps.obj$rbar,
                   thetab=leaps.obj$thetab,
                   first=as.integer(leaps.obj$first),
                   last=as.integer(leaps.obj$last),
                   vorder=as.integer(leaps.obj$vorder),
                   tol=leaps.obj$tol,rss=leaps.obj$rss,
                   bound=leaps.obj$bound,nvmax=as.integer(leaps.obj$nvmax),
                   ress=leaps.obj$ress,ir=as.integer(leaps.obj$ir),
                   nbest=as.integer(leaps.obj$nbest),
                   lopt=matrix(as.integer(leaps.obj$lopt),
                   ncol=nbest),il=as.integer(leaps.obj$il),
                   wk=numeric(dimwk),dimwk=as.integer(dimwk),
                   ier=as.integer(0), PACKAGE="leaps")
  rval$dimwk<-rval$wk<-NULL
  rval$xnames<-leaps.obj$xnames
  rval$method<-c("'sequential replacement'",leaps.obj$method)
  rval$force.in<-leaps.obj$force.in
  rval$force.out<-leaps.obj$force.out
  rval$sserr<-leaps.obj$sserr
  rval$intercept<-leaps.obj$intercept
  rval$lindep<-leaps.obj$lindep
  rval$reorder<-leaps.obj$reorder
  rval$nullrss<-leaps.obj$nullrss
  rval$nn<-leaps.obj$nn
  class(rval)<-"regsubsets"
  if(rval$ier!=0)
      warning(paste("SEQREP returned error code",rval$ier))
  rval
}



print.regsubsets<-function(x,...){
    ll<-x #CMD check
    cat("Subset selection object\n")
    if (!is.null(ll$call)) {
        cat("Call: ")
        print(ll$call)
    }
    cat(ll$np-ll$intercept)
    cat(" Variables ")
    if (ll$intercept) cat(" (and intercept)")
    cat("\n")
    fmat<-cbind(ll$force.in,ll$force.out)
    colnames(fmat)<-c("Forced in","Forced out")
    rownames(fmat)<-ll$xnames
    print(fmat[-1,])
    cat(ll$nbest)
    cat(" subsets of each size up to ")
    cat(ll$nvmax-ll$intercept)
    cat("\n")
    cat("Selection Algorithm: ")
    if (is.null(ll$method)) cat(" not done") else cat(ll$method)
    cat("\n")
    invisible(NULL)
}

print.summary.regsubsets<-function(x,...){
    print(x$obj)
    print(x$outmat)
}

summary.regsubsets<-function(object,all.best=TRUE,matrix=TRUE,matrix.logical=FALSE,
                             df=NULL,...){
    ll<-object #CMD check
    triangle<-function(k) {j<-k-1;1+j*(j+1)/2}
    nmodl<-ll$nbest*ll$nvmax
    if(all.best) nshow<-ll$nbest else nshow<-1
    if (!is.null(df)) n1<-df else n1<-ll$nn-ll$intercept
    outmat<-NULL
    rmat<-NULL
    rnames<-NULL	
    outnames<-NULL
    rsqvec<-NULL
    cpvec<-NULL
    adjr2vec<-NULL
    bicvec<-NULL
    rssvec<-NULL
    sigma2<-ll$sserr/(n1+ll$intercept-ll$last)
    for (i in ll$first:min(ll$last,ll$nvmax)){
        if(!matrix) outmat<-NULL
        for(j in 1:nshow){
            if (ll$ress[i,j]>=1e35) next
            if ((j>1) &
                (all(ll$lopt[triangle(i):(triangle(i+1)-1),j-1]==
                     ll$lopt[triangle(i):(triangle(i+1)-1),j])))
                next
            rline<-rep(FALSE,ll$np)
            rline[ll$lopt[triangle(i):(triangle(i+1)-1),j]]<-TRUE
            outnames<-c(outnames,paste(i-ll$intercept," (",j,")"))
            rnames<-c(rnames,as.character(i-ll$intercept))
            rmat<-rbind(rmat,rline)
            vr<-ll$ress[i,j]/ll$nullrss
            rssvec<-c(rssvec,ll$ress[i,j])
            rsqvec<-c(rsqvec,1-vr)
            adjr2vec<-c(adjr2vec,1-vr*n1/(n1+ll$intercept-i))
            cpvec<-c(cpvec,ll$ress[i,j]/sigma2-(n1+ll$intercept-2*i))
            bicvec<-c(bicvec,(n1+ll$intercept)*log(vr)+i*log(n1+ll$intercept))
        }
    }
    rownames(rmat)<-rnames
    cn<-ll$xnames
    colnames(rmat)<-cn
    reorder<-if (is.null(ll$reorder)) 1:NCOL(rmat) else c(1,1+ll$reorder)
    rmat<-rmat[,order(reorder),drop=FALSE]
    if (matrix){
        if (!matrix.logical)
            outmat<-ifelse(rmat,"*"," ")
        else
            outmat<-rmat
        rownames(outmat)<-outnames
        if (ll$intercept) outmat<-outmat[,-1,drop=FALSE]
    }
    rval<-list(which=rmat,rsq=rsqvec,rss=rssvec,adjr2=adjr2vec,
               cp=cpvec,bic=bicvec,outmat=outmat,obj=ll)
    class(rval)<-"summary.regsubsets"
    rval
}
  

regsubsets<-function(x,...){
  UseMethod("regsubsets",x)
}

regsubsets.default<-function(x,y,weights=rep(1,length(y)),nbest=1,
                             nvmax=8,force.in=NULL,force.out=NULL,
                             intercept=TRUE,
                             method=c("exhaustive","backward","forward","seqrep"),
                             really.big=FALSE,nested=(nbest==1),...)
{
    
    a<-leaps.setup(x,y,wt=weights,nbest=nbest,nvmax=nvmax,
                   force.in=force.in,force.out=force.out,
                   intercept=intercept)
    switch(1+pmatch(method[1],
                    c("exhaustive","backward","forward","seqrep"),
                    nomatch=0),
           stop(paste("Ambiguous or unrecognised method name :",method)),
           leaps.exhaustive(a,really.big=really.big),
           leaps.backward(a,nested),
           leaps.forward(a,nested),
           leaps.seqrep(a))
}



regsubsets.biglm<-function(x,nbest=1,
                             nvmax=8,force.in=NULL,
                             method=c("exhaustive","backward","forward","seqrep"),
                             really.big=FALSE,nested=(nbest==1),...)
{
    
    a<-leaps.from.biglm(x,nbest=nbest,nvmax=nvmax,
                   force.in=force.in)
    switch(1+pmatch(method[1],
                    c("exhaustive","backward","forward","seqrep"),
                    nomatch=0),
           stop(paste("Ambiguous or unrecognised method name :",method)),
           leaps.exhaustive(a,really.big=really.big),
           leaps.backward(a,nested),
           leaps.forward(a,nested),
           leaps.seqrep(a))
}


regsubsets.formula<-function(x,data,weights=NULL,nbest=1,nvmax=8,force.in=NULL,
                             force.out=NULL,intercept=TRUE,
                             method=c("exhaustive","backward","forward","seqrep"),
                             really.big=FALSE,nested=(nbest==1),...){
  formula<-x
  rm(x)
  mm<-match.call()
  mm$formula<-formula
  mm$x<-NULL
  mm$nbest<-mm$nvmax<-mm$force.in<-mm$force.out<-NULL
  mm$intercept<-mm$method<-mm$really.big<-NULL
  mm[[1]]<-as.name("model.frame")
  mm<-eval(mm,sys.frame(sys.parent()))
  x<-model.matrix(terms(formula,data=data),mm)[,-1]
  y<-model.extract(mm,"response")
  wt<-model.extract(mm,"weights")
  if (is.null(wt))
      wt<-rep(1,length(y))
  else
      wt<-weights
  a<-leaps.setup(x,y,wt=wt,nbest=nbest,nvmax=nvmax, force.in=force.in,
                 force.out=force.out, intercept=intercept)
  rval<-switch(1+pmatch(method[1],
                        c("exhaustive","backward","forward","seqrep"),
                        nomatch=0),
               stop(paste("Ambiguous or unrecognised method name :",method)),
               leaps.exhaustive(a,really.big),
               leaps.backward(a,nested),
               leaps.forward(a,nested),
               leaps.seqrep(a))
  rval$call<-sys.call(sys.parent())
  rval
}

  
leaps<-function(x,y,wt=rep(1,NROW(x)),int=TRUE,method=c("Cp","adjr2","r2"),
                nbest=10,names=NULL,df=NROW(x),strictly.compatible=TRUE){
    if (!is.logical(int))
        stop("int should be TRUE or FALSE")
    if (!is.null(names))
        colnames(x)<-names
    method<-method[1]
    if (pmatch(method,c("Cp","adjr2","r2"),nomatch=0)==0)
        stop("Ambiguous or unrecognised method name")
    if (strictly.compatible){
        if (NCOL(x)>31)
            stop("leaps does not allow more than 31 variables; use regsubsets()")
        if (is.null(names))
            colnames(x)<-c(as.character(1:9),LETTERS)[1:NCOL(x)]
    }
    a<-leaps.setup(x,y,wt=wt,nbest=nbest,nvmax=NCOL(x)+int,
                   intercept=int,warn.dep=FALSE)
    if (strictly.compatible & any(a$lindep))
        stop("leaps requires full-rank design matrix; use regsubsets()")
    b<-leaps.exhaustive(a)
    d<-summary(b)
    rval<-list(which=d$which)
    if (int)
        rval$which<-rval$which[,-1,drop=FALSE]
    rval$label<-colnames(d$which)
    rval$size<-as.numeric(rownames(d$which))+int
    if (pmatch(method,c("Cp"),nomatch=0)==1){
        rval$Cp<-d$cp
    }
    if (pmatch(method,c("r2"),nomatch=0)==1){
        rval$r2<-d$rsq
    }
    if (pmatch(method,c("adjr2"),nomatch=0)==1){
        rval$adjr2<-d$adjr2
    }
    rval
}

vcov.regsubsets<-function(object, id, ...){
  betas<-coef(object,id, vcov=TRUE)
  if (length(id)==1)
    attr(betas,"vcov")
  else
    lapply(betas, function(beta) attr(beta,"vcov"))
}

coef.regsubsets<-function(object, id,vcov=FALSE,...){
   s<-summary(object)
   invars<-s$which[id,,drop=FALSE]
   betas<-vector("list",length(id))
   for(i in 1:length(id)){
     thismodel<-which(invars[i,])
     qr<-.Fortran("REORDR", np=as.integer(object$np), nrbar=as.integer(object$nrbar),
                  vorder=as.integer(object$vorder),
                  d=as.double(object$d), rbar=as.double(object$rbar), thetab=as.double(object$thetab),
                  rss=as.double(object$rss), tol=as.double(object$tol), list=as.integer(thismodel),
                  n=as.integer(length(thismodel)),pos1=1L, ier=integer(1))
     beta<-.Fortran("REGCF",np=as.integer(qr$np), nrbar=as.integer(qr$nrbar),
                    d=as.double(qr$d), rbar=as.double(qr$rbar), thetab=as.double(qr$thetab),tol=as.double(qr$tol),
                    beta=numeric(length(thismodel)), nreq=as.integer(length(thismodel)), ier=numeric(1))$beta
     names(beta)<-object$xnames[qr$vorder[1:qr$n]]
     reorder<-order(qr$vorder[1:qr$n])
     beta<-beta[reorder]
     if(vcov){
       p<-length(thismodel)
       R<-diag(qr$np)
       R[row(R)>col(R)]<-qr$rbar
       R<-t(R)
       R<-sqrt(qr$d)*R
       R<-R[1:p,1:p,drop=FALSE]
       R<-chol2inv(R)
       dimnames(R)<-list(object$xnames[qr$vorder[1:p]],object$xnames[qr$vorder[1:p]])
       V<-R*s$rss[id[i]]/(object$nn-p)
       V<-V[reorder,reorder]
       attr(beta,"vcov")<-V
     }
     betas[[i]]<-beta
   }
   if(length(id)==1)
     beta
   else
     betas
}

Try the leaps package in your browser

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

leaps documentation built on Jan. 17, 2020, 1:09 a.m.