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
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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With