固定効果ロジットの尤度比検定の関数を書いた

二項ロジット(ロジスティック回帰)でモデルの比較を行うとき、あるモデルが別のモデルのネストになってるならば尤度比検定を行うのが一つの方法だと思うんですが、これをRでやろうと思うとlrtestというパッケージですごく簡単にできます。

ただし、固定効果ロジットを使うとglmで個人のダミー含める推定だとバイアスがかかるので、clogitやbifeというパッケージを使います。でもlrtestはこれらのパッケージに対応してないっぽい。

固定効果ロジットについてはこちらの記事
keita43a.hatenablog.com

そこで、関数を書いてみた。
bifeのバイアス補正を使って固定効果を入れたロジスティック回帰を推定したかったので、bifeのオブジェクトを引数にとる関数である。

# Define LRtest function for bife
                lrtest.bife = function(obj1,obj2){
                   loglik1 = obj1$logl_info$loglik   # Loglikelihood of object 1 (simpler model)
                   loglik2 = obj2$logl_info$loglik  # loglikelihood of object 2 (complicated model)
                   lr = -2*(loglik1 - loglik2)             # Likelihood ratio statistics
                   n1 = length(obj1$par$beta) + length(obj1$par$alpha)   # number of parameter of object 1
                   n2 = length(obj2$par$beta) + length(obj2$par$alpha)  # number of parameter of object 2
                   df = n2-n1 # Degree of freedom
                   pval = pchisq(lr,df,lower.tail = FALSE)
                   star = ifelse(pval<=0.001,"***",
                                 ifelse(pval<=0.01,"**",
                                        ifelse(pval<=0.05,"*","")))
                   return(list(num.coef = c(n1,n2),
                               loglik = c(loglik1,loglik2),
                               lik.ratio = lr, 
                               df.chi = df,
                               pval = pval,
                               star = star))
                }

num.coefがそれぞれのモデルで推定された係数の数。loglikが対数尤度関数の値、lik.ratioが対数尤度比(統計量), df.chiがカイ二乗分布の自由度(=係数の数の差), pvalがp値, starはp値に対応して0.05, 0.01, 0.001の順に*が増える。

注意が必要なのは、関数の引数のうち、１つ目(obj1)はシンプルな方のモデル(帰無仮説)、２つ目が複雑な(変数の多い)モデルである。

bifeパッケージに入っているPSIDというデータで試してみる。
女性の労働参加についてのパネルデータで、9年間の年ごとの値が記録されている。
LFPが0か1をとり、女性の労働参加を表す。INCHは夫の所得、AGEは年齢、KID1は0~2歳の子供の数、KID2は3~5歳の子供の数を表す。
ここでは、仮にモデル1とモデル2の差を、KID2を含めるかどうかとする。

# Package load
library(bife)

# Load data
df1 = bife::psid

# Model 1: 
model1 = bife(LFP ~ INCH + AGE + KID1|ID,
              data = df1,
              model = "logit",
              bias_corr = "ana")
summary(model1)

# Model 2: add KID2
model2 = bife(LFP ~ INCH + AGE + KID1 + KID2|ID, 
              data = df1,
              model = "logit",
              bias_corr = "ana")
summary(model2)


lrtest.bife(model1,model2)

結果は以下のようになる。

> lrtest.bife(model1,model2)
$num.coef
[1] 667 668

$loglik
[1] -3072.623 -3044.168

$lik.ratio
[1] 56.91012

$df.chi
[1] 1

$pval
[1] 4.56193e-14

$star
[1] "***"