Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Calculating BIC manually for lm object

Tags:

r

lm

I can't quite figure out why I'm having trouble calculating the Bayesian Information Criterion correctly, and was hoping someone could point me in the right direction.

I am doing this because I am trying to calculate the BIC manually (for plm objects, which don't seem to have an established routine associated with them). I took the formula from the Wikipedia page, which gives the formula for the BIC in terms of the Residual Sum of Squares, rather than the Log Likelihood.

y<-rnorm(100)
x<-rnorm(100)
m.test<-lm(y ~ x) 

n<-100 
rss<-sum(m.test$residuals^2) 
k<-3 
bic.mine<-n*log(rss/n)+k*log(n) #formula from wikipedia
bic.stats<-BIC(m.test) #using stats package
abs(bic.mine-bic.stats) #mine is off!

Running the code a bunch of times, I realize that the difference between the BIC I obtain and the BIC obtained from the stats package is constant, so my suspicion is that I'm missing some kind of scaling factor. Is that right? Thanks in advance.

Edit: Thanks for all the comments. I tried to implement the suggestions and post an answer, but I'm still off by a constant. Revised code below.

y<-rnorm(100)
x<-rnorm(100)
m.test<-lm(y ~ x) 

n<-100
res<-m.test$residuals
rss<-sum(res^2) 
k<-3; df<-n-k; w<-rep(1,N) #params, dfs, weights
ll<-0.5 * (sum(log(w)) - n *
             (log(2 * pi) + 1 - log(n) + log(sum(w * res^2))))
ll.stats<-logLik(m.test)
abs(ll.stats-ll)==0 #same, prob is not here

bic.mine<-n*log(rss/n)+k*log(n) #formula from wikipedia
bic.exact<- -2 * ll + log(n) * df #suggestions from comments
bic.stats<-BIC(m.test) #using stats package
abs(bic.mine-bic.stats) #mine was off
abs(bic.exact-bic.stats) #this is still off, though
like image 823
daanoo Avatar asked Feb 01 '16 13:02

daanoo


1 Answers

Thanks to help from the commenters, here's the answer:

y<-rnorm(100)
x<-rnorm(100)
m<-lm(y ~ x) 

To get the BIC or AIC, you first need the associated log likelihood.

Calculating the log likelihood requires a vector of residuals, the number of observations in the data, and a vector of weights (if applicable)

res<-m$residuals
n<-nrow(m$model)    
w<-rep(1,n) #not applicable

ll<-0.5 * (sum(log(w)) - n * (log(2 * pi) + 1 - log(n) + log(sum(w * res^2))))
ll-logLik(m)==0 #TRUE

Calculating the BIC or AIC requires ll, and additionally requires the df associated with the calculation of the log likelihood, which is equal to the original number of parameters being estimated plus 1.

k.original<-length(m$coefficients)
df.ll<-k.original+1 
bic<- -2 * ll + log(n) * df.ll
bic-BIC(m)==0 #TRUE
aic<- -2 * ll + 2 * df.ll
aic-AIC(m)==0 #TRUE
like image 192
daanoo Avatar answered Oct 12 '22 00:10

daanoo