# R/leaps.R In leaps: Regression Subset Selection

#### Documented in coef.regsubsetsleapsleaps.backwardleaps.exhaustiveleaps.forwardleaps.seqrepleaps.setupprint.regsubsetsprint.summary.regsubsetsregsubsetsregsubsets.biglmregsubsets.defaultregsubsets.formulasummary.regsubsetsvcov.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,
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,
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),
if (initr\$ier!=0)
warning(paste("INITR returned error code",initr\$ier))
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),
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,
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,
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),
if (initr\$ier!=0)
warning(paste("INITR returned error code",initr\$ier))
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),
class(rval)<-"regsubsets"
invisible(rval)
}

warn.extra<-function(obj){
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,
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\$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,
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\$nn<-leaps.obj\$nn
class(rval)<-"regsubsets"
if(nested) {
rval\$nbest<-1
rval\$lopt[]<-rval\$vorder[rval\$lopt]
}
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),
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\$nn<-leaps.obj\$nn
if(nested) {
rval\$nbest<-1
rval\$lopt[]<-rval\$vorder[rval\$lopt]
}
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),
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\$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
bicvec<-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)
rsqvec<-c(rsqvec,1-vr)
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]
}
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
}

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]
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
}
}
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),
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<-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.