plm()の結果から、対数尤度を取得する。

Introduction
正規誤差モデルにおける最尤法
Conclusion

Introduction

天下のStata様とは違い、Rの各パッケージは各々の開発者の単騎行動によりつくられているため、痒い所に手が届かない場合も少なくない*1。

パネルデータ分析のパッケージに関してもそれは言えて、Rで固定効果モデルをplm()で推定した際、対数尤度やAIC,BICは推定されないのが玉に瑕である。
なので、自分で作ろう、という話。

◆参考URL

正規誤差モデルにおける最尤法

正規誤差の対数尤度関数

参考URLにあるとおり平均、誤差が $(\mu, \sigma^2)$ である正規分布における対数尤度関数は

$ln L ( \mu , \sigma^2; x) = -\frac{n}{2}(ln(2\pi + \sigma^2)) - \frac{1}{2\sigma^2} \sum_{i=1}^n (x_i-\mu)^2$

となる。
OLSやFixed Effects Modelのような正規誤差を仮定する回帰分析における予測値 $\hat{y}_i$ をこの $\mu$ に置き換えてしまえば、
該当の回帰分析における対数尤度が求められる。

OLSの場合

OLSを推定する際の定番であるlm()関数においては、戻り値にlogLik()関数を適用すると対数尤度が得られる
（分析対象としてはpanelrパッケージのWageDataを用いる）

ols1 <- lm(lwage ~ occ + wks,data=WageData)
logLik(ols1)
#'log Lik.' -2459.692 (df=4)

これを上の方法で、求めてみる

resid_to_logLik <- function(ols_res){
  resid <- residuals( ols_res)
  sgm <- summary(ols_res)$sigma
  n <- length(resid)
  llik <- -1 * (n/2) * (log(2*pi) + log( sgm**2)) - (1 / (2*(sgm**2))) * sum( resid ** 2)  
  return(llik)
} #function

resid_to_logLik(ols1)
#-2459.692

なるほど一致している。

plmの場合

固定効果モデルも内実はOLSなので、基本的にOLSと同じ方法をとればよい。
うえの、対数尤度関数を実際に計算されるためには、予測値 $\hat{y}_i$ （と実測値の差）と、分散推定値 $\sigma^2$ を求める必要がある。

予測値と実測値の乖離は残差として取得が可能である。
分散推定値の $\sigma^2$ については、一般によく知られているように、残差平方和／自由度で求めることができる。
（参考：統計学の基礎(12.18)）

OLSではsummary()$sigmaによって誤差分散（の平方根）が取得できるが、それもdevianceを自由度で引いたもので求めている（→参考URL③参照）。
したがって、以下のように実装すればplmでも対数尤度を求められる。

#関数定義
plm_logLik <- function( plm_res){
  df_v <- plm_res$df.residual #自由度取得
  rss_v <- sum(plm_res$residuals ** 2) #残差平方和
  est_sgm2 <- rss_v / df_v #分散推定値
  n <- length( plm_res$residuals)
  llik <- -1 * (n/2) * (log(2*pi) + log( est_sgm2))
 - (1 / (2*(est_sgm2))) * sum( plm_res$residuals ** 2)  
  return(list( logLik = llik, df= df_v))
}

#実行
plm1 <- plm( lwage~ occ + wks , model = "within" ,
 data = WageData , index=c("id", "t"))

plm_logLik(plm1)
#$logLik
#[1] 9.912123