[R] Calculation of BIC done by leaps-package

Jan Henckens jan at henckens.de
Sun Dec 26 15:47:31 CET 2010


Hi Folks,

I've got a question concerning the calculation of the Schwarz-Criterion 
(BIC) done by summary.regsubsets() of the leaps-package:

Using regsubsets() to perform subset-selection I receive an regsubsets 
object that can be summarized by summary.regsubsets(). After this 
operation the resulting summary contains a vector of BIC-values 
representing models of size i=1,...,K.

My problem is that I can't reproduce the calculation of these BIC 
values. I already tried to use extractAIC(...,k=log(n)), 
AIC(...,k=log(n)) and manual calculation using the RSS-vector but none 
matches the calculation done by the summary-function. I already checked 
for constants that could be the reason for the differences but i found 
out, that the values vary apart of adding a constant term.


The source code of the leaps-package states the package calculates the 
BIC this way:

bicvec<-c(bicvec,(n1+ll$intercept)*log(vr)+i*log(n1+ll$intercept))

with:

## number of observations - Intercept:
n1<-ll$nn-ll$intercept
## fraction of sum of squared residulas model i
## and sum of squared residuals null model, I
## just can't understand why the vector ll$ress
## is subscripted double
vr<-ll$ress[i,j]/ll$nullrss
## maximum number of variables
i

^^ This seems to match the calculation done by extractAIC but it doesn't!

Maybe anyone can tell me about the reason of the variation of the 
BIC-values?

Best regards,
Jan Henckens



### Minimal Example:
require(leaps)
bridge <- 
read.table("http://www.stat.tamu.edu/~sheather/book/docs/datasets/bridge.txt", 
header=TRUE)
fmla.full <- formula(Time ~ .)
(lm.model <- summary(regsubsets(fmla.full,data=bridge,weights=NULL, 
intercept=TRUE, method="forward")))
lm.model$bic
### The first two models constructed via lm():
extractAIC(lm(Time~Dwgs,data=bridge),k=log(nrow(bridge)))
extractAIC(lm(Time~Dwgs+Case,data=bridge),k=log(nrow(bridge)))

or see

http://www.henckens.de/min_example.R



-- 
jan.henckens | jöllenbecker str. 58 | 33613 bielefeld | germany
tel 0521-5251970



More information about the R-help mailing list