固定効果ロジットの尤度比検定の関数を書いた
二項ロジット(ロジスティック回帰)でモデルの比較を行うとき、あるモデルが別のモデルのネストになってるならば尤度比検定を行うのが一つの方法だと思うんですが、これを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の順に*が増える。
注意が必要なのは、関数の引数のうち、1つ目(obj1)はシンプルな方のモデル(帰無仮説)、2つ目が複雑な(変数の多い)モデルである。
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] "***"
つまり、モデル1が正しいという帰無仮説のもとで、KID2を含めたモデル2との比較検定を行った結果、帰無仮説は棄却されたということになります。
ポイントは、model1とmodel2は両方とも個人固定効果(ID)を含んだモデルという点です。